Cholesky.h
Go to the documentation of this file.
1 /*!
2  *
3  *
4  * \brief Cholesky Decompositions for a Matrix A = LL^T
5  *
6  *
7  *
8  *
9  * \author O. Krause
10  * \date 2012
11  *
12  *
13  * \par Copyright 1995-2015 Shark Development Team
14  *
15  * <BR><HR>
16  * This file is part of Shark.
17  * <http://image.diku.dk/shark/>
18  *
19  * Shark is free software: you can redistribute it and/or modify
20  * it under the terms of the GNU Lesser General Public License as published
21  * by the Free Software Foundation, either version 3 of the License, or
22  * (at your option) any later version.
23  *
24  * Shark is distributed in the hope that it will be useful,
25  * but WITHOUT ANY WARRANTY; without even the implied warranty of
26  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27  * GNU Lesser General Public License for more details.
28  *
29  * You should have received a copy of the GNU Lesser General Public License
30  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
31  *
32  */
33 
34 #ifndef SHARK_LINALG_CHOLESKY_H
35 #define SHARK_LINALG_CHOLESKY_H
36 
37 #include <shark/LinAlg/Base.h>
39 
40 namespace shark{ namespace blas{
41 
42 /**
43  * \ingroup shark_globals
44  *
45  * @{
46  */
47 
48 /*!
49  * \brief Lower triangular Cholesky decomposition.
50  *
51  * Given an \f$ m \times m \f$ symmetric positive definite matrix
52  * \f$A\f$, compute the lower triangular matrix \f$L\f$ such that
53  * \f$A = LL^T \f$.
54  * An exception is thrown if the matrix is not positive definite.
55  * If you suspect the matrix to be positive semi-definite, use
56  * pivotingCholeskyDecomposition instead
57  *
58  * \param A \f$ m \times m \f$ matrix, which must be symmetric and positive definite
59  * \param L \f$ m \times m \f$ matrix, which stores the Cholesky factor
60  * \return none
61  *
62  *
63  */
64 template<class MatrixT,class MatrixL>
68 ){
69  SIZE_CHECK(A().size1() == A().size2());
70  size_t m = A().size1();
71  ensure_size(L,m, m);
72  L().clear();
73  for(std::size_t i = 0; i != m; ++i){
74  for(std::size_t j = 0; j <= i; ++j){
75  L()(i,j) = A()(i,j);
76  }
77  }
78  if(kernels::potrf<lower>(L()) != 0){
79  throw SHARKEXCEPTION("[Cholesky Decomposition] The Matrix is not positive definite");
80  }
81 }
82 
83 
84 /// \brief Updates a covariance factor by a rank one update
85 ///
86 /// Let \f$ A=LL^T \f$ be a matrix with its lower cholesky factor. Assume we want to update
87 /// A using a simple rank-one update \f$ A = \alpha A+ \beta vv^T \f$. This invalidates L and
88 /// it needs to be recomputed which is O(n^3). instead we can update the factorisation
89 /// directly by performing a similar, albeit more complex algorithm on L, which can be done
90 /// in O(L^2).
91 ///
92 /// Alpha is not required to be positive, but if it is not negative, one has to be carefull
93 /// that the update would keep A positive definite. Otherwise the decomposition does not
94 /// exist anymore and an exception is thrown.
95 ///
96 /// \param L the lower cholesky factor to be updated
97 /// \param v the update vector
98 /// \param alpha the scaling factor, must be positive.
99 /// \param beta the update factor. it Can be positive or negative
100 template<class Matrix,class Vector>
103  vector_expression<Vector> const& v,
104  double alpha, double beta
105 ){
106  //implementation blatantly stolen from Eigen
107  std::size_t n = v().size();
108  blas::vector<double> temp = v();
109  double betaPrime = 1;
110  double a = std::sqrt(alpha);
111  for(std::size_t j=0; j != n; ++j)
112  {
113  double Ljj = a*L()(j,j);
114  double dj = Ljj*Ljj;
115  double wj = temp(j);
116  double swj2 = beta*wj*wj;
117  double gamma = dj*betaPrime + swj2;
118 
119  double x = dj + swj2/betaPrime;
120  if (x <= 0.0)
121  throw SHARKEXCEPTION("[choleskyUpdate] update makes matrix indefinite, no update available");
122  double nLjj = std::sqrt(x);
123  L()(j,j) = nLjj;
124  betaPrime += swj2/dj;
125 
126  // Update the terms of L
127  if(j+1 <n)
128  {
129  subrange(column(L,j),j+1,n) *= a;
130  noalias(subrange(temp,j+1,n)) -= (wj/Ljj) * subrange(column(L,j),j+1,n);
131  if(gamma == 0)
132  continue;
133  subrange(column(L,j),j+1,n) *= nLjj/Ljj;
134  noalias(subrange(column(L,j),j+1,n))+= (nLjj * beta*wj/gamma)*subrange(temp,j+1,n);
135  }
136  }
137 }
138 
139 /*!
140  * \brief Lower triangular Cholesky decomposition with full pivoting performed in place.
141  *
142  * Given an \f$ m \times m \f$ symmetric positive semi-definite matrix
143  * \f$A\f$, compute the lower triangular matrix \f$L\f$ and permutation Matrix P such that
144  * \f$P^TAP = LL^T \f$. If matrix A has rank(A) = k, the first k columns of A hold the full
145  * decomposition, while the rest of the matrix is zero.
146  * This method is slower than the cholesky decomposition without pivoting but numerically more
147  * stable. The diagonal elements are ordered such that i > j => L(i,i) >= L(j,j)
148  *
149  * The implementation used here is described in the working paper
150  * "LAPACK-Style Codes for Level 2 and 3 Pivoted Cholesky Factorizations"
151  * http://www.netlib.org/lapack/lawnspdf/lawn161.pdf
152  *
153  * The computation is carried out in place this means A is destroied and replaced by L.
154  *
155  *
156  * \param Lref \f$ m \times m \f$ matrix, which must be symmetric and positive definite. It is replaced by L in the end.
157  * \param P The pivoting matrix
158  * \return The rank of the matrix A
159  */
160 template<class MatrixL>
164 );
165 
166 /*!
167  * \brief Lower triangular Cholesky decomposition with full pivoting
168  *
169  * Given an \f$ m \times m \f$ symmetric positive semi-definite matrix
170  * \f$A\f$, compute the lower triangular matrix \f$L\f$ and permutation Matrix P such that
171  * \f$P^TAP = LL^T \f$. If matrix A has rank(A) = k, the first k columns of A hold the full
172  * decomposition, while the rest of the matrix is zero.
173  * This method is slower than the cholesky decomposition without pivoting but numerically more
174  * stable. The diagonal elements are ordered such that i > j => L(i,i) >= L(j,j)
175  *
176  * The implementation used here is described in the working paper
177  * "LAPACK-Style Codes for Level 2 and 3 Pivoted Cholesky Factorizations"
178  * http://www.netlib.org/lapack/lawnspdf/lawn161.pdf
179  *
180  *
181  * \param A \f$ m \times m \f$ matrix, which must be symmetric and positive definite
182  * \param P The pivoting matrix
183  * \param L \f$ m \times m \f$ matrix, which stores the Cholesky factor
184  * \return The rank of the matrix A
185  *
186  *
187  */
188 template<class MatrixA,class MatrixL>
193 ){
194  //ensure sizes are correct
195  SIZE_CHECK(A().size1() == A().size2());
196  size_t m = A().size1();
197  ensure_size(L,m,m);
198  noalias(L()) = A;
200 }
201 
202 /** @}*/
203 }}
204 
205 //implementation of the template functions
206 #include "Impl/Cholesky.inl"
207 
208 #endif