TooN 2.1
|
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_LAPACK_CHOLESKY_H 00032 #define TOON_INCLUDE_LAPACK_CHOLESKY_H 00033 00034 #include <TooN/TooN.h> 00035 00036 #include <TooN/lapack.h> 00037 00038 #include <assert.h> 00039 00040 namespace TooN { 00041 00042 00043 /** 00044 Decomposes a positive-semidefinite symmetric matrix A (such as a covariance) into L*L^T, where L is lower-triangular. 00045 Also can compute A = S*S^T, with S lower triangular. The LDL^T form is faster to compute than the class Cholesky decomposition. 00046 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. 00047 Also efficiently computes det(A) and rank(A). 00048 It can be used as follows: 00049 @code 00050 // Declare some matrices. 00051 Matrix<3> A = ...; // we'll pretend it is pos-def 00052 Matrix<2,3> M; 00053 Matrix<2> B; 00054 Vector<3> y = make_Vector(2,3,4); 00055 // create the Cholesky decomposition of A 00056 Cholesky<3> chol(A); 00057 // compute x = A^-1 * y 00058 x = cholA.backsub(y); 00059 //compute A^-1 00060 Matrix<3> Ainv = cholA.get_inverse(); 00061 @endcode 00062 @ingroup gDecomps 00063 00064 Cholesky decomposition of a symmetric matrix. 00065 Only the lower half of the matrix is considered 00066 This uses the non-sqrt version of the decomposition 00067 giving symmetric M = L*D*L.T() where the diagonal of L contains ones 00068 @param Size the size of the matrix 00069 @param Precision the precision of the entries in the matrix and its decomposition 00070 **/ 00071 template <int Size, typename Precision=DefaultPrecision> 00072 class Lapack_Cholesky { 00073 public: 00074 00075 Lapack_Cholesky(){} 00076 00077 template<class P2, class B2> 00078 Lapack_Cholesky(const Matrix<Size, Size, P2, B2>& m) 00079 : my_cholesky(m), my_cholesky_lapack(m) { 00080 SizeMismatch<Size,Size>::test(m.num_rows(), m.num_cols()); 00081 do_compute(); 00082 } 00083 00084 /// Constructor for Size=Dynamic 00085 Lapack_Cholesky(int size) : my_cholesky(size,size), my_cholesky_lapack(size,size) {} 00086 00087 template<class P2, class B2> void compute(const Matrix<Size, Size, P2, B2>& m){ 00088 SizeMismatch<Size,Size>::test(m.num_rows(), m.num_cols()); 00089 SizeMismatch<Size,Size>::test(m.num_rows(), my_cholesky.num_rows()); 00090 my_cholesky_lapack=m; 00091 do_compute(); 00092 } 00093 00094 00095 00096 void do_compute(){ 00097 FortranInteger N = my_cholesky.num_rows(); 00098 FortranInteger info; 00099 potrf_("L", &N, my_cholesky_lapack.my_data, &N, &info); 00100 for (int i=0;i<N;i++) { 00101 int j; 00102 for (j=0;j<=i;j++) { 00103 my_cholesky[i][j]=my_cholesky_lapack[j][i]; 00104 } 00105 // LAPACK does not set upper triangle to zero, 00106 // must be done here 00107 for (;j<N;j++) { 00108 my_cholesky[i][j]=0; 00109 } 00110 } 00111 assert(info >= 0); 00112 if (info > 0) { 00113 my_rank = info-1; 00114 } else { 00115 my_rank = N; 00116 } 00117 } 00118 00119 int rank() const { return my_rank; } 00120 00121 template <int Size2, typename P2, typename B2> 00122 Vector<Size, Precision> backsub (const Vector<Size2, P2, B2>& v) const { 00123 SizeMismatch<Size,Size2>::test(my_cholesky.num_cols(), v.size()); 00124 00125 Vector<Size, Precision> result(v); 00126 FortranInteger N=my_cholesky.num_rows(); 00127 FortranInteger NRHS=1; 00128 FortranInteger info; 00129 potrs_("L", &N, &NRHS, my_cholesky_lapack.my_data, &N, result.my_data, &N, &info); 00130 assert(info==0); 00131 return result; 00132 } 00133 00134 template <int Size2, int Cols2, typename P2, typename B2> 00135 Matrix<Size, Cols2, Precision, ColMajor> backsub (const Matrix<Size2, Cols2, P2, B2>& m) const { 00136 SizeMismatch<Size,Size2>::test(my_cholesky.num_cols(), m.num_rows()); 00137 00138 Matrix<Size, Cols2, Precision, ColMajor> result(m); 00139 FortranInteger N=my_cholesky.num_rows(); 00140 FortranInteger NRHS=m.num_cols(); 00141 FortranInteger info; 00142 potrs_("L", &N, &NRHS, my_cholesky_lapack.my_data, &N, result.my_data, &N, &info); 00143 assert(info==0); 00144 return result; 00145 } 00146 00147 template <int Size2, typename P2, typename B2> 00148 Precision mahalanobis(const Vector<Size2, P2, B2>& v) const { 00149 return v * backsub(v); 00150 } 00151 00152 Matrix<Size,Size,Precision> get_L() const { 00153 return my_cholesky; 00154 } 00155 00156 Precision determinant() const { 00157 Precision det = my_cholesky[0][0]; 00158 for (int i=1; i<my_cholesky.num_rows(); i++) 00159 det *= my_cholesky[i][i]; 00160 return det*det; 00161 } 00162 00163 Matrix<> get_inverse() const { 00164 Matrix<Size, Size, Precision> M(my_cholesky.num_rows(),my_cholesky.num_rows()); 00165 M=my_cholesky_lapack; 00166 FortranInteger N = my_cholesky.num_rows(); 00167 FortranInteger info; 00168 potri_("L", &N, M.my_data, &N, &info); 00169 assert(info == 0); 00170 for (int i=1;i<N;i++) { 00171 for (int j=0;j<i;j++) { 00172 M[i][j]=M[j][i]; 00173 } 00174 } 00175 return M; 00176 } 00177 00178 private: 00179 Matrix<Size,Size,Precision> my_cholesky; 00180 Matrix<Size,Size,Precision> my_cholesky_lapack; 00181 FortranInteger my_rank; 00182 }; 00183 00184 00185 } 00186 00187 #endif