cholesky_base.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /* Cholesky Decomposition
24  *
25  * Copyright (C) 2000 Thomas Walter
26  *
27  * 03 May 2000: Modified for GSL by Brian Gough
28  * 29 Jul 2005: Additions by Gerard Jungman
29  * Copyright (C) 2000,2001, 2002, 2003, 2005, 2007 Brian Gough,
30  * Gerard Jungman
31  *
32  * This is free software; you can redistribute it and/or modify it
33  * under the terms of the GNU General Public License as published by the
34  * Free Software Foundation; either version 3, or (at your option) any
35  * later version.
36  *
37  * This source is distributed in the hope that it will be useful, but
38  * WITHOUT ANY WARRANTY; without even the implied warranty of
39  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
40  * General Public License for more details.
41  */
42 /** \file cholesky_base.h
43  \brief File defining Cholesky decomposition
44 */
45 
46 #ifdef DOXYGEN
47 namespace o2scl_linalg {
48 #endif
49 
50  /** \brief Compute the in-place Cholesky decomposition of a symmetric
51  positive-definite square matrix
52 
53  On input, the upper triangular part of A is ignored (only the
54  lower triangular part and diagonal are used). On output,
55  the diagonal and lower triangular part contain the matrix L
56  and the upper triangular part contains L^T.
57 
58  If the matrix is not positive-definite, the error handler
59  will be called, unless \c err_on_fail is false, in which
60  case a non-zero value will be returned.
61  */
62  template<class mat_t> int cholesky_decomp(const size_t M, mat_t &A,
63  bool err_on_fail=true) {
64 
65  size_t i,j,k;
66 
67  /* [GSL] Do the first 2 rows explicitly. It is simple, and faster.
68  And one can return if the matrix has only 1 or 2 rows.
69  */
70 
71  double A_00=O2SCL_IX2(A,0,0);
72 
73  // AWS: The GSL version stores GSL_NAN in L_00 and then throws
74  // an error if A_00 <= 0. We throw the error first and then
75  // the square root should always be safe?
76 
77  if (A_00<=0.0) {
78  if (err_on_fail) {
79  O2SCL_ERR2("Matrix not positive definite (A[0][0]<=0) in ",
80  "cholesky_decomp().",o2scl::exc_einval);
81  } else {
82  return 1;
83  }
84  }
85 
86  double L_00=sqrt(A_00);
87  O2SCL_IX2(A,0,0)=L_00;
88 
89  if (M>1) {
90  double A_10=O2SCL_IX2(A,1,0);
91  double A_11=O2SCL_IX2(A,1,1);
92 
93  double L_10=A_10/L_00;
94  double diag=A_11-L_10*L_10;
95 
96  if (diag<=0.0) {
97  if (err_on_fail) {
98  O2SCL_ERR2("Matrix not positive definite (diag<=0 for 2x2) in ",
99  "cholesky_decomp().",o2scl::exc_einval);
100  } else {
101  return 2;
102  }
103  }
104  double L_11=sqrt(diag);
105 
106  O2SCL_IX2(A,1,0)=L_10;
107  O2SCL_IX2(A,1,1)=L_11;
108  }
109 
110  for (k=2;k<M;k++) {
111  double A_kk=O2SCL_IX2(A,k,k);
112 
113  for (i=0;i<k;i++) {
114  double sum=0.0;
115 
116  double A_ki=O2SCL_IX2(A,k,i);
117  double A_ii=O2SCL_IX2(A,i,i);
118 
119  // AWS: Should change to use a form of ddot() here
120  if (i>0) {
121  sum=0.0;
122  for(j=0;j<i;j++) {
123  sum+=O2SCL_IX2(A,i,j)*O2SCL_IX2(A,k,j);
124  }
125  }
126 
127  A_ki=(A_ki-sum)/A_ii;
128  O2SCL_IX2(A,k,i)=A_ki;
129  }
130 
131  {
132  double sum=dnrm2_subrow(A,k,0,k);
133  double diag=A_kk-sum*sum;
134 
135  if (diag<=0.0) {
136  if (err_on_fail) {
137  O2SCL_ERR2("Matrix not positive definite (diag<=0) in ",
138  "cholesky_decomp().",o2scl::exc_einval);
139  } else {
140  return 3;
141  }
142  }
143 
144  double L_kk=sqrt(diag);
145 
146  O2SCL_IX2(A,k,k)=L_kk;
147  }
148  }
149 
150  /* [GSL] Now copy the transposed lower triangle to the upper
151  triangle, the diagonal is common.
152  */
153 
154  for (i=1;i<M;i++) {
155  for (j=0;j<i;j++) {
156  double A_ij=O2SCL_IX2(A,i,j);
157  O2SCL_IX2(A,j,i)=A_ij;
158  }
159  }
160 
161  return 0;
162  }
163 
164  /** \brief Compute the determinant of a matrix from its Cholesky decomposition
165 
166  */
167  template<class mat_t>
168  double cholesky_det(const size_t M, const mat_t &A) {
169 
170  double det=1.0;
171 
172  for (size_t i=0;i<M;i++) {
173  det*=O2SCL_IX2(A,i,i);
174  }
175 
176  return det;
177  }
178 
179  /** \brief Compute the logarithm of the absolute value of the
180  determinant of a matrix from its Cholesky decomposition
181  */
182  template<class mat_t>
183  double cholesky_lndet(const size_t M, const mat_t &A) {
184 
185  double lndet = 0.0;
186 
187  for (size_t i = 0; i < M; i++) {
188  lndet+=log(fabs(O2SCL_IX2(A,i,i)));
189  }
190 
191  return lndet;
192  }
193 
194  /** \brief Solve a symmetric positive-definite linear system after a
195  Cholesky decomposition
196 
197  Given the Cholesky decomposition of a matrix A in \c LLT,
198  this function solves the system <tt>A*x=b</tt>.
199  */
200  template<class mat_t, class vec_t, class vec2_t>
201  void cholesky_solve(const size_t N, const mat_t &LLT,
202  const vec_t &b, vec2_t &x) {
203 
204  // [GSL] Copy x <- b
205  o2scl::vector_copy(N,b,x);
206 
207  // [GSL] Solve for c using forward-substitution, L c=b
208  o2scl_cblas::dtrsv(o2scl_cblas::o2cblas_RowMajor,
209  o2scl_cblas::o2cblas_Lower,
210  o2scl_cblas::o2cblas_NoTrans,
211  o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
212 
213  // [GSL] Perform back-substitution,U x=c
214  o2scl_cblas::dtrsv(o2scl_cblas::o2cblas_RowMajor,
215  o2scl_cblas::o2cblas_Upper,
216  o2scl_cblas::o2cblas_NoTrans,
217  o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
218 
219  return;
220  }
221 
222  /** \brief Solve a linear system in place using a Cholesky decomposition
223  */
224  template<class mat_t, class vec_t>
225  void cholesky_svx(const size_t N, const mat_t &LLT, vec_t &x) {
226 
227  // [GSL] Solve for c using forward-substitution, L c=b
228  o2scl_cblas::dtrsv(o2scl_cblas::o2cblas_RowMajor,
229  o2scl_cblas::o2cblas_Lower,
230  o2scl_cblas::o2cblas_NoTrans,
231  o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
232 
233  // [GSL] Perform back-substitution, U x=c
234  o2scl_cblas::dtrsv(o2scl_cblas::o2cblas_RowMajor,
235  o2scl_cblas::o2cblas_Upper,
236  o2scl_cblas::o2cblas_NoTrans,
237  o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
238 
239  return;
240  }
241 
242  /** \brief Compute the inverse of a symmetric positive definite matrix
243  given the Cholesky decomposition
244 
245  Given a Cholesky decomposition produced by \ref cholesky_decomp(),
246  this function returns the inverse of that matrix in \c LLT.
247  */
248  template<class mat_t> void cholesky_invert(const size_t N, mat_t &LLT) {
249 
250  size_t i, j;
251  double sum;
252 
253  // [GSL] invert the lower triangle of LLT
254  for (i=0;i<N;++i) {
255 
256  j=N-i-1;
257 
258  O2SCL_IX2(LLT,j,j)=1.0/O2SCL_IX2(LLT,j,j);
259  double ajj=-O2SCL_IX2(LLT,j,j);
260 
261  if (j<N-1) {
262 
263  // This section is just the equivalent of dtrmv() for a part of
264  // the matrix LLT.
265  {
266 
267  size_t ix=N-j-2;
268  for (size_t ii=N-j-1;ii>0 && ii--;) {
269  double temp=0.0;
270  const size_t j_min=0;
271  const size_t j_max=ii;
272  size_t jx=j_min;
273  for (size_t jj=j_min;jj<j_max;jj++) {
274  temp+=O2SCL_IX2(LLT,jx+j+1,j)*O2SCL_IX2(LLT,ii+j+1,jj+j+1);
275  jx++;
276  }
277  O2SCL_IX2(LLT,ix+j+1,j)=temp+O2SCL_IX2(LLT,ix+j+1,j)*
278  O2SCL_IX2(LLT,ii+j+1,ii+j+1);
279  ix--;
280  }
281 
282  }
283 
284  o2scl_cblas::dscal_subcol(LLT,j+1,j,N,ajj);
285 
286  }
287  }
288 
289  /*
290  [GSL] The lower triangle of LLT now contains L^{-1}. Now compute
291  A^{-1}=L^{-t} L^{-1}
292 
293  The (ij) element of A^{-1} is column i of L^{-1} dotted into
294  column j of L^{-1}
295  */
296 
297  for (i=0;i<N;++i) {
298 
299  for (j=i+1;j<N;++j) {
300 
301  // [GSL] Compute Ainv_{ij}=sum_k Linv_{ki} Linv_{kj}.
302 
303  // AWS: Should change to use a form of ddot() here
304  sum=0.0;
305  for(size_t k=j;k<N;k++) {
306  sum+=O2SCL_IX2(LLT,k,i)*O2SCL_IX2(LLT,k,j);
307  }
308 
309  // [GSL] Store in upper triangle
310  O2SCL_IX2(LLT,i,j)=sum;
311  }
312 
313  // [GSL] now compute the diagonal element
314 
315  // AWS: Should change to use a form of ddot() here
316  sum=0.0;
317  for(size_t k=i;k<N;k++) {
318  sum+=O2SCL_IX2(LLT,k,i)*O2SCL_IX2(LLT,k,i);
319  }
320 
321  O2SCL_IX2(LLT,i,i)=sum;
322  }
323 
324  // [GSL] Copy the transposed upper triangle to the lower triangle
325 
326  for (j=1;j<N;j++) {
327  for (i=0;i<j;i++) {
328  O2SCL_IX2(LLT,j,i)=O2SCL_IX2(LLT,i,j);
329  }
330  }
331 
332  return;
333  }
334 
335  /** \brief Cholesky decomposition with unit-diagonal triangular parts.
336 
337  A = L D L^T, where diag(L) = (1,1,...,1).
338  Upon exit, A contains L and L^T as for Cholesky, and
339  the diagonal of A is (1,1,...,1). The vector Dis set
340  to the diagonal elements of the diagonal matrix D.
341  */
342  template<class mat_t, class vec_t>
343  int cholesky_decomp_unit(const size_t N, mat_t &A, vec_t &D) {
344 
345  size_t i, j;
346 
347  // [GSL] Initial Cholesky decomposition
348  int stat_chol=cholesky_decomp(N,A);
349 
350  // [GSL] Calculate D from diagonal part of initial Cholesky
351  for(i=0;i<N;++i) {
352  const double C_ii=O2SCL_IX2(A,i,i);
353  O2SCL_IX(D,i)=C_ii*C_ii;
354  }
355 
356  // [GSL] Multiply initial Cholesky by 1/sqrt(D) on the right
357  for(i=0;i<N;++i) {
358  for(j=0;j<N;++j) {
359  O2SCL_IX2(A,i,j)=O2SCL_IX2(A,i,j)/sqrt(O2SCL_IX(D,j));
360  }
361  }
362 
363  /* [GSL] Because the initial Cholesky contained both L and
364  transpose(L), the result of the multiplication is not symmetric
365  anymore; but the lower triangle _is_ correct. Therefore we
366  reflect it to the upper triangle and declare victory.
367  */
368  for(i=0;i<N;++i) {
369  for(j=i+1;j<N;++j) {
370  O2SCL_IX2(A,i,j)=O2SCL_IX2(A,j,i);
371  }
372  }
373 
374  return stat_chol;
375  }
376 
377 #ifdef DOXYGEN
378 }
379 #endif
o2scl_cblas::dnrm2_subrow
double dnrm2_subrow(const mat_t &M, const size_t ir, const size_t ic, const size_t N)
Compute the norm of a subrow of a matrix.
Definition: cblas_base.h:1568
o2scl_linalg
The namespace for linear algebra classes and functions.
Definition: bidiag.h:36
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
o2scl_cblas::dscal_subcol
void dscal_subcol(mat_t &A, const size_t ir, const size_t ic, const size_t M, const double alpha)
Compute for a subcolumn of a matrix.
Definition: cblas_base.h:1399
o2scl::vector_copy
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:124
o2scl_linalg::cholesky_svx
void cholesky_svx(const size_t N, const mat_t &LLT, vec_t &x)
Solve a linear system in place using a Cholesky decomposition.
Definition: cholesky_base.h:225
o2scl_linalg::cholesky_det
double cholesky_det(const size_t M, const mat_t &A)
Compute the determinant of a matrix from its Cholesky decomposition.
Definition: cholesky_base.h:168
o2scl_linalg::cholesky_decomp_unit
int cholesky_decomp_unit(const size_t N, mat_t &A, vec_t &D)
Cholesky decomposition with unit-diagonal triangular parts.
Definition: cholesky_base.h:343
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl_linalg::cholesky_decomp
int cholesky_decomp(const size_t M, mat_t &A, bool err_on_fail=true)
Compute the in-place Cholesky decomposition of a symmetric positive-definite square matrix.
Definition: cholesky_base.h:62
o2scl_linalg::cholesky_lndet
double cholesky_lndet(const size_t M, const mat_t &A)
Compute the logarithm of the absolute value of the determinant of a matrix from its Cholesky decompos...
Definition: cholesky_base.h:183
o2scl_cblas::dtrsv
void dtrsv(const enum o2cblas_order order, const enum o2cblas_uplo Uplo, const enum o2cblas_transpose TransA, const enum o2cblas_diag Diag, const size_t M, const size_t N, const mat_t &A, vec_t &X)
Compute .
Definition: cblas_base.h:335
o2scl_linalg::cholesky_solve
void cholesky_solve(const size_t N, const mat_t &LLT, const vec_t &b, vec2_t &x)
Solve a symmetric positive-definite linear system after a Cholesky decomposition.
Definition: cholesky_base.h:201
o2scl_linalg::cholesky_invert
void cholesky_invert(const size_t N, mat_t &LLT)
Compute the inverse of a symmetric positive definite matrix given the Cholesky decomposition.
Definition: cholesky_base.h:248

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).