File: cholesky.c

package info (click to toggle)
gsl-doc 2.3-1
  • links: PTS
  • area: non-free
  • in suites: buster
  • size: 27,748 kB
  • ctags: 15,177
  • sloc: ansic: 235,014; sh: 11,585; makefile: 925
file content (324 lines) | stat: -rw-r--r-- 8,781 bytes parent folder | download | duplicates (4)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
/* multilarge_nlinear/cholesky.c
 * 
 * Copyright (C) 2016 Patrick Alken
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/*
 * This module calculates the solution of the normal equations least squares
 * system:
 *
 * [ J^T J + mu D^T D ] p = -J^T f
 *
 * using the modified Cholesky decomposition.
 */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_multilarge_nlinear.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_permutation.h>

#include "common.c"

typedef struct
{
  gsl_matrix *JTJ;           /* J^T J */
  gsl_matrix *work_JTJ;      /* copy of J^T J */
  gsl_vector *rhs;           /* -J^T f, size p */
  gsl_permutation *perm;     /* permutation matrix for modified Cholesky */
  gsl_vector *work3p;        /* workspace, size 3*p */
  gsl_vector *workn;         /* workspace, size n */
  double mu;                 /* current regularization parameter */
} cholesky_state_t;

static void *cholesky_alloc (const size_t n, const size_t p);
static int cholesky_init(const void * vtrust_state, void * vstate);
static int cholesky_presolve(const double mu, const void * vtrust_state, void * vstate);
static int cholesky_solve(const gsl_vector * g, gsl_vector *x,
                          const  void * vtrust_state, void *vstate);
static int cholesky_rcond(double * rcond, const gsl_matrix * JTJ, void * vstate);
static int cholesky_covar(const gsl_matrix * JTJ, gsl_matrix * covar, void * vstate);
static int cholesky_solve_rhs(const gsl_vector * b, gsl_vector *x, cholesky_state_t *state);
static int cholesky_regularize(const double mu, const gsl_vector * diag, gsl_matrix * A,
                               cholesky_state_t * state);

static void *
cholesky_alloc (const size_t n, const size_t p)
{
  cholesky_state_t *state;

  state = calloc(1, sizeof(cholesky_state_t));
  if (state == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate cholesky state", GSL_ENOMEM);
    }

  state->JTJ = gsl_matrix_alloc(p, p);
  if (state->JTJ == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for JTJ", GSL_ENOMEM);
    }

  state->work_JTJ = gsl_matrix_alloc(p, p);
  if (state->work_JTJ == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for JTJ workspace",
                      GSL_ENOMEM);
    }

  state->rhs = gsl_vector_alloc(p);
  if (state->rhs == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for rhs", GSL_ENOMEM);
    }

  state->perm = gsl_permutation_alloc(p);
  if (state->perm == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for perm", GSL_ENOMEM);
    }

  state->work3p = gsl_vector_alloc(3 * p);
  if (state->work3p == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for work3p", GSL_ENOMEM);
    }

  state->workn = gsl_vector_alloc(n);
  if (state->workn == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for workn", GSL_ENOMEM);
    }

  state->mu = -1.0;

  return state;
}

static void
cholesky_free(void *vstate)
{
  cholesky_state_t *state = (cholesky_state_t *) vstate;

  if (state->JTJ)
    gsl_matrix_free(state->JTJ);

  if (state->work_JTJ)
    gsl_matrix_free(state->work_JTJ);

  if (state->rhs)
    gsl_vector_free(state->rhs);

  if (state->perm)
    gsl_permutation_free(state->perm);

  if (state->work3p)
    gsl_vector_free(state->work3p);

  if (state->workn)
    gsl_vector_free(state->workn);

  free(state);
}

static int
cholesky_init(const void * vtrust_state, void * vstate)
{
  const gsl_multilarge_nlinear_trust_state *trust_state =
    (const gsl_multilarge_nlinear_trust_state *) vtrust_state;
  cholesky_state_t *state = (cholesky_state_t *) vstate;

  /* store J^T J normal equations matrix */
  gsl_matrix_tricpy('L', 1, state->JTJ, trust_state->JTJ);

  return GSL_SUCCESS;
}

/*
cholesky_presolve()
  Compute the modified Cholesky decomposition of J^T J + mu D^T D.
Modified Cholesky is used in case mu = 0 and there are rounding
errors in forming J^T J which could lead to an indefinite matrix.

Inputs: mu     - LM parameter
        vstate - workspace

Notes:
1) On output, state->work_JTJ contains the Cholesky decomposition of
J^T J + mu D^T D
*/

static int
cholesky_presolve(const double mu, const void * vtrust_state, void * vstate)
{
  const gsl_multilarge_nlinear_trust_state *trust_state =
    (const gsl_multilarge_nlinear_trust_state *) vtrust_state;
  cholesky_state_t *state = (cholesky_state_t *) vstate;
  gsl_matrix *JTJ = state->work_JTJ;
  const gsl_vector *diag = trust_state->diag;
  int status;

  /* copy lower triangle of A to workspace */
  gsl_matrix_tricpy('L', 1, JTJ, state->JTJ);

  /* augment normal equations: A -> A + mu D^T D */
  status = cholesky_regularize(mu, diag, JTJ, state);
  if (status)
    return status;

  /* compute modified Cholesky decomposition */
  status = gsl_linalg_mcholesky_decomp(JTJ, state->perm, NULL);
  if (status)
    return status;

  state->mu = mu;

  return GSL_SUCCESS;
}

/*
cholesky_solve()
  Compute (J^T J + mu D^T D) x = -g

where g = J^T f

Inputs: g      - right hand side vector g, size p
        x      - (output) solution vector
        vstate - cholesky workspace
*/

static int
cholesky_solve(const gsl_vector * g, gsl_vector *x,
               const void * vtrust_state, void *vstate)
{
  cholesky_state_t *state = (cholesky_state_t *) vstate;
  int status;

  status = cholesky_solve_rhs(g, x, state);
  if (status)
    return status;

  /* reverse direction to go downhill */
  gsl_vector_scale(x, -1.0);

  (void) vtrust_state;

  return GSL_SUCCESS;
}

static int
cholesky_rcond(double * rcond, const gsl_matrix * JTJ, void * vstate)
{
  int status;
  cholesky_state_t *state = (cholesky_state_t *) vstate;
  double rcond_JTJ;

  /* its possible the current Cholesky decomposition is from the previous
   * iteration so do a new one to be sure we use the right Jacobian */

  /* copy lower triangle of JTJ to workspace */
  gsl_matrix_tricpy('L', 1, state->work_JTJ, JTJ);

  /* compute modified Cholesky decomposition */
  status = gsl_linalg_mcholesky_decomp(state->work_JTJ, state->perm, NULL);
  if (status)
    return status;

  status = gsl_linalg_mcholesky_rcond(state->work_JTJ, state->perm, &rcond_JTJ, state->work3p);
  if (status == GSL_SUCCESS)
    *rcond = sqrt(rcond_JTJ);

  return status;
}

static int
cholesky_covar(const gsl_matrix * JTJ, gsl_matrix * covar, void * vstate)
{
  int status;
  cholesky_state_t *state = (cholesky_state_t *) vstate;

  /* its possible the current Cholesky decomposition is from the previous
   * iteration so do a new one to be sure we use the right Jacobian */

  /* copy lower triangle of JTJ to workspace */
  gsl_matrix_tricpy('L', 1, state->work_JTJ, JTJ);

  /* compute modified Cholesky decomposition */
  status = gsl_linalg_mcholesky_decomp(state->work_JTJ, state->perm, NULL);
  if (status)
    return status;

  status = gsl_linalg_mcholesky_invert(state->work_JTJ, state->perm, covar);
  if (status)
    return status;

  return GSL_SUCCESS;
}

/* solve: (J^T J + mu D^T D) x = b */
static int
cholesky_solve_rhs(const gsl_vector * b, gsl_vector *x, cholesky_state_t *state)
{
  int status;
  gsl_matrix *JTJ = state->work_JTJ;

  status = gsl_linalg_mcholesky_solve(JTJ, state->perm, b, x);
  if (status)
    return status;

  return GSL_SUCCESS;
}

/* A <- A + mu D^T D */
static int
cholesky_regularize(const double mu, const gsl_vector * diag, gsl_matrix * A,
                    cholesky_state_t * state)
{
  (void) state;

  if (mu != 0.0)
    {
      size_t i;

      for (i = 0; i < diag->size; ++i)
        {
          double di = gsl_vector_get(diag, i);
          double *Aii = gsl_matrix_ptr(A, i, i);
          *Aii += mu * di * di;
        }
    }

  return GSL_SUCCESS;
}

static const gsl_multilarge_nlinear_solver cholesky_type =
{
  "cholesky",
  cholesky_alloc,
  cholesky_init,
  cholesky_presolve,
  cholesky_solve,
  cholesky_rcond,
  cholesky_covar,
  cholesky_free
};

const gsl_multilarge_nlinear_solver *gsl_multilarge_nlinear_solver_cholesky = &cholesky_type;