v0.9.0
cholesky.hpp
Go to the documentation of this file.
1 /** -*- c++ -*- \file cholesky.hpp \brief cholesky decomposition */
2 /*
3  - begin : 2005-08-24
4  - copyright : (C) 2005 by Gunter Winkler, Konstantin Kutzkow
5  - email : guwi17@gmx.de
6 
7  This library is free software; you can redistribute it and/or
8  modify it under the terms of the GNU Lesser General Public
9  License as published by the Free Software Foundation; either
10  version 2.1 of the License, or (at your option) any later version.
11 
12  This library is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15  Lesser General Public License for more details.
16 
17  You should have received a copy of the GNU Lesser General Public
18  License along with this library; if not, write to the Free Software
19  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
20 
21 */
22 
23 #ifndef _H_CHOLESKY_HPP_
24 #define _H_CHOLESKY_HPP_
25 
26 
27 #include <cassert>
28 
29 #include <boost/numeric/ublas/vector.hpp>
30 #include <boost/numeric/ublas/vector_proxy.hpp>
31 
32 #include <boost/numeric/ublas/matrix.hpp>
33 #include <boost/numeric/ublas/matrix_proxy.hpp>
34 
35 #include <boost/numeric/ublas/vector_expression.hpp>
36 #include <boost/numeric/ublas/matrix_expression.hpp>
37 
38 #include <boost/numeric/ublas/triangular.hpp>
39 
40 //namespace ublas = boost::numeric::ublas;
41 
42 
43 /** \brief decompose the symmetric positive definit matrix A into product L L^T.
44  *
45  * \param MATRIX type of input matrix
46  * \param TRIA type of lower triangular output matrix
47  * \param A square symmetric positive definite input matrix (only the lower triangle is accessed)
48  * \param L lower triangular output matrix
49  * \return nonzero if decompositon fails (the value ist 1 + the numer of the failing row)
50  */
51 template < class MATRIX, class TRIA >
52 size_t cholesky_decompose(const MATRIX& A, TRIA& L)
53 {
54  using namespace ublas;
55 
56  //typedef typename MATRIX::value_type T;
57 
58  assert( A.size1() == A.size2() );
59  assert( A.size1() == L.size1() );
60  assert( A.size2() == L.size2() );
61 
62  const size_t n = A.size1();
63 
64  for (size_t k=0 ; k < n; k++) {
65 
66  double qL_kk = A(k,k) - inner_prod( project( row(L, k), range(0, k) ),
67  project( row(L, k), range(0, k) ) );
68 
69  if (qL_kk <= 0) {
70  return 1 + k;
71  } else {
72  double L_kk = sqrt( qL_kk );
73  L(k,k) = L_kk;
74 
75  matrix_column<TRIA> cLk(L, k);
76  project( cLk, range(k+1, n) )
77  = ( project( column(A, k), range(k+1, n) )
78  - prod( project(L, range(k+1, n), range(0, k)),
79  project(row(L, k), range(0, k) ) ) ) / L_kk;
80  }
81  }
82  return 0;
83 }
84 
85 
86 /** \brief decompose the symmetric positive definit matrix A into product L L^T.
87  *
88  * \param MATRIX type of matrix A
89  * \param A input: square symmetric positive definite matrix (only the lower triangle is accessed)
90  * \param A output: the lower triangle of A is replaced by the cholesky factor
91  * \return nonzero if decompositon fails (the value ist 1 + the numer of the failing row)
92  */
93 template < class MATRIX >
94 size_t cholesky_decompose(MATRIX& A)
95 {
96  using namespace ublas;
97 
98  //typedef typename MATRIX::value_type T;
99 
100  const MATRIX& A_c(A);
101 
102  const size_t n = A.size1();
103 
104  for (size_t k=0 ; k < n; k++) {
105 
106  double qL_kk = A_c(k,k) - inner_prod( project( row(A_c, k), range(0, k) ),
107  project( row(A_c, k), range(0, k) ) );
108 
109  if (qL_kk <= 0) {
110  return 1 + k;
111  } else {
112  double L_kk = sqrt( qL_kk );
113 
114  matrix_column<MATRIX> cLk(A, k);
115  project( cLk, range(k+1, n) )
116  = ( project( column(A_c, k), range(k+1, n) )
117  - prod( project(A_c, range(k+1, n), range(0, k)),
118  project(row(A_c, k), range(0, k) ) ) ) / L_kk;
119  A(k,k) = L_kk;
120  }
121  }
122  return 0;
123 }
124 
125 #if 0
126  using namespace ublas;
127 
128  // Operations:
129  // n * (n - 1) / 2 + n = n * (n + 1) / 2 multiplications,
130  // n * (n - 1) / 2 additions
131 
132  // Dense (proxy) case
133  template<class E1, class E2>
134  void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
135  lower_tag, column_major_tag) {
136  std::cout << " is_lc ";
137  typedef typename E2::size_type size_type;
138  typedef typename E2::difference_type difference_type;
139  typedef typename E2::value_type value_type;
140 
141  BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
142  BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
143  size_type size = e2 ().size ();
144  for (size_type n = 0; n < size; ++ n) {
145 #ifndef BOOST_UBLAS_SINGULAR_CHECK
146  BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
147 #else
148  if (e1 () (n, n) == value_type/*zero*/())
149  singular ().raise ();
150 #endif
151  value_type t = e2 () (n) / e1 () (n, n);
152  e2 () (n) = t;
153  if (t != value_type/*zero*/()) {
154  project( e2 (), range(n+1, size) )
155  .minus_assign( t * project( column( e1 (), n), range(n+1, size) ) );
156  }
157  }
158  }
159 #endif
160 
161 
162 
163 
164 
165 /** \brief decompose the symmetric positive definit matrix A into product L L^T.
166  *
167  * \param MATRIX type of matrix A
168  * \param A input: square symmetric positive definite matrix (only the lower triangle is accessed)
169  * \param A output: the lower triangle of A is replaced by the cholesky factor
170  * \return nonzero if decompositon fails (the value ist 1 + the numer of the failing row)
171  */
172 template < class MATRIX >
174 {
175  using namespace ublas;
176 
177  typedef typename MATRIX::value_type T;
178 
179  // read access to a const matrix is faster
180  const MATRIX& A_c(A);
181 
182  const size_t n = A.size1();
183 
184  for (size_t k=0 ; k < n; k++) {
185 
186  double qL_kk = A_c(k,k) - inner_prod( project( row( A_c, k ), range(0, k) ),
187  project( row( A_c, k ), range(0, k) ) );
188 
189  if (qL_kk <= 0) {
190  return 1 + k;
191  } else {
192  double L_kk = sqrt( qL_kk );
193 
194  // aktualisieren
195  for (size_t i = k+1; i < A.size1(); ++i) {
196  T* Aik = A.find_element(i, k);
197 
198  if (Aik != 0) {
199  *Aik = ( *Aik - inner_prod( project( row( A_c, k ), range(0, k) ),
200  project( row( A_c, i ), range(0, k) ) ) ) / L_kk;
201  }
202  }
203 
204  A(k,k) = L_kk;
205  }
206  }
207 
208  return 0;
209 }
210 
211 
212 
213 
214 /** \brief solve system L L^T x = b inplace
215  *
216  * \param L a triangular matrix
217  * \param x input: right hand side b; output: solution x
218  */
219 template < class TRIA, class VEC >
220 void
221 cholesky_solve(const TRIA& L, VEC& x, ublas::lower)
222 {
223  using namespace ublas;
224 // ::inplace_solve(L, x, lower_tag(), typename TRIA::orientation_category () );
225  inplace_solve(L, x, lower_tag() );
226  inplace_solve(trans(L), x, upper_tag());
227 }
228 
229 
230 #endif
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
Definition: cholesky.hpp:52
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
size_t incomplete_cholesky_decompose(MATRIX &A)
decompose the symmetric positive definit matrix A into product L L^T.
Definition: cholesky.hpp:173