1 /* multifit/multireg.c
2 *
3 * Copyright (C) 2015 Patrick Alken
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 /*
21 * References:
22 *
23 * [1] P. C. Hansen & D. P. O'Leary, "The use of the L-curve in
24 * the regularization of discrete ill-posed problems", SIAM J. Sci.
25 * Comput. 14 (1993), pp. 1487-1503.
26 *
27 * [2] P. C. Hansen, "Discrete Inverse Problems: Insight and Algorithms,"
28 * SIAM Press, 2010.
29 */
30
31 #include <config.h>
32 #include <gsl/gsl_errno.h>
33 #include <gsl/gsl_multifit.h>
34 #include <gsl/gsl_blas.h>
35 #include <gsl/gsl_vector.h>
36 #include <gsl/gsl_matrix.h>
37 #include <gsl/gsl_linalg.h>
38
39 #include "linear_common.c"
40
41 int
gsl_multifit_linear_solve(const double lambda,const gsl_matrix * X,const gsl_vector * y,gsl_vector * c,double * rnorm,double * snorm,gsl_multifit_linear_workspace * work)42 gsl_multifit_linear_solve (const double lambda,
43 const gsl_matrix * X,
44 const gsl_vector * y,
45 gsl_vector * c,
46 double *rnorm,
47 double *snorm,
48 gsl_multifit_linear_workspace * work)
49 {
50 size_t rank;
51 int status;
52
53 status = multifit_linear_solve(X, y, GSL_DBL_EPSILON, lambda, &rank, c,
54 rnorm, snorm, work);
55
56 return status;
57 } /* gsl_multifit_linear_solve() */
58
59 /*
60 gsl_multifit_linear_applyW()
61 Apply weight matrix to (X,y) LS system
62
63 Inputs: X - least squares matrix n-by-p
64 w - weight vector n-by-1 or NULL for W = I
65 y - right hand side n-by-1
66 WX - (output) sqrt(W) X, n-by-p
67 Wy - (output) sqrt(W) y, n-by-1
68
69 Notes:
70 1) If w = NULL, on output WX = X and Wy = y
71 2) It is allowed for WX = X and Wy = y for in-place transform
72 */
73
74 int
gsl_multifit_linear_applyW(const gsl_matrix * X,const gsl_vector * w,const gsl_vector * y,gsl_matrix * WX,gsl_vector * Wy)75 gsl_multifit_linear_applyW(const gsl_matrix * X,
76 const gsl_vector * w,
77 const gsl_vector * y,
78 gsl_matrix * WX,
79 gsl_vector * Wy)
80 {
81 const size_t n = X->size1;
82 const size_t p = X->size2;
83
84 if (n != y->size)
85 {
86 GSL_ERROR("y vector does not match X", GSL_EBADLEN);
87 }
88 else if (w != NULL && n != w->size)
89 {
90 GSL_ERROR("weight vector does not match X", GSL_EBADLEN);
91 }
92 else if (n != WX->size1 || p != WX->size2)
93 {
94 GSL_ERROR("WX matrix dimensions do not match X", GSL_EBADLEN);
95 }
96 else if (n != Wy->size)
97 {
98 GSL_ERROR("Wy vector must be length n", GSL_EBADLEN);
99 }
100 else
101 {
102 size_t i;
103
104 /* copy WX = X; Wy = y if distinct pointers */
105 if (WX != X)
106 gsl_matrix_memcpy(WX, X);
107 if (Wy != y)
108 gsl_vector_memcpy(Wy, y);
109
110 if (w != NULL)
111 {
112 /* construct WX = sqrt(W) X and Wy = sqrt(W) y */
113 for (i = 0; i < n; ++i)
114 {
115 double wi = gsl_vector_get(w, i);
116 double swi;
117 gsl_vector_view row = gsl_matrix_row(WX, i);
118 double *yi = gsl_vector_ptr(Wy, i);
119
120 if (wi < 0.0)
121 wi = 0.0;
122
123 swi = sqrt(wi);
124 gsl_vector_scale(&row.vector, swi);
125 *yi *= swi;
126 }
127 }
128
129 return GSL_SUCCESS;
130 }
131 }
132
133 /*
134 gsl_multifit_linear_wstdform1()
135 Using regularization matrix
136 L = diag(l_1,l_2,...,l_p), transform to Tikhonov standard form:
137
138 X~ = sqrt(W) X L^{-1}
139 y~ = sqrt(W) y
140 c~ = L c
141
142 Inputs: L - Tikhonov matrix as a vector of diagonal elements p-by-1;
143 or NULL for L = I
144 X - least squares matrix n-by-p
145 y - right hand side vector n-by-1
146 w - weight vector n-by-1; or NULL for W = I
147 Xs - least squares matrix in standard form X~ n-by-p
148 ys - right hand side vector in standard form y~ n-by-1
149 work - workspace
150
151 Return: success/error
152
153 Notes:
154 1) It is allowed for X = Xs and y = ys
155 */
156
157 int
gsl_multifit_linear_wstdform1(const gsl_vector * L,const gsl_matrix * X,const gsl_vector * w,const gsl_vector * y,gsl_matrix * Xs,gsl_vector * ys,gsl_multifit_linear_workspace * work)158 gsl_multifit_linear_wstdform1 (const gsl_vector * L,
159 const gsl_matrix * X,
160 const gsl_vector * w,
161 const gsl_vector * y,
162 gsl_matrix * Xs,
163 gsl_vector * ys,
164 gsl_multifit_linear_workspace * work)
165 {
166 const size_t n = X->size1;
167 const size_t p = X->size2;
168
169 if (n > work->nmax || p > work->pmax)
170 {
171 GSL_ERROR("observation matrix larger than workspace", GSL_EBADLEN);
172 }
173 else if (L != NULL && p != L->size)
174 {
175 GSL_ERROR("L vector does not match X", GSL_EBADLEN);
176 }
177 else if (n != y->size)
178 {
179 GSL_ERROR("y vector does not match X", GSL_EBADLEN);
180 }
181 else if (w != NULL && n != w->size)
182 {
183 GSL_ERROR("weight vector does not match X", GSL_EBADLEN);
184 }
185 else if (n != Xs->size1 || p != Xs->size2)
186 {
187 GSL_ERROR("Xs matrix dimensions do not match X", GSL_EBADLEN);
188 }
189 else if (n != ys->size)
190 {
191 GSL_ERROR("ys vector must be length n", GSL_EBADLEN);
192 }
193 else
194 {
195 int status = GSL_SUCCESS;
196
197 /* compute Xs = sqrt(W) X and ys = sqrt(W) y */
198 status = gsl_multifit_linear_applyW(X, w, y, Xs, ys);
199 if (status)
200 return status;
201
202 if (L != NULL)
203 {
204 size_t j;
205
206 /* construct X~ = sqrt(W) X * L^{-1} matrix */
207 for (j = 0; j < p; ++j)
208 {
209 gsl_vector_view Xj = gsl_matrix_column(Xs, j);
210 double lj = gsl_vector_get(L, j);
211
212 if (lj == 0.0)
213 {
214 GSL_ERROR("L matrix is singular", GSL_EDOM);
215 }
216
217 gsl_vector_scale(&Xj.vector, 1.0 / lj);
218 }
219 }
220
221 return status;
222 }
223 }
224
225 /*
226 gsl_multifit_linear_stdform1()
227 Using regularization matrix L = diag(l_1,l_2,...,l_p),
228 and W = I, transform to Tikhonov standard form:
229
230 X~ = X L^{-1}
231 y~ = y
232 c~ = L c
233
234 Inputs: L - Tikhonov matrix as a vector of diagonal elements p-by-1
235 X - least squares matrix n-by-p
236 y - right hand side vector n-by-1
237 Xs - least squares matrix in standard form X~ n-by-p
238 ys - right hand side vector in standard form y~ n-by-1
239 work - workspace
240
241 Return: success/error
242
243 Notes:
244 1) It is allowed for X = Xs
245 */
246
247 int
gsl_multifit_linear_stdform1(const gsl_vector * L,const gsl_matrix * X,const gsl_vector * y,gsl_matrix * Xs,gsl_vector * ys,gsl_multifit_linear_workspace * work)248 gsl_multifit_linear_stdform1 (const gsl_vector * L,
249 const gsl_matrix * X,
250 const gsl_vector * y,
251 gsl_matrix * Xs,
252 gsl_vector * ys,
253 gsl_multifit_linear_workspace * work)
254 {
255 int status;
256
257 status = gsl_multifit_linear_wstdform1(L, X, NULL, y, Xs, ys, work);
258
259 return status;
260 }
261
262 int
gsl_multifit_linear_L_decomp(gsl_matrix * L,gsl_vector * tau)263 gsl_multifit_linear_L_decomp (gsl_matrix * L, gsl_vector * tau)
264 {
265 const size_t m = L->size1;
266 const size_t p = L->size2;
267 int status;
268
269 if (tau->size != GSL_MIN(m, p))
270 {
271 GSL_ERROR("tau vector must be min(m,p)", GSL_EBADLEN);
272 }
273 else if (m >= p)
274 {
275 /* square or tall L matrix */
276 status = gsl_linalg_QR_decomp(L, tau);
277 return status;
278 }
279 else
280 {
281 /* more columns than rows, compute qr(L^T) */
282 gsl_matrix_view LTQR = gsl_matrix_view_array(L->data, p, m);
283 gsl_matrix *LT = gsl_matrix_alloc(p, m);
284
285 /* XXX: use temporary storage due to difficulties in transforming
286 * a rectangular matrix in-place */
287 gsl_matrix_transpose_memcpy(LT, L);
288 gsl_matrix_memcpy(<QR.matrix, LT);
289 gsl_matrix_free(LT);
290
291 status = gsl_linalg_QR_decomp(<QR.matrix, tau);
292
293 return status;
294 }
295 }
296
297 /*
298 gsl_multifit_linear_wstdform2()
299 Using regularization matrix L which is m-by-p, transform to Tikhonov
300 standard form. This routine is separated into two cases:
301
302 Case 1: m >= p, here we can use the QR decomposition of L = QR, and note
303 that ||L c|| = ||R c|| where R is p-by-p. Therefore,
304
305 X~ = X R^{-1} is n-by-p
306 y~ = y is n-by-1
307 c~ is p-by-1
308 M is not used
309
310 Case 2: m < p
311
312 X~ is (n - p + m)-by-m
313 y~ is (n - p + m)-by-1
314 c~ is m-by-1
315 M is n-by-p (workspace)
316
317 Inputs: LQR - output from gsl_multifit_linear_L_decomp()
318 Ltau - output from gsl_multifit_linear_L_decomp()
319 X - least squares matrix n-by-p
320 w - weight vector n-by-1; or NULL for W = I
321 y - right hand side vector n-by-1
322 Xs - (output) least squares matrix in standard form
323 case 1: n-by-p
324 case 2: (n - p + m)-by-m
325 ys - (output) right hand side vector in standard form
326 case 1: n-by-1
327 case 2: (n - p + m)-by-1
328 M - (output) workspace matrix needed to reconstruct solution vector
329 case 1: not used
330 case 2: n-by-p
331 work - workspace
332
333 Return: success/error
334
335 Notes:
336 1) If m >= p, on output:
337 Xs = X R^{-1}
338 ys = y
339
340 2) If m < p, on output:
341 M(:,1:pm) contains QR decomposition of A * K_o, needed to reconstruct
342 solution vector, where pm = p - m; M(:,p) contains Householder scalars
343 */
344
345 int
gsl_multifit_linear_wstdform2(const gsl_matrix * LQR,const gsl_vector * Ltau,const gsl_matrix * X,const gsl_vector * w,const gsl_vector * y,gsl_matrix * Xs,gsl_vector * ys,gsl_matrix * M,gsl_multifit_linear_workspace * work)346 gsl_multifit_linear_wstdform2 (const gsl_matrix * LQR,
347 const gsl_vector * Ltau,
348 const gsl_matrix * X,
349 const gsl_vector * w,
350 const gsl_vector * y,
351 gsl_matrix * Xs,
352 gsl_vector * ys,
353 gsl_matrix * M,
354 gsl_multifit_linear_workspace * work)
355 {
356 const size_t m = LQR->size1;
357 const size_t n = X->size1;
358 const size_t p = X->size2;
359
360 if (n > work->nmax || p > work->pmax)
361 {
362 GSL_ERROR("observation matrix larger than workspace", GSL_EBADLEN);
363 }
364 else if (p != LQR->size2)
365 {
366 GSL_ERROR("LQR and X matrices have different numbers of columns", GSL_EBADLEN);
367 }
368 else if (n != y->size)
369 {
370 GSL_ERROR("y vector does not match X", GSL_EBADLEN);
371 }
372 else if (w != NULL && n != w->size)
373 {
374 GSL_ERROR("weights vector must be length n", GSL_EBADLEN);
375 }
376 else if (m >= p) /* square or tall L matrix */
377 {
378 /* the sizes of Xs and ys depend on whether m >= p or m < p */
379 if (n != Xs->size1 || p != Xs->size2)
380 {
381 GSL_ERROR("Xs matrix must be n-by-p", GSL_EBADLEN);
382 }
383 else if (n != ys->size)
384 {
385 GSL_ERROR("ys vector must have length n", GSL_EBADLEN);
386 }
387 else
388 {
389 int status;
390 size_t i;
391 gsl_matrix_const_view R = gsl_matrix_const_submatrix(LQR, 0, 0, p, p);
392
393 /* compute Xs = sqrt(W) X and ys = sqrt(W) y */
394 status = gsl_multifit_linear_applyW(X, w, y, Xs, ys);
395 if (status)
396 return status;
397
398 /* compute X~ = X R^{-1} using QR decomposition of L */
399 for (i = 0; i < n; ++i)
400 {
401 gsl_vector_view v = gsl_matrix_row(Xs, i);
402
403 /* solve: R^T y = X_i */
404 gsl_blas_dtrsv(CblasUpper, CblasTrans, CblasNonUnit, &R.matrix, &v.vector);
405 }
406
407 return GSL_SUCCESS;
408 }
409 }
410 else /* L matrix with m < p */
411 {
412 const size_t pm = p - m;
413 const size_t npm = n - pm;
414
415 /*
416 * This code closely follows section 2.6.1 of Hansen's
417 * "Regularization Tools" manual
418 */
419
420 if (npm != Xs->size1 || m != Xs->size2)
421 {
422 GSL_ERROR("Xs matrix must be (n-p+m)-by-m", GSL_EBADLEN);
423 }
424 else if (npm != ys->size)
425 {
426 GSL_ERROR("ys vector must be of length (n-p+m)", GSL_EBADLEN);
427 }
428 else if (n != M->size1 || p != M->size2)
429 {
430 GSL_ERROR("M matrix must be n-by-p", GSL_EBADLEN);
431 }
432 else
433 {
434 int status;
435 gsl_matrix_view A = gsl_matrix_submatrix(work->A, 0, 0, n, p);
436 gsl_vector_view b = gsl_vector_subvector(work->t, 0, n);
437
438 gsl_matrix_view LTQR = gsl_matrix_view_array(LQR->data, p, m); /* qr(L^T) */
439 gsl_matrix_view Rp = gsl_matrix_view_array(LQR->data, m, m); /* R factor of L^T */
440 gsl_vector_const_view LTtau = gsl_vector_const_subvector(Ltau, 0, m);
441
442 /*
443 * M(:,1:p-m) will hold QR decomposition of A K_o; M(:,p) will hold
444 * Householder scalars
445 */
446 gsl_matrix_view MQR = gsl_matrix_submatrix(M, 0, 0, n, pm);
447 gsl_vector_view Mtau = gsl_matrix_subcolumn(M, p - 1, 0, GSL_MIN(n, pm));
448
449 gsl_matrix_view AKo, AKp, HqTAKp;
450 gsl_vector_view v;
451 size_t i;
452
453 /* compute A = sqrt(W) X and b = sqrt(W) y */
454 status = gsl_multifit_linear_applyW(X, w, y, &A.matrix, &b.vector);
455 if (status)
456 return status;
457
458 /* compute: A <- A K = [ A K_p ; A K_o ] */
459 gsl_linalg_QR_matQ(<QR.matrix, <tau.vector, &A.matrix);
460 AKp = gsl_matrix_submatrix(&A.matrix, 0, 0, n, m);
461 AKo = gsl_matrix_submatrix(&A.matrix, 0, m, n, pm);
462
463 /* compute QR decomposition [H,T] = qr(A * K_o) and store in M */
464 gsl_matrix_memcpy(&MQR.matrix, &AKo.matrix);
465 gsl_linalg_QR_decomp(&MQR.matrix, &Mtau.vector);
466
467 /* AKp currently contains A K_p; apply H^T from the left to get H^T A K_p */
468 gsl_linalg_QR_QTmat(&MQR.matrix, &Mtau.vector, &AKp.matrix);
469
470 /* the last npm rows correspond to H_q^T A K_p */
471 HqTAKp = gsl_matrix_submatrix(&AKp.matrix, pm, 0, npm, m);
472
473 /* solve: Xs R_p^T = H_q^T A K_p for Xs */
474 gsl_matrix_memcpy(Xs, &HqTAKp.matrix);
475 for (i = 0; i < npm; ++i)
476 {
477 gsl_vector_view x = gsl_matrix_row(Xs, i);
478 gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, &Rp.matrix, &x.vector);
479 }
480
481 /*
482 * compute: ys = H_q^T b; this is equivalent to computing
483 * the last q elements of H^T b (q = npm)
484 */
485 v = gsl_vector_subvector(&b.vector, pm, npm);
486 gsl_linalg_QR_QTvec(&MQR.matrix, &Mtau.vector, &b.vector);
487 gsl_vector_memcpy(ys, &v.vector);
488
489 return GSL_SUCCESS;
490 }
491 }
492 }
493
494 int
gsl_multifit_linear_stdform2(const gsl_matrix * LQR,const gsl_vector * Ltau,const gsl_matrix * X,const gsl_vector * y,gsl_matrix * Xs,gsl_vector * ys,gsl_matrix * M,gsl_multifit_linear_workspace * work)495 gsl_multifit_linear_stdform2 (const gsl_matrix * LQR,
496 const gsl_vector * Ltau,
497 const gsl_matrix * X,
498 const gsl_vector * y,
499 gsl_matrix * Xs,
500 gsl_vector * ys,
501 gsl_matrix * M,
502 gsl_multifit_linear_workspace * work)
503 {
504 int status;
505
506 status = gsl_multifit_linear_wstdform2(LQR, Ltau, X, NULL, y, Xs, ys, M, work);
507
508 return status;
509 }
510
511 /*
512 gsl_multifit_linear_genform1()
513 Backtransform regularized solution vector using matrix
514 L = diag(L)
515 */
516
517 int
gsl_multifit_linear_genform1(const gsl_vector * L,const gsl_vector * cs,gsl_vector * c,gsl_multifit_linear_workspace * work)518 gsl_multifit_linear_genform1 (const gsl_vector * L,
519 const gsl_vector * cs,
520 gsl_vector * c,
521 gsl_multifit_linear_workspace * work)
522 {
523 if (L->size > work->pmax)
524 {
525 GSL_ERROR("L vector does not match workspace", GSL_EBADLEN);
526 }
527 else if (L->size != cs->size)
528 {
529 GSL_ERROR("cs vector does not match L", GSL_EBADLEN);
530 }
531 else if (L->size != c->size)
532 {
533 GSL_ERROR("c vector does not match L", GSL_EBADLEN);
534 }
535 else
536 {
537 /* compute true solution vector c = L^{-1} c~ */
538 gsl_vector_memcpy(c, cs);
539 gsl_vector_div(c, L);
540
541 return GSL_SUCCESS;
542 }
543 }
544
545 /*
546 gsl_multifit_linear_wgenform2()
547 Backtransform regularized solution vector in standard form to recover
548 original vector
549
550 Inputs: LQR - output from gsl_multifit_linear_L_decomp()
551 Ltau - output from gsl_multifit_linear_L_decomp()
552 X - original least squares matrix n-by-p
553 w - original weight vector n-by-1 or NULL for W = I
554 y - original rhs vector n-by-1
555 cs - standard form solution vector
556 c - (output) original solution vector p-by-1
557 M - matrix computed by gsl_multifit_linear_wstdform2()
558 work - workspace
559 */
560
561 int
gsl_multifit_linear_wgenform2(const gsl_matrix * LQR,const gsl_vector * Ltau,const gsl_matrix * X,const gsl_vector * w,const gsl_vector * y,const gsl_vector * cs,const gsl_matrix * M,gsl_vector * c,gsl_multifit_linear_workspace * work)562 gsl_multifit_linear_wgenform2 (const gsl_matrix * LQR,
563 const gsl_vector * Ltau,
564 const gsl_matrix * X,
565 const gsl_vector * w,
566 const gsl_vector * y,
567 const gsl_vector * cs,
568 const gsl_matrix * M,
569 gsl_vector * c,
570 gsl_multifit_linear_workspace * work)
571 {
572 const size_t m = LQR->size1;
573 const size_t n = X->size1;
574 const size_t p = X->size2;
575
576 if (n > work->nmax || p > work->pmax)
577 {
578 GSL_ERROR("X matrix does not match workspace", GSL_EBADLEN);
579 }
580 else if (p != LQR->size2)
581 {
582 GSL_ERROR("LQR matrix does not match X", GSL_EBADLEN);
583 }
584 else if (p != c->size)
585 {
586 GSL_ERROR("c vector does not match X", GSL_EBADLEN);
587 }
588 else if (w != NULL && n != w->size)
589 {
590 GSL_ERROR("w vector does not match X", GSL_EBADLEN);
591 }
592 else if (n != y->size)
593 {
594 GSL_ERROR("y vector does not match X", GSL_EBADLEN);
595 }
596 else if (m >= p) /* square or tall L matrix */
597 {
598 if (p != cs->size)
599 {
600 GSL_ERROR("cs vector must be length p", GSL_EBADLEN);
601 }
602 else
603 {
604 int s;
605 gsl_matrix_const_view R = gsl_matrix_const_submatrix(LQR, 0, 0, p, p); /* R factor of L */
606
607 /* solve R c = cs for true solution c, using QR decomposition of L */
608 gsl_vector_memcpy(c, cs);
609 s = gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, &R.matrix, c);
610
611 return s;
612 }
613 }
614 else /* rectangular L matrix with m < p */
615 {
616 if (m != cs->size)
617 {
618 GSL_ERROR("cs vector must be length m", GSL_EBADLEN);
619 }
620 else if (n != M->size1 || p != M->size2)
621 {
622 GSL_ERROR("M matrix must be size n-by-p", GSL_EBADLEN);
623 }
624 else
625 {
626 int status;
627 const size_t pm = p - m;
628 gsl_matrix_view A = gsl_matrix_submatrix(work->A, 0, 0, n, p);
629 gsl_vector_view b = gsl_vector_subvector(work->t, 0, n);
630 gsl_matrix_view Rp = gsl_matrix_view_array(LQR->data, m, m); /* R_p */
631 gsl_matrix_view LTQR = gsl_matrix_view_array(LQR->data, p, m);
632 gsl_vector_const_view LTtau = gsl_vector_const_subvector(Ltau, 0, m);
633 gsl_matrix_const_view MQR = gsl_matrix_const_submatrix(M, 0, 0, n, pm);
634 gsl_vector_const_view Mtau = gsl_matrix_const_subcolumn(M, p - 1, 0, GSL_MIN(n, pm));
635 gsl_matrix_const_view To = gsl_matrix_const_submatrix(&MQR.matrix, 0, 0, pm, pm);
636 gsl_vector_view workp = gsl_vector_subvector(work->xt, 0, p);
637 gsl_vector_view v1, v2;
638
639 /* compute A = sqrt(W) X and b = sqrt(W) y */
640 status = gsl_multifit_linear_applyW(X, w, y, &A.matrix, &b.vector);
641 if (status)
642 return status;
643
644 /* initialize c to zero */
645 gsl_vector_set_zero(c);
646
647 /* compute c = L_inv cs = K_p R_p^{-T} cs */
648
649 /* set c(1:m) = R_p^{-T} cs */
650 v1 = gsl_vector_subvector(c, 0, m);
651 gsl_vector_memcpy(&v1.vector, cs);
652 gsl_blas_dtrsv(CblasUpper, CblasTrans, CblasNonUnit, &Rp.matrix, &v1.vector);
653
654 /* c <- K R_p^{-T} cs = [ K_p R_p^{_T} cs ; 0 ] */
655 gsl_linalg_QR_Qvec(<QR.matrix, <tau.vector, c);
656
657 /* compute: b1 = b - A L_inv cs */
658 gsl_blas_dgemv(CblasNoTrans, -1.0, &A.matrix, c, 1.0, &b.vector);
659
660 /* compute: b2 = H^T b1 */
661 gsl_linalg_QR_QTvec(&MQR.matrix, &Mtau.vector, &b.vector);
662
663 /* compute: b3 = T_o^{-1} b2 */
664 v1 = gsl_vector_subvector(&b.vector, 0, pm);
665 gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, &To.matrix, &v1.vector);
666
667 /* compute: b4 = K_o b3 */
668 gsl_vector_set_zero(&workp.vector);
669 v2 = gsl_vector_subvector(&workp.vector, m, pm);
670 gsl_vector_memcpy(&v2.vector, &v1.vector);
671 gsl_linalg_QR_Qvec(<QR.matrix, <tau.vector, &workp.vector);
672
673 /* final solution vector */
674 gsl_vector_add(c, &workp.vector);
675
676 return GSL_SUCCESS;
677 }
678 }
679 }
680
681 int
gsl_multifit_linear_genform2(const gsl_matrix * LQR,const gsl_vector * Ltau,const gsl_matrix * X,const gsl_vector * y,const gsl_vector * cs,const gsl_matrix * M,gsl_vector * c,gsl_multifit_linear_workspace * work)682 gsl_multifit_linear_genform2 (const gsl_matrix * LQR,
683 const gsl_vector * Ltau,
684 const gsl_matrix * X,
685 const gsl_vector * y,
686 const gsl_vector * cs,
687 const gsl_matrix * M,
688 gsl_vector * c,
689 gsl_multifit_linear_workspace * work)
690 {
691 int status;
692
693 status = gsl_multifit_linear_wgenform2(LQR, Ltau, X, NULL, y, cs, M, c, work);
694
695 return status;
696 }
697
698 /*
699 gsl_multifit_linear_lreg()
700 Calculate regularization parameters to use in L-curve
701 analysis
702
703 Inputs: smin - smallest singular value of LS system
704 smax - largest singular value of LS system > 0
705 reg_param - (output) vector of regularization parameters
706 derived from singular values
707
708 Return: success/error
709 */
710
711 int
gsl_multifit_linear_lreg(const double smin,const double smax,gsl_vector * reg_param)712 gsl_multifit_linear_lreg (const double smin, const double smax,
713 gsl_vector * reg_param)
714 {
715 if (smax <= 0.0)
716 {
717 GSL_ERROR("smax must be positive", GSL_EINVAL);
718 }
719 else
720 {
721 const size_t N = reg_param->size;
722
723 /* smallest regularization parameter */
724 const double smin_ratio = 16.0 * GSL_DBL_EPSILON;
725 const double new_smin = GSL_MAX(smin, smax*smin_ratio);
726 double ratio;
727 size_t i;
728
729 gsl_vector_set(reg_param, N - 1, new_smin);
730
731 /* ratio so that reg_param(1) = s(1) */
732 ratio = pow(smax / new_smin, 1.0 / ((double)N - 1.0));
733
734 /* calculate the regularization parameters */
735 for (i = N - 1; i > 0 && i--; )
736 {
737 double rp1 = gsl_vector_get(reg_param, i + 1);
738 gsl_vector_set(reg_param, i, ratio * rp1);
739 }
740
741 return GSL_SUCCESS;
742 }
743 }
744
745 /*
746 gsl_multifit_linear_lcurve()
747 Calculate L-curve using regularization parameters estimated
748 from singular values of least squares matrix
749
750 Inputs: y - right hand side vector
751 reg_param - (output) vector of regularization parameters
752 derived from singular values
753 rho - (output) vector of residual norms ||y - X c||
754 eta - (output) vector of solution norms ||lambda c||
755 work - workspace
756
757 Return: success/error
758
759 Notes:
760 1) SVD of X must be computed first by calling multifit_linear_svd();
761 work->n and work->p are initialized by this function
762 */
763
764 int
gsl_multifit_linear_lcurve(const gsl_vector * y,gsl_vector * reg_param,gsl_vector * rho,gsl_vector * eta,gsl_multifit_linear_workspace * work)765 gsl_multifit_linear_lcurve (const gsl_vector * y,
766 gsl_vector * reg_param,
767 gsl_vector * rho, gsl_vector * eta,
768 gsl_multifit_linear_workspace * work)
769 {
770 const size_t n = y->size;
771 const size_t N = rho->size; /* number of points on L-curve */
772
773 if (n != work->n)
774 {
775 GSL_ERROR("y vector does not match workspace", GSL_EBADLEN);
776 }
777 else if (N < 3)
778 {
779 GSL_ERROR ("at least 3 points are needed for L-curve analysis",
780 GSL_EBADLEN);
781 }
782 else if (N != eta->size)
783 {
784 GSL_ERROR ("size of rho and eta vectors do not match",
785 GSL_EBADLEN);
786 }
787 else if (reg_param->size != eta->size)
788 {
789 GSL_ERROR ("size of reg_param and eta vectors do not match",
790 GSL_EBADLEN);
791 }
792 else
793 {
794 int status = GSL_SUCCESS;
795 const size_t p = work->p;
796
797 size_t i, j;
798
799 gsl_matrix_view A = gsl_matrix_submatrix(work->A, 0, 0, n, p);
800 gsl_vector_view S = gsl_vector_subvector(work->S, 0, p);
801 gsl_vector_view xt = gsl_vector_subvector(work->xt, 0, p);
802 gsl_vector_view workp = gsl_matrix_subcolumn(work->QSI, 0, 0, p);
803 gsl_vector_view workp2 = gsl_vector_subvector(work->D, 0, p); /* D isn't used for regularized problems */
804
805 const double smax = gsl_vector_get(&S.vector, 0);
806 const double smin = gsl_vector_get(&S.vector, p - 1);
807
808 double dr; /* residual error from projection */
809 double normy = gsl_blas_dnrm2(y);
810 double normUTy;
811
812 /* compute projection xt = U^T y */
813 gsl_blas_dgemv (CblasTrans, 1.0, &A.matrix, y, 0.0, &xt.vector);
814
815 normUTy = gsl_blas_dnrm2(&xt.vector);
816 dr = normy*normy - normUTy*normUTy;
817
818 /* calculate regularization parameters */
819 gsl_multifit_linear_lreg(smin, smax, reg_param);
820
821 for (i = 0; i < N; ++i)
822 {
823 double lambda = gsl_vector_get(reg_param, i);
824 double lambda_sq = lambda * lambda;
825
826 for (j = 0; j < p; ++j)
827 {
828 double sj = gsl_vector_get(&S.vector, j);
829 double xtj = gsl_vector_get(&xt.vector, j);
830 double f = sj / (sj*sj + lambda_sq);
831
832 gsl_vector_set(&workp.vector, j, f * xtj);
833 gsl_vector_set(&workp2.vector, j, (1.0 - sj*f) * xtj);
834 }
835
836 gsl_vector_set(eta, i, gsl_blas_dnrm2(&workp.vector));
837 gsl_vector_set(rho, i, gsl_blas_dnrm2(&workp2.vector));
838 }
839
840 if (n > p && dr > 0.0)
841 {
842 /* add correction to residual norm (see eqs 6-7 of [1]) */
843 for (i = 0; i < N; ++i)
844 {
845 double rhoi = gsl_vector_get(rho, i);
846 double *ptr = gsl_vector_ptr(rho, i);
847
848 *ptr = sqrt(rhoi*rhoi + dr);
849 }
850 }
851
852 /* restore D to identity matrix */
853 gsl_vector_set_all(work->D, 1.0);
854
855 return status;
856 }
857 } /* gsl_multifit_linear_lcurve() */
858
859 /*
860 gsl_multifit_linear_lcurvature()
861 Calculate curvature of previously computed L-curve as a function of
862 regularization parameter
863
864 Inputs: y - right hand side vector
865 reg_param - vector of regularization parameters, length N
866 rho - vector of residual norms ||y - X c||, length N
867 eta - vector of solution norms ||c||, length N
868 kappa - (output) vector of curvature values, length N
869 work - workspace
870
871 Return: success/error
872
873 Notes:
874 1) SVD of X must be computed first by calling multifit_linear_svd();
875 work->n and work->p are initialized by this function
876
877 2) Vectors reg_param, rho, eta must be computed by calling gsl_multifit_linear_lcurve()
878 */
879
880 int
gsl_multifit_linear_lcurvature(const gsl_vector * y,const gsl_vector * reg_param,const gsl_vector * rho,const gsl_vector * eta,gsl_vector * kappa,gsl_multifit_linear_workspace * work)881 gsl_multifit_linear_lcurvature (const gsl_vector * y,
882 const gsl_vector * reg_param,
883 const gsl_vector * rho,
884 const gsl_vector * eta,
885 gsl_vector * kappa,
886 gsl_multifit_linear_workspace * work)
887 {
888 const size_t n = y->size;
889 const size_t N = rho->size; /* number of points on L-curve */
890
891 if (n != work->n)
892 {
893 GSL_ERROR("y vector does not match workspace", GSL_EBADLEN);
894 }
895 else if (N != eta->size)
896 {
897 GSL_ERROR ("size of rho and eta vectors do not match",
898 GSL_EBADLEN);
899 }
900 else if (reg_param->size != N)
901 {
902 GSL_ERROR ("size of reg_param and rho vectors do not match",
903 GSL_EBADLEN);
904 }
905 else if (kappa->size != N)
906 {
907 GSL_ERROR ("size of reg_param and rho vectors do not match",
908 GSL_EBADLEN);
909 }
910 else
911 {
912 int status = GSL_SUCCESS;
913 const size_t p = work->p;
914 gsl_matrix_view U = gsl_matrix_submatrix(work->A, 0, 0, n, p);
915 gsl_vector_view S = gsl_vector_subvector(work->S, 0, p);
916 gsl_vector_view beta = gsl_vector_subvector(work->xt, 0, p);
917 size_t i;
918
919 /* compute projection beta = U^T y */
920 gsl_blas_dgemv (CblasTrans, 1.0, &U.matrix, y, 0.0, &beta.vector);
921
922 for (i = 0; i < N; ++i)
923 {
924 double lambda = gsl_vector_get(reg_param, i);
925 double lambda_sq = lambda * lambda;
926 double eta_i = gsl_vector_get(eta, i);
927 double rho_i = gsl_vector_get(rho, i);
928 double phi_i = 0.0, dphi_i = 0.0;
929 double psi_i = 0.0, dpsi_i = 0.0;
930 double deta_i, ddeta_i, drho_i, ddrho_i;
931 double dlogeta_i, ddlogeta_i, dlogrho_i, ddlogrho_i;
932 double kappa_i;
933 size_t j;
934
935 for (j = 0; j < p; ++j)
936 {
937 double beta_j = gsl_vector_get(&beta.vector, j);
938 double s_j = gsl_vector_get(&S.vector, j);
939 double sj_sq = s_j * s_j;
940 double f_j = sj_sq / (sj_sq + lambda_sq); /* filter factor */
941 double onemf_j = 1.0 - f_j;
942 double f1_j = -2.0 * f_j * onemf_j / lambda;
943 double f2_j = -f1_j * (3.0 - 4.0 * f_j) / lambda;
944 double xi_j = beta_j / s_j;
945
946 phi_i += f_j * f1_j * xi_j * xi_j;
947 psi_i += onemf_j * f1_j * beta_j * beta_j;
948
949 dphi_i += (f1_j * f1_j + f_j * f2_j) * xi_j * xi_j;
950 dpsi_i += (-f1_j * f1_j + onemf_j * f2_j) * beta_j * beta_j;
951 }
952
953 deta_i = phi_i / eta_i;
954 drho_i = -psi_i / rho_i;
955 ddeta_i = dphi_i / eta_i - deta_i * (deta_i / eta_i);
956 ddrho_i = -dpsi_i / rho_i - drho_i * (drho_i / rho_i);
957
958 dlogeta_i = deta_i / eta_i;
959 dlogrho_i = drho_i / rho_i;
960 ddlogeta_i = ddeta_i / eta_i - dlogeta_i * dlogeta_i;
961 ddlogrho_i = ddrho_i / rho_i - dlogrho_i * dlogrho_i;
962
963 kappa_i = (dlogrho_i * ddlogeta_i - ddlogrho_i * dlogeta_i) /
964 pow(dlogrho_i * dlogrho_i + dlogeta_i * dlogeta_i, 1.5);
965 gsl_vector_set(kappa, i, kappa_i);
966 }
967
968 return status;
969 }
970 }
971
972 /*
973 gsl_multifit_linear_lcorner()
974 Determine point on L-curve of maximum curvature. For each
975 set of 3 points on the L-curve, the circle which passes through
976 the 3 points is computed. The radius of the circle is then used
977 as an estimate of the curvature at the middle point. The point
978 with maximum curvature is then selected.
979
980 Inputs: rho - vector of residual norms ||A x - b||
981 eta - vector of solution norms ||L x||
982 idx - (output) index i such that
983 (log(rho(i)),log(eta(i))) is the point of
984 maximum curvature
985
986 Return: success/error
987 */
988
989 int
gsl_multifit_linear_lcorner(const gsl_vector * rho,const gsl_vector * eta,size_t * idx)990 gsl_multifit_linear_lcorner(const gsl_vector *rho,
991 const gsl_vector *eta,
992 size_t *idx)
993 {
994 const size_t n = rho->size;
995
996 if (n < 3)
997 {
998 GSL_ERROR ("at least 3 points are needed for L-curve analysis",
999 GSL_EBADLEN);
1000 }
1001 else if (n != eta->size)
1002 {
1003 GSL_ERROR ("size of rho and eta vectors do not match",
1004 GSL_EBADLEN);
1005 }
1006 else
1007 {
1008 int s = GSL_SUCCESS;
1009 size_t i;
1010 double x1, y1; /* first point of triangle on L-curve */
1011 double x2, y2; /* second point of triangle on L-curve */
1012 double rmin = -1.0; /* minimum radius of curvature */
1013
1014 /* initial values */
1015 x1 = log(gsl_vector_get(rho, 0));
1016 y1 = log(gsl_vector_get(eta, 0));
1017
1018 x2 = log(gsl_vector_get(rho, 1));
1019 y2 = log(gsl_vector_get(eta, 1));
1020
1021 for (i = 1; i < n - 1; ++i)
1022 {
1023 /*
1024 * The points (x1,y1), (x2,y2), (x3,y3) are the previous,
1025 * current, and next point on the L-curve. We will find
1026 * the circle which fits these 3 points and take its radius
1027 * as an estimate of the curvature at this point.
1028 */
1029 double x3 = log(gsl_vector_get(rho, i + 1));
1030 double y3 = log(gsl_vector_get(eta, i + 1));
1031
1032 double x21 = x2 - x1;
1033 double y21 = y2 - y1;
1034 double x31 = x3 - x1;
1035 double y31 = y3 - y1;
1036 double h21 = x21*x21 + y21*y21;
1037 double h31 = x31*x31 + y31*y31;
1038 double d = fabs(2.0 * (x21*y31 - x31*y21));
1039 double r = sqrt(h21*h31*((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2))) / d;
1040
1041 /* if d =~ 0 then there are nearly colinear points */
1042 if (gsl_finite(r))
1043 {
1044 /* check for smallest radius of curvature */
1045 if (r < rmin || rmin < 0.0)
1046 {
1047 rmin = r;
1048 *idx = i;
1049 }
1050 }
1051
1052 /* update previous/current L-curve values */
1053 x1 = x2;
1054 y1 = y2;
1055 x2 = x3;
1056 y2 = y3;
1057 }
1058
1059 /* check if a minimum radius was found */
1060 if (rmin < 0.0)
1061 {
1062 /* possibly co-linear points */
1063 GSL_ERROR("failed to find minimum radius", GSL_EINVAL);
1064 }
1065
1066 return s;
1067 }
1068 } /* gsl_multifit_linear_lcorner() */
1069
1070 /*
1071 gsl_multifit_linear_lcorner2()
1072 Determine point on L-curve (lambda^2, ||c||^2) of maximum curvature.
1073 For each set of 3 points on the L-curve, the circle which passes through
1074 the 3 points is computed. The radius of the circle is then used
1075 as an estimate of the curvature at the middle point. The point
1076 with maximum curvature is then selected.
1077
1078 This routine is based on the paper
1079
1080 M. Rezghi and S. M. Hosseini, "A new variant of L-curve for Tikhonov
1081 regularization", J. Comp. App. Math., 231 (2009).
1082
1083 Inputs: reg_param - vector of regularization parameters
1084 eta - vector of solution norms ||L x||
1085 idx - (output) index i such that
1086 (lambda(i)^2,eta(i)^2) is the point of
1087 maximum curvature
1088
1089 Return: success/error
1090 */
1091
1092 int
gsl_multifit_linear_lcorner2(const gsl_vector * reg_param,const gsl_vector * eta,size_t * idx)1093 gsl_multifit_linear_lcorner2(const gsl_vector *reg_param,
1094 const gsl_vector *eta,
1095 size_t *idx)
1096 {
1097 const size_t n = reg_param->size;
1098
1099 if (n < 3)
1100 {
1101 GSL_ERROR ("at least 3 points are needed for L-curve analysis",
1102 GSL_EBADLEN);
1103 }
1104 else if (n != eta->size)
1105 {
1106 GSL_ERROR ("size of reg_param and eta vectors do not match",
1107 GSL_EBADLEN);
1108 }
1109 else
1110 {
1111 int s = GSL_SUCCESS;
1112 size_t i;
1113 double x1, y1; /* first point of triangle on L-curve */
1114 double x2, y2; /* second point of triangle on L-curve */
1115 double rmin = -1.0; /* minimum radius of curvature */
1116
1117 /* initial values */
1118 x1 = gsl_vector_get(reg_param, 0);
1119 x1 *= x1;
1120 y1 = gsl_vector_get(eta, 0);
1121 y1 *= y1;
1122
1123 x2 = gsl_vector_get(reg_param, 1);
1124 x2 *= x2;
1125 y2 = gsl_vector_get(eta, 1);
1126 y2 *= y2;
1127
1128 for (i = 1; i < n - 1; ++i)
1129 {
1130 /*
1131 * The points (x1,y1), (x2,y2), (x3,y3) are the previous,
1132 * current, and next point on the L-curve. We will find
1133 * the circle which fits these 3 points and take its radius
1134 * as an estimate of the curvature at this point.
1135 */
1136 double lamip1 = gsl_vector_get(reg_param, i + 1);
1137 double etaip1 = gsl_vector_get(eta, i + 1);
1138 double x3 = lamip1 * lamip1;
1139 double y3 = etaip1 * etaip1;
1140
1141 double x21 = x2 - x1;
1142 double y21 = y2 - y1;
1143 double x31 = x3 - x1;
1144 double y31 = y3 - y1;
1145 double h21 = x21*x21 + y21*y21;
1146 double h31 = x31*x31 + y31*y31;
1147 double d = fabs(2.0 * (x21*y31 - x31*y21));
1148 double r = sqrt(h21*h31*((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2))) / d;
1149
1150 /* if d =~ 0 then there are nearly colinear points */
1151 if (gsl_finite(r))
1152 {
1153 /* check for smallest radius of curvature */
1154 if (r < rmin || rmin < 0.0)
1155 {
1156 rmin = r;
1157 *idx = i;
1158 }
1159 }
1160
1161 /* update previous/current L-curve values */
1162 x1 = x2;
1163 y1 = y2;
1164 x2 = x3;
1165 y2 = y3;
1166 }
1167
1168 /* check if a minimum radius was found */
1169 if (rmin < 0.0)
1170 {
1171 /* possibly co-linear points */
1172 GSL_ERROR("failed to find minimum radius", GSL_EINVAL);
1173 }
1174
1175 return s;
1176 }
1177 } /* gsl_multifit_linear_lcorner2() */
1178
1179 #define GSL_MULTIFIT_MAXK 100
1180
1181 /*
1182 gsl_multifit_linear_L()
1183 Compute discrete approximation to derivative operator of order
1184 k on a regular grid of p points, ie: L is (p-k)-by-p
1185 */
1186
1187 int
gsl_multifit_linear_Lk(const size_t p,const size_t k,gsl_matrix * L)1188 gsl_multifit_linear_Lk(const size_t p, const size_t k, gsl_matrix *L)
1189 {
1190 if (p <= k)
1191 {
1192 GSL_ERROR("p must be larger than derivative order", GSL_EBADLEN);
1193 }
1194 else if (k >= GSL_MULTIFIT_MAXK - 1)
1195 {
1196 GSL_ERROR("derivative order k too large", GSL_EBADLEN);
1197 }
1198 else if (p - k != L->size1 || p != L->size2)
1199 {
1200 GSL_ERROR("L matrix must be (p-k)-by-p", GSL_EBADLEN);
1201 }
1202 else
1203 {
1204 double c_data[GSL_MULTIFIT_MAXK];
1205 gsl_vector_view cv = gsl_vector_view_array(c_data, k + 1);
1206 size_t i, j;
1207
1208 /* zeroth derivative */
1209 if (k == 0)
1210 {
1211 gsl_matrix_set_identity(L);
1212 return GSL_SUCCESS;
1213 }
1214
1215 gsl_matrix_set_zero(L);
1216
1217 gsl_vector_set_zero(&cv.vector);
1218 gsl_vector_set(&cv.vector, 0, -1.0);
1219 gsl_vector_set(&cv.vector, 1, 1.0);
1220
1221 for (i = 1; i < k; ++i)
1222 {
1223 double cjm1 = 0.0;
1224
1225 for (j = 0; j < k + 1; ++j)
1226 {
1227 double cj = gsl_vector_get(&cv.vector, j);
1228
1229 gsl_vector_set(&cv.vector, j, cjm1 - cj);
1230 cjm1 = cj;
1231 }
1232 }
1233
1234 /* build L, the c_i are the entries on the diagonals */
1235 for (i = 0; i < k + 1; ++i)
1236 {
1237 gsl_vector_view v = gsl_matrix_superdiagonal(L, i);
1238 double ci = gsl_vector_get(&cv.vector, i);
1239
1240 gsl_vector_set_all(&v.vector, ci);
1241 }
1242
1243 return GSL_SUCCESS;
1244 }
1245 } /* gsl_multifit_linear_Lk() */
1246
1247 /*
1248 gsl_multifit_linear_Lsobolev()
1249 Construct Sobolev smoothing norm operator
1250
1251 L = [ a_0 I; a_1 L_1; a_2 L_2; ...; a_k L_k ]
1252
1253 by computing the Cholesky factor of L^T L
1254
1255 Inputs: p - number of columns of L
1256 kmax - maximum derivative order (< p)
1257 alpha - vector of weights; alpha_k multiplies L_k, size kmax + 1
1258 L - (output) upper triangular Sobolev matrix p-by-p,
1259 stored in upper triangle
1260 work - workspace
1261
1262 Notes:
1263 1) work->Q is used to store intermediate L_k matrices
1264 */
1265
1266 int
gsl_multifit_linear_Lsobolev(const size_t p,const size_t kmax,const gsl_vector * alpha,gsl_matrix * L,gsl_multifit_linear_workspace * work)1267 gsl_multifit_linear_Lsobolev(const size_t p, const size_t kmax,
1268 const gsl_vector *alpha, gsl_matrix *L,
1269 gsl_multifit_linear_workspace *work)
1270 {
1271 if (p > work->pmax)
1272 {
1273 GSL_ERROR("p is larger than workspace", GSL_EBADLEN);
1274 }
1275 else if (p <= kmax)
1276 {
1277 GSL_ERROR("p must be larger than derivative order", GSL_EBADLEN);
1278 }
1279 else if (kmax + 1 != alpha->size)
1280 {
1281 GSL_ERROR("alpha must be size kmax + 1", GSL_EBADLEN);
1282 }
1283 else if (p != L->size1)
1284 {
1285 GSL_ERROR("L matrix is wrong size", GSL_EBADLEN);
1286 }
1287 else if (L->size1 != L->size2)
1288 {
1289 GSL_ERROR("L matrix is not square", GSL_ENOTSQR);
1290 }
1291 else
1292 {
1293 int s;
1294 size_t j, k;
1295 gsl_vector_view d = gsl_matrix_diagonal(L);
1296 const double alpha0 = gsl_vector_get(alpha, 0);
1297
1298 /* initialize L to alpha0^2 I */
1299 gsl_matrix_set_zero(L);
1300 gsl_vector_add_constant(&d.vector, alpha0 * alpha0);
1301
1302 for (k = 1; k <= kmax; ++k)
1303 {
1304 gsl_matrix_view Lk = gsl_matrix_submatrix(work->Q, 0, 0, p - k, p);
1305 double ak = gsl_vector_get(alpha, k);
1306
1307 /* compute a_k L_k */
1308 s = gsl_multifit_linear_Lk(p, k, &Lk.matrix);
1309 if (s)
1310 return s;
1311
1312 gsl_matrix_scale(&Lk.matrix, ak);
1313
1314 /* LTL += L_k^T L_k */
1315 gsl_blas_dsyrk(CblasLower, CblasTrans, 1.0, &Lk.matrix, 1.0, L);
1316 }
1317
1318 s = gsl_linalg_cholesky_decomp(L);
1319 if (s)
1320 return s;
1321
1322 /* copy Cholesky factor to upper triangle and zero out bottom */
1323 gsl_matrix_transpose_tricpy(CblasLower, CblasUnit, L, L);
1324
1325 for (j = 0; j < p; ++j)
1326 {
1327 for (k = 0; k < j; ++k)
1328 gsl_matrix_set(L, j, k, 0.0);
1329 }
1330
1331 return GSL_SUCCESS;
1332 }
1333 }
1334