TooN 2.1
Cholesky.h
00001 // -*- c++ -*-
00002 
00003 //     Copyright (C) 2009 Tom Drummond (twd20@cam.ac.uk)
00004 //
00005 // This file is part of the TooN Library.  This library is free
00006 // software; you can redistribute it and/or modify it under the
00007 // terms of the GNU General Public License as published by the
00008 // Free Software Foundation; either version 2, or (at your option)
00009 // any later version.
00010 
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014 // GNU General Public License for more details.
00015 
00016 // You should have received a copy of the GNU General Public License along
00017 // with this library; see the file COPYING.  If not, write to the Free
00018 // Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
00019 // USA.
00020 
00021 // As a special exception, you may use this file as part of a free software
00022 // library without restriction.  Specifically, if other files instantiate
00023 // templates or use macros or inline functions from this file, or you compile
00024 // this file and link it with other files to produce an executable, this
00025 // file does not by itself cause the resulting executable to be covered by
00026 // the GNU General Public License.  This exception does not however
00027 // invalidate any other reasons why the executable file might be covered by
00028 // the GNU General Public License.
00029 
00030 
00031 #ifndef TOON_INCLUDE_CHOLESKY_H
00032 #define TOON_INCLUDE_CHOLESKY_H
00033 
00034 #include <TooN/TooN.h>
00035 
00036 namespace TooN {
00037 
00038 
00039 /**
00040 Decomposes a positive-semidefinite symmetric matrix A (such as a covariance) into L*D*L^T, where L is lower-triangular and D is diagonal.
00041 Also can compute the classic A = L*L^T, with L lower triangular.  The LDL^T form is faster to compute than the classical Cholesky decomposition. 
00042 Use get_unscaled_L() and get_D() to access the individual matrices of L*D*L^T decomposition. Use get_L() to access the lower triangular matrix of the classic Cholesky decomposition L*L^T.
00043 The decomposition can be used to compute A^-1*x, A^-1*M, M*A^-1*M^T, and A^-1 itself, though the latter rarely needs to be explicitly represented.
00044 Also efficiently computes det(A) and rank(A).
00045 It can be used as follows:
00046 @code
00047 // Declare some matrices.
00048 Matrix<3> A = ...; // we'll pretend it is pos-def
00049 Matrix<2,3> M;
00050 Matrix<2> B;
00051 Vector<3> y = make_Vector(2,3,4);
00052 // create the Cholesky decomposition of A
00053 Cholesky<3> chol(A);
00054 // compute x = A^-1 * y
00055 x = cholA.backsub(y);
00056 //compute A^-1
00057 Matrix<3> Ainv = cholA.get_inverse();
00058 @endcode
00059 @ingroup gDecomps
00060 
00061 Cholesky decomposition of a symmetric matrix.
00062 Only the lower half of the matrix is considered
00063 This uses the non-sqrt version of the decomposition
00064 giving symmetric M = L*D*L.T() where the diagonal of L contains ones
00065 @param Size the size of the matrix
00066 @param Precision the precision of the entries in the matrix and its decomposition
00067 **/
00068 template <int Size=Dynamic, class Precision=DefaultPrecision>
00069 class Cholesky {
00070 public:
00071     Cholesky(){}
00072 
00073     /// Construct the Cholesky decomposition of a matrix. This initialises the class, and
00074     /// performs the decomposition immediately.
00075     /// Run time is O(N^3)
00076     template<class P2, class B2>
00077     Cholesky(const Matrix<Size, Size, P2, B2>& m)
00078         : my_cholesky(m) {
00079         SizeMismatch<Size,Size>::test(m.num_rows(), m.num_cols());
00080         do_compute();
00081     }
00082     
00083     /// Constructor for Size=Dynamic
00084     Cholesky(int size) : my_cholesky(size,size) {}
00085 
00086 
00087     /// Compute the LDL^T decomposition of another matrix.
00088     /// Run time is O(N^3)
00089     template<class P2, class B2> void compute(const Matrix<Size, Size, P2, B2>& m){
00090         SizeMismatch<Size,Size>::test(m.num_rows(), m.num_cols());
00091         SizeMismatch<Size,Size>::test(m.num_rows(), my_cholesky.num_rows());
00092         my_cholesky=m;
00093         do_compute();
00094     }
00095     
00096     private:
00097     void do_compute() {
00098         int size=my_cholesky.num_rows();
00099         for(int col=0; col<size; col++){
00100             Precision inv_diag = 1;
00101             for(int row=col; row < size; row++){
00102                 // correct for the parts of cholesky already computed
00103                 Precision val = my_cholesky(row,col);
00104                 for(int col2=0; col2<col; col2++){
00105                     // val-=my_cholesky(col,col2)*my_cholesky(row,col2)*my_cholesky(col2,col2);
00106                     val-=my_cholesky(col2,col)*my_cholesky(row,col2);
00107                 }
00108                 if(row==col){
00109                     // this is the diagonal element so don't divide
00110                     my_cholesky(row,col)=val;
00111                     if(val == 0){
00112                         my_rank = row;
00113                         return;
00114                     }
00115                     inv_diag=1/val;
00116                 } else {
00117                     // cache the value without division in the upper half
00118                     my_cholesky(col,row)=val;
00119                     // divide my the diagonal element for all others
00120                     my_cholesky(row,col)=val*inv_diag;
00121                 }
00122             }
00123         }
00124         my_rank = size;
00125     }
00126 
00127     public:
00128 
00129     /// Compute x = A^-1*v
00130     /// Run time is O(N^2)
00131     template<int Size2, class P2, class B2>
00132     Vector<Size, Precision> backsub (const Vector<Size2, P2, B2>& v) const {
00133         int size=my_cholesky.num_rows();
00134         SizeMismatch<Size,Size2>::test(size, v.size());
00135 
00136         // first backsub through L
00137         Vector<Size, Precision> y(size);
00138         for(int i=0; i<size; i++){
00139             Precision val = v[i];
00140             for(int j=0; j<i; j++){
00141                 val -= my_cholesky(i,j)*y[j];
00142             }
00143             y[i]=val;
00144         }
00145         
00146         // backsub through diagonal
00147         for(int i=0; i<size; i++){
00148             y[i]/=my_cholesky(i,i);
00149         }
00150 
00151         // backsub through L.T()
00152         Vector<Size,Precision> result(size);
00153         for(int i=size-1; i>=0; i--){
00154             Precision val = y[i];
00155             for(int j=i+1; j<size; j++){
00156                 val -= my_cholesky(j,i)*result[j];
00157             }
00158             result[i]=val;
00159         }
00160 
00161         return result;
00162     }
00163 
00164     /**overload
00165     */
00166     template<int Size2, int C2, class P2, class B2>
00167     Matrix<Size, C2, Precision> backsub (const Matrix<Size2, C2, P2, B2>& m) const {
00168         int size=my_cholesky.num_rows();
00169         SizeMismatch<Size,Size2>::test(size, m.num_rows());
00170 
00171         // first backsub through L
00172         Matrix<Size, C2, Precision> y(size, m.num_cols());
00173         for(int i=0; i<size; i++){
00174             Vector<C2, Precision> val = m[i];
00175             for(int j=0; j<i; j++){
00176                 val -= my_cholesky(i,j)*y[j];
00177             }
00178             y[i]=val;
00179         }
00180         
00181         // backsub through diagonal
00182         for(int i=0; i<size; i++){
00183             y[i]*=(1/my_cholesky(i,i));
00184         }
00185 
00186         // backsub through L.T()
00187         Matrix<Size,C2,Precision> result(size, m.num_cols());
00188         for(int i=size-1; i>=0; i--){
00189             Vector<C2,Precision> val = y[i];
00190             for(int j=i+1; j<size; j++){
00191                 val -= my_cholesky(j,i)*result[j];
00192             }
00193             result[i]=val;
00194         }
00195         return result;
00196     }
00197 
00198 
00199     /// Compute A^-1 and store in M
00200     /// Run time is O(N^3)
00201     // easy way to get inverse - could be made more efficient
00202     Matrix<Size,Size,Precision> get_inverse(){
00203         Matrix<Size,Size,Precision>I(Identity(my_cholesky.num_rows()));
00204         return backsub(I);
00205     }
00206     
00207     ///Compute the determinant.
00208     Precision determinant(){
00209         Precision answer=my_cholesky(0,0);
00210         for(int i=1; i<my_cholesky.num_rows(); i++){
00211             answer*=my_cholesky(i,i);
00212         }
00213         return answer;
00214     }
00215 
00216     template <int Size2, typename P2, typename B2>
00217     Precision mahalanobis(const Vector<Size2, P2, B2>& v) const {
00218         return v * backsub(v);
00219     }
00220 
00221     Matrix<Size,Size,Precision> get_unscaled_L() const {
00222         Matrix<Size,Size,Precision> m(my_cholesky.num_rows(),
00223                           my_cholesky.num_rows());
00224         m=Identity;
00225         for (int i=1;i<my_cholesky.num_rows();i++) {
00226             for (int j=0;j<i;j++) {
00227                 m(i,j)=my_cholesky(i,j);
00228             }
00229         }
00230         return m;
00231     }
00232             
00233     Matrix<Size,Size,Precision> get_D() const {
00234         Matrix<Size,Size,Precision> m(my_cholesky.num_rows(),
00235                           my_cholesky.num_rows());
00236         m=Zeros;
00237         for (int i=0;i<my_cholesky.num_rows();i++) {
00238             m(i,i)=my_cholesky(i,i);
00239         }
00240         return m;
00241     }
00242     
00243     Matrix<Size,Size,Precision> get_L() const {
00244         using std::sqrt;
00245         Matrix<Size,Size,Precision> m(my_cholesky.num_rows(),
00246                           my_cholesky.num_rows());
00247         m=Zeros;
00248         for (int j=0;j<my_cholesky.num_cols();j++) {
00249             Precision sqrtd=sqrt(my_cholesky(j,j));
00250             m(j,j)=sqrtd;
00251             for (int i=j+1;i<my_cholesky.num_rows();i++) {
00252                 m(i,j)=my_cholesky(i,j)*sqrtd;
00253             }
00254         }
00255         return m;
00256     }
00257 
00258     int rank() const { return my_rank; }
00259 
00260 private:
00261     Matrix<Size,Size,Precision> my_cholesky;
00262     int my_rank;
00263 };
00264 
00265 
00266 }
00267 
00268 #endif