1 /* Cholesky Decomposition
2 *
3 * Copyright (C) 2000 Thomas Walter
4 * Copyright (C) 2000, 2001, 2002, 2003, 2005, 2007 Brian Gough, Gerard Jungman
5 * Copyright (C) 2016, 2019 Patrick Alken
6 *
7 * This is free software; you can redistribute it and/or modify it
8 * under the terms of the GNU General Public License as published by the
9 * Free Software Foundation; either version 3, or (at your option) any
10 * later version.
11 *
12 * This source is distributed in the hope that it will be useful, but WITHOUT
13 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 * for more details.
16 *
17 * 03 May 2000: Modified for GSL by Brian Gough
18 * 29 Jul 2005: Additions by Gerard Jungman
19 * 04 Mar 2016: Change Cholesky algorithm to gaxpy version by Patrick Alken
20 * 26 May 2019: implement recursive Cholesky with Level 3 BLAS by Patrick Alken
21 */
22
23 /*
24 * Cholesky decomposition of a symmetric positive definite matrix.
25 *
26 * This algorithm does:
27 * A = L * L'
28 * with
29 * L := lower left triangle matrix
30 * L' := the transposed form of L.
31 *
32 */
33
34 #include <config.h>
35
36 #include <gsl/gsl_math.h>
37 #include <gsl/gsl_errno.h>
38 #include <gsl/gsl_vector.h>
39 #include <gsl/gsl_matrix.h>
40 #include <gsl/gsl_blas.h>
41 #include <gsl/gsl_linalg.h>
42
43 #include "recurse.h"
44
45 static double cholesky_norm1(const gsl_matrix * LLT, gsl_vector * work);
46 static int cholesky_Ainv(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params);
47 static int cholesky_decomp_L2 (gsl_matrix * A);
48 static int cholesky_decomp_L3 (gsl_matrix * A);
49
50 /*
51 In GSL 2.2, we decided to modify the behavior of the Cholesky decomposition
52 to store the Cholesky factor in the lower triangle, and store the original
53 matrix in the upper triangle. Previous versions stored the Cholesky factor in
54 both places. The routine gsl_linalg_cholesky_decomp1 was added for the new
55 behavior, and gsl_linalg_cholesky_decomp is maintained for backward compatibility.
56 It will be removed in a future release.
57 */
58
59 int
gsl_linalg_cholesky_decomp(gsl_matrix * A)60 gsl_linalg_cholesky_decomp (gsl_matrix * A)
61 {
62 int status;
63
64 status = gsl_linalg_cholesky_decomp1(A);
65 if (status == GSL_SUCCESS)
66 {
67 gsl_matrix_transpose_tricpy(CblasLower, CblasUnit, A, A);
68 }
69
70 return status;
71 }
72
73 /*
74 gsl_linalg_cholesky_decomp1()
75 Perform Cholesky decomposition of a symmetric positive
76 definite matrix using lower triangle using Level 3 BLAS algorithm.
77
78 Inputs: A - (input) symmetric, positive definite matrix
79 (output) lower triangle contains Cholesky factor
80
81 Return: success/error
82
83 Notes:
84 1) original matrix is saved in upper triangle on output
85 */
86
87 int
gsl_linalg_cholesky_decomp1(gsl_matrix * A)88 gsl_linalg_cholesky_decomp1 (gsl_matrix * A)
89 {
90 const size_t N = A->size1;
91
92 if (N != A->size2)
93 {
94 GSL_ERROR("Cholesky decomposition requires square matrix", GSL_ENOTSQR);
95 }
96 else
97 {
98 /* save original matrix in upper triangle for later rcond calculation */
99 gsl_matrix_transpose_tricpy(CblasLower, CblasUnit, A, A);
100
101 return cholesky_decomp_L3(A);
102 }
103 }
104
105 int
gsl_linalg_cholesky_solve(const gsl_matrix * LLT,const gsl_vector * b,gsl_vector * x)106 gsl_linalg_cholesky_solve (const gsl_matrix * LLT,
107 const gsl_vector * b,
108 gsl_vector * x)
109 {
110 if (LLT->size1 != LLT->size2)
111 {
112 GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
113 }
114 else if (LLT->size1 != b->size)
115 {
116 GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
117 }
118 else if (LLT->size2 != x->size)
119 {
120 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
121 }
122 else
123 {
124 int status;
125
126 /* copy x <- b */
127 gsl_vector_memcpy (x, b);
128
129 status = gsl_linalg_cholesky_svx(LLT, x);
130
131 return status;
132 }
133 }
134
135 int
gsl_linalg_cholesky_svx(const gsl_matrix * LLT,gsl_vector * x)136 gsl_linalg_cholesky_svx (const gsl_matrix * LLT,
137 gsl_vector * x)
138 {
139 if (LLT->size1 != LLT->size2)
140 {
141 GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
142 }
143 else if (LLT->size2 != x->size)
144 {
145 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
146 }
147 else
148 {
149 /* solve for c using forward-substitution, L c = b */
150 gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasNonUnit, LLT, x);
151
152 /* perform back-substitution, L^T x = c */
153 gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LLT, x);
154
155 return GSL_SUCCESS;
156 }
157 }
158
159 int
gsl_linalg_cholesky_solve_mat(const gsl_matrix * LLT,const gsl_matrix * B,gsl_matrix * X)160 gsl_linalg_cholesky_solve_mat (const gsl_matrix * LLT,
161 const gsl_matrix * B,
162 gsl_matrix * X)
163 {
164 if (LLT->size1 != LLT->size2)
165 {
166 GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
167 }
168 else if (LLT->size1 != B->size1)
169 {
170 GSL_ERROR ("matrix size must match B size", GSL_EBADLEN);
171 }
172 else if (LLT->size2 != X->size1)
173 {
174 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
175 }
176 else
177 {
178 int status;
179
180 /* copy X <- B */
181 gsl_matrix_memcpy (X, B);
182
183 status = gsl_linalg_cholesky_svx_mat(LLT, X);
184
185 return status;
186 }
187 }
188
189 int
gsl_linalg_cholesky_svx_mat(const gsl_matrix * LLT,gsl_matrix * X)190 gsl_linalg_cholesky_svx_mat (const gsl_matrix * LLT,
191 gsl_matrix * X)
192 {
193 if (LLT->size1 != LLT->size2)
194 {
195 GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
196 }
197 else if (LLT->size2 != X->size1)
198 {
199 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
200 }
201 else
202 {
203 /* solve for C using forward-substitution, L C = B */
204 gsl_blas_dtrsm (CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit, 1.0,
205 LLT, X);
206
207 /* perform back-substitution, L^T X = C */
208 gsl_blas_dtrsm (CblasLeft, CblasLower, CblasTrans, CblasNonUnit, 1.0,
209 LLT, X);
210
211 return GSL_SUCCESS;
212 }
213 }
214
215 /*
216 gsl_linalg_cholesky_invert()
217 Compute the inverse of a symmetric positive definite matrix in
218 Cholesky form.
219
220 Inputs: LLT - matrix in cholesky form on input
221 A^{-1} = L^{-t} L^{-1} on output
222
223 Return: success or error
224 */
225
226 int
gsl_linalg_cholesky_invert(gsl_matrix * LLT)227 gsl_linalg_cholesky_invert(gsl_matrix * LLT)
228 {
229 if (LLT->size1 != LLT->size2)
230 {
231 GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
232 }
233 else
234 {
235 int status;
236
237 /* invert the lower triangle of LLT */
238 status = gsl_linalg_tri_invert(CblasLower, CblasNonUnit, LLT);
239 if (status)
240 return status;
241
242 /* compute A^{-1} = L^{-T} L^{-1} */
243 status = gsl_linalg_tri_LTL(LLT);
244 if (status)
245 return status;
246
247 /* copy lower triangle to upper */
248 gsl_matrix_transpose_tricpy(CblasLower, CblasUnit, LLT, LLT);
249
250 return GSL_SUCCESS;
251 }
252 }
253
254 int
gsl_linalg_cholesky_decomp_unit(gsl_matrix * A,gsl_vector * D)255 gsl_linalg_cholesky_decomp_unit(gsl_matrix * A, gsl_vector * D)
256 {
257 const size_t N = A->size1;
258 size_t i, j;
259
260 /* initial Cholesky */
261 int stat_chol = gsl_linalg_cholesky_decomp1(A);
262
263 if(stat_chol == GSL_SUCCESS)
264 {
265 /* calculate D from diagonal part of initial Cholesky */
266 for(i = 0; i < N; ++i)
267 {
268 const double C_ii = gsl_matrix_get(A, i, i);
269 gsl_vector_set(D, i, C_ii*C_ii);
270 }
271
272 /* multiply initial Cholesky by 1/sqrt(D) on the right */
273 for(i = 0; i < N; ++i)
274 {
275 for(j = 0; j < N; ++j)
276 {
277 gsl_matrix_set(A, i, j, gsl_matrix_get(A, i, j) / sqrt(gsl_vector_get(D, j)));
278 }
279 }
280
281 /* Because the initial Cholesky contained both L and transpose(L),
282 the result of the multiplication is not symmetric anymore;
283 but the lower triangle _is_ correct. Therefore we reflect
284 it to the upper triangle and declare victory.
285 */
286 for(i = 0; i < N; ++i)
287 for(j = i + 1; j < N; ++j)
288 gsl_matrix_set(A, i, j, gsl_matrix_get(A, j, i));
289 }
290
291 return stat_chol;
292 }
293
294 /*
295 gsl_linalg_cholesky_scale()
296 This function computes scale factors diag(S), such that
297
298 diag(S) A diag(S)
299
300 has a condition number within a factor N of the matrix
301 with the smallest condition number over all possible
302 diagonal scalings. See Corollary 7.6 of:
303
304 N. J. Higham, Accuracy and Stability of Numerical Algorithms (2nd Edition),
305 SIAM, 2002.
306
307 Inputs: A - symmetric positive definite matrix
308 S - (output) scale factors, S_i = 1 / sqrt(A_ii)
309 */
310
311 int
gsl_linalg_cholesky_scale(const gsl_matrix * A,gsl_vector * S)312 gsl_linalg_cholesky_scale(const gsl_matrix * A, gsl_vector * S)
313 {
314 const size_t M = A->size1;
315 const size_t N = A->size2;
316
317 if (M != N)
318 {
319 GSL_ERROR("A is not a square matrix", GSL_ENOTSQR);
320 }
321 else if (N != S->size)
322 {
323 GSL_ERROR("S must have length N", GSL_EBADLEN);
324 }
325 else
326 {
327 size_t i;
328
329 /* compute S_i = 1/sqrt(A_{ii}) */
330 for (i = 0; i < N; ++i)
331 {
332 double Aii = gsl_matrix_get(A, i, i);
333
334 if (Aii <= 0.0)
335 gsl_vector_set(S, i, 1.0); /* matrix not positive definite */
336 else
337 gsl_vector_set(S, i, 1.0 / sqrt(Aii));
338 }
339
340 return GSL_SUCCESS;
341 }
342 }
343
344 /*
345 gsl_linalg_cholesky_scale_apply()
346 This function applies scale transformation to A:
347
348 A <- diag(S) A diag(S)
349
350 Inputs: A - (input/output)
351 on input, symmetric positive definite matrix
352 on output, diag(S) * A * diag(S) in lower triangle
353 S - (input) scale factors
354 */
355
356 int
gsl_linalg_cholesky_scale_apply(gsl_matrix * A,const gsl_vector * S)357 gsl_linalg_cholesky_scale_apply(gsl_matrix * A, const gsl_vector * S)
358 {
359 const size_t M = A->size1;
360 const size_t N = A->size2;
361
362 if (M != N)
363 {
364 GSL_ERROR("A is not a square matrix", GSL_ENOTSQR);
365 }
366 else if (N != S->size)
367 {
368 GSL_ERROR("S must have length N", GSL_EBADLEN);
369 }
370 else
371 {
372 size_t i, j;
373
374 /* compute: A <- diag(S) A diag(S) using lower triangle */
375 for (j = 0; j < N; ++j)
376 {
377 double sj = gsl_vector_get(S, j);
378
379 for (i = j; i < N; ++i)
380 {
381 double si = gsl_vector_get(S, i);
382 double *Aij = gsl_matrix_ptr(A, i, j);
383 *Aij *= si * sj;
384 }
385 }
386
387 return GSL_SUCCESS;
388 }
389 }
390
391 int
gsl_linalg_cholesky_decomp2(gsl_matrix * A,gsl_vector * S)392 gsl_linalg_cholesky_decomp2(gsl_matrix * A, gsl_vector * S)
393 {
394 const size_t M = A->size1;
395 const size_t N = A->size2;
396
397 if (M != N)
398 {
399 GSL_ERROR("cholesky decomposition requires square matrix", GSL_ENOTSQR);
400 }
401 else if (N != S->size)
402 {
403 GSL_ERROR("S must have length N", GSL_EBADLEN);
404 }
405 else
406 {
407 int status;
408
409 /* compute scaling factors to reduce cond(A) */
410 status = gsl_linalg_cholesky_scale(A, S);
411 if (status)
412 return status;
413
414 /* apply scaling factors */
415 status = gsl_linalg_cholesky_scale_apply(A, S);
416 if (status)
417 return status;
418
419 /* compute Cholesky decomposition of diag(S) A diag(S) */
420 status = gsl_linalg_cholesky_decomp1(A);
421 if (status)
422 return status;
423
424 return GSL_SUCCESS;
425 }
426 }
427
428 int
gsl_linalg_cholesky_svx2(const gsl_matrix * LLT,const gsl_vector * S,gsl_vector * x)429 gsl_linalg_cholesky_svx2 (const gsl_matrix * LLT,
430 const gsl_vector * S,
431 gsl_vector * x)
432 {
433 if (LLT->size1 != LLT->size2)
434 {
435 GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
436 }
437 else if (LLT->size2 != S->size)
438 {
439 GSL_ERROR ("matrix size must match S", GSL_EBADLEN);
440 }
441 else if (LLT->size2 != x->size)
442 {
443 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
444 }
445 else
446 {
447 /* b~ = diag(S) b */
448 gsl_vector_mul(x, S);
449
450 /* Solve for c using forward-substitution, L c = b~ */
451 gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasNonUnit, LLT, x);
452
453 /* Perform back-substitution, L^T x~ = c */
454 gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LLT, x);
455
456 /* compute original solution vector x = S x~ */
457 gsl_vector_mul(x, S);
458
459 return GSL_SUCCESS;
460 }
461 }
462
463 int
gsl_linalg_cholesky_solve2(const gsl_matrix * LLT,const gsl_vector * S,const gsl_vector * b,gsl_vector * x)464 gsl_linalg_cholesky_solve2 (const gsl_matrix * LLT,
465 const gsl_vector * S,
466 const gsl_vector * b,
467 gsl_vector * x)
468 {
469 if (LLT->size1 != LLT->size2)
470 {
471 GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
472 }
473 else if (LLT->size1 != S->size)
474 {
475 GSL_ERROR ("matrix size must match S size", GSL_EBADLEN);
476 }
477 else if (LLT->size1 != b->size)
478 {
479 GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
480 }
481 else if (LLT->size2 != x->size)
482 {
483 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
484 }
485 else
486 {
487 int status;
488
489 /* Copy x <- b */
490 gsl_vector_memcpy (x, b);
491
492 status = gsl_linalg_cholesky_svx2(LLT, S, x);
493
494 return status;
495 }
496 }
497
498 int
gsl_linalg_cholesky_rcond(const gsl_matrix * LLT,double * rcond,gsl_vector * work)499 gsl_linalg_cholesky_rcond (const gsl_matrix * LLT, double * rcond,
500 gsl_vector * work)
501 {
502 const size_t M = LLT->size1;
503 const size_t N = LLT->size2;
504
505 if (M != N)
506 {
507 GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
508 }
509 else if (work->size != 3 * N)
510 {
511 GSL_ERROR ("work vector must have length 3*N", GSL_EBADLEN);
512 }
513 else
514 {
515 int status;
516 double Anorm = cholesky_norm1(LLT, work); /* ||A||_1 */
517 double Ainvnorm; /* ||A^{-1}||_1 */
518
519 *rcond = 0.0;
520
521 /* don't continue if matrix is singular */
522 if (Anorm == 0.0)
523 return GSL_SUCCESS;
524
525 /* estimate ||A^{-1}||_1 */
526 status = gsl_linalg_invnorm1(N, cholesky_Ainv, (void *) LLT, &Ainvnorm, work);
527
528 if (status)
529 return status;
530
531 if (Ainvnorm != 0.0)
532 *rcond = (1.0 / Anorm) / Ainvnorm;
533
534 return GSL_SUCCESS;
535 }
536 }
537
538 /* compute 1-norm of original matrix, stored in upper triangle of LLT;
539 * diagonal entries have to be reconstructed */
540 static double
cholesky_norm1(const gsl_matrix * LLT,gsl_vector * work)541 cholesky_norm1(const gsl_matrix * LLT, gsl_vector * work)
542 {
543 const size_t N = LLT->size1;
544 double max = 0.0;
545 size_t i, j;
546
547 for (j = 0; j < N; ++j)
548 {
549 double sum = 0.0;
550 gsl_vector_const_view lj = gsl_matrix_const_subrow(LLT, j, 0, j + 1);
551 double Ajj;
552
553 /* compute diagonal (j,j) entry of A */
554 gsl_blas_ddot(&lj.vector, &lj.vector, &Ajj);
555
556 for (i = 0; i < j; ++i)
557 {
558 double *wi = gsl_vector_ptr(work, i);
559 double Aij = gsl_matrix_get(LLT, i, j);
560 double absAij = fabs(Aij);
561
562 sum += absAij;
563 *wi += absAij;
564 }
565
566 gsl_vector_set(work, j, sum + fabs(Ajj));
567 }
568
569 for (i = 0; i < N; ++i)
570 {
571 double wi = gsl_vector_get(work, i);
572 max = GSL_MAX(max, wi);
573 }
574
575 return max;
576 }
577
578 /* x := A^{-1} x = A^{-t} x, A = L L^T */
579 static int
cholesky_Ainv(CBLAS_TRANSPOSE_t TransA,gsl_vector * x,void * params)580 cholesky_Ainv(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params)
581 {
582 int status;
583 gsl_matrix * A = (gsl_matrix * ) params;
584
585 (void) TransA; /* unused parameter warning */
586
587 /* compute L^{-1} x */
588 status = gsl_blas_dtrsv(CblasLower, CblasNoTrans, CblasNonUnit, A, x);
589 if (status)
590 return status;
591
592 /* compute L^{-t} x */
593 status = gsl_blas_dtrsv(CblasLower, CblasTrans, CblasNonUnit, A, x);
594 if (status)
595 return status;
596
597 return GSL_SUCCESS;
598 }
599
600 /*
601 cholesky_decomp_L2()
602 Perform Cholesky decomposition of a symmetric positive
603 definite matrix using lower triangle
604
605 Inputs: A - (input) symmetric, positive definite matrix
606 (output) lower triangle contains Cholesky factor
607
608 Return: success/error
609
610 Notes:
611 1) Based on algorithm 4.2.1 (Gaxpy Cholesky) of Golub and
612 Van Loan, Matrix Computations (4th ed), using Level 2 BLAS.
613 */
614
615 static int
cholesky_decomp_L2(gsl_matrix * A)616 cholesky_decomp_L2 (gsl_matrix * A)
617 {
618 const size_t N = A->size1;
619
620 if (N != A->size2)
621 {
622 GSL_ERROR("Cholesky decomposition requires square matrix", GSL_ENOTSQR);
623 }
624 else
625 {
626 size_t j;
627
628 for (j = 0; j < N; ++j)
629 {
630 double ajj;
631 gsl_vector_view v = gsl_matrix_subcolumn(A, j, j, N - j); /* A(j:n,j) */
632
633 if (j > 0)
634 {
635 gsl_vector_view w = gsl_matrix_subrow(A, j, 0, j); /* A(j,1:j-1)^T */
636 gsl_matrix_view m = gsl_matrix_submatrix(A, j, 0, N - j, j); /* A(j:n,1:j-1) */
637
638 gsl_blas_dgemv(CblasNoTrans, -1.0, &m.matrix, &w.vector, 1.0, &v.vector);
639 }
640
641 ajj = gsl_matrix_get(A, j, j);
642
643 if (ajj <= 0.0)
644 {
645 GSL_ERROR("matrix is not positive definite", GSL_EDOM);
646 }
647
648 ajj = sqrt(ajj);
649 gsl_vector_scale(&v.vector, 1.0 / ajj);
650 }
651
652 return GSL_SUCCESS;
653 }
654 }
655
656 /*
657 cholesky_decomp_L3()
658 Perform Cholesky decomposition of a symmetric positive
659 definite matrix using Level 3 BLAS.
660
661 Inputs: A - (input) symmetric, positive definite matrix in lower triangle
662 (output) lower triangle contains Cholesky factor
663
664 Return: success/error
665
666 Notes:
667 1) Based on ReLAPACK recursive block Cholesky algorithm using Level 3 BLAS
668
669 2) 28 May 2019: performed several benchmark tests of this recursive variant
670 against the right-looking block variant from LAPACK. This recursive variant
671 performed faster in all cases, so it is now the default algorithm.
672 */
673
674 static int
cholesky_decomp_L3(gsl_matrix * A)675 cholesky_decomp_L3 (gsl_matrix * A)
676 {
677 const size_t N = A->size1;
678
679 if (N != A->size2)
680 {
681 GSL_ERROR("Cholesky decomposition requires square matrix", GSL_ENOTSQR);
682 }
683 else if (N <= CROSSOVER_CHOLESKY)
684 {
685 /* use unblocked Level 2 algorithm */
686 return cholesky_decomp_L2(A);
687 }
688 else
689 {
690 /*
691 * partition matrix:
692 *
693 * A11 A12
694 * A21 A22
695 *
696 * where A11 is N1-by-N1
697 */
698 int status;
699 const size_t N1 = GSL_LINALG_SPLIT(N);
700 const size_t N2 = N - N1;
701 gsl_matrix_view A11 = gsl_matrix_submatrix(A, 0, 0, N1, N1);
702 gsl_matrix_view A21 = gsl_matrix_submatrix(A, N1, 0, N2, N1);
703 gsl_matrix_view A22 = gsl_matrix_submatrix(A, N1, N1, N2, N2);
704
705 /* recursion on A11 */
706 status = cholesky_decomp_L3(&A11.matrix);
707 if (status)
708 return status;
709
710 /* A21 = A21 * L11^{-T} */
711 gsl_blas_dtrsm(CblasRight, CblasLower, CblasTrans, CblasNonUnit, 1.0, &A11.matrix, &A21.matrix);
712
713 /* A22 -= L21 L21^T */
714 gsl_blas_dsyrk(CblasLower, CblasNoTrans, -1.0, &A21.matrix, 1.0, &A22.matrix);
715
716 /* recursion on A22 */
717 status = cholesky_decomp_L3(&A22.matrix);
718 if (status)
719 return status;
720
721 return GSL_SUCCESS;
722 }
723 }
724
725