00001 // -*- c++ -*- 00002 00003 // Copyright (C) 2005,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 #ifndef __SYMEIGEN_H 00031 #define __SYMEIGEN_H 00032 00033 #include <iostream> 00034 #include <cassert> 00035 #include <cmath> 00036 #include <TooN/lapack.h> 00037 00038 #include <TooN/TooN.h> 00039 00040 namespace TooN { 00041 00042 00043 namespace Internal{ 00044 ///Default condition number for SymEigen::backsub, SymEigen::get_pinv and SymEigen::get_inv_diag 00045 static const double symeigen_condition_no=1e9; 00046 00047 ///@internal 00048 ///@brief Compute eigensystems for sizes > 2 00049 ///Helper struct for computing eigensystems, to allow for specialization on 00050 ///2x2 matrices. 00051 ///@ingroup gInternal 00052 template <int Size> struct ComputeSymEigen { 00053 00054 ///@internal 00055 ///Compute an eigensystem. 00056 ///@param m Input matrix (assumed to be symmetric) 00057 ///@param evectors Eigen vector output 00058 ///@param evalues Eigen values output 00059 template<typename P, typename B> 00060 static inline void compute(const Matrix<Size,Size,P, B>& m, Matrix<Size,Size,P> & evectors, Vector<Size, P>& evalues) { 00061 evectors = m; 00062 int N = evalues.size(); 00063 int lda = evalues.size(); 00064 int info; 00065 int lwork=-1; 00066 P size; 00067 00068 // find out how much space fortran needs 00069 syev_((char*)"V",(char*)"U",&N,&evectors[0][0],&lda,&evalues[0], &size,&lwork,&info); 00070 lwork = int(size); 00071 Vector<Dynamic, P> WORK(lwork); 00072 00073 // now compute the decomposition 00074 syev_((char*)"V",(char*)"U",&N,&evectors[0][0],&lda,&evalues[0], &WORK[0],&lwork,&info); 00075 00076 if(info!=0){ 00077 std::cerr << "In SymEigen<"<<Size<<">: " << info 00078 << " off-diagonal elements of an intermediate tridiagonal form did not converge to zero." << std::endl 00079 << "M = " << m << std::endl; 00080 } 00081 } 00082 }; 00083 00084 ///@internal 00085 ///@brief Compute 2x2 eigensystems 00086 ///Helper struct for computing eigensystems, specialized on 2x2 matrices. 00087 ///@ingroup gInternal 00088 template <> struct ComputeSymEigen<2> { 00089 00090 ///@internal 00091 ///Compute an eigensystem. 00092 ///@param m Input matrix (assumed to be symmetric) 00093 ///@param eig Eigen vector output 00094 ///@param ev Eigen values output 00095 template<typename P, typename B> 00096 static inline void compute(const Matrix<2,2,P,B>& m, Matrix<2,2,P>& eig, Vector<2, P>& ev) { 00097 double trace = m[0][0] + m[1][1]; 00098 double det = m[0][0]*m[1][1] - m[0][1]*m[1][0]; 00099 double disc = trace*trace - 4 * det; 00100 assert(disc>=0); 00101 using std::sqrt; 00102 double root_disc = sqrt(disc); 00103 ev[0] = 0.5 * (trace - root_disc); 00104 ev[1] = 0.5 * (trace + root_disc); 00105 double a = m[0][0] - ev[0]; 00106 double b = m[0][1]; 00107 double magsq = a*a + b*b; 00108 if (magsq == 0) { 00109 eig[0][0] = 1.0; 00110 eig[0][1] = 0; 00111 } else { 00112 eig[0][0] = -b; 00113 eig[0][1] = a; 00114 eig[0] *= 1.0/sqrt(magsq); 00115 } 00116 eig[1][0] = -eig[0][1]; 00117 eig[1][1] = eig[0][0]; 00118 } 00119 }; 00120 00121 }; 00122 00123 /** 00124 Performs eigen decomposition of a matrix. 00125 Real symmetric (and hence square matrices) can be decomposed into 00126 \f[M = U \times \Lambda \times U^T\f] 00127 where \f$U\f$ is an orthogonal matrix (and hence \f$U^T = U^{-1}\f$) whose columns 00128 are the eigenvectors of \f$M\f$ and \f$\Lambda\f$ is a diagonal matrix whose entries 00129 are the eigenvalues of \f$M\f$. These quantities are often of use directly, and can 00130 be obtained as follows: 00131 @code 00132 // construct M 00133 Matrix<3> M(3,3); 00134 M[0]=makeVector(4,0,2); 00135 M[1]=makeVector(0,5,3); 00136 M[2]=makeVector(2,3,6); 00137 00138 // create the eigen decomposition of M 00139 SymEigen<3> eigM(M); 00140 cout << "A=" << M << endl; 00141 cout << "(E,v)=eig(A)" << endl; 00142 // print the smallest eigenvalue 00143 cout << "v[0]=" << eigM.get_evalues()[0] << endl; 00144 // print the associated eigenvector 00145 cout << "E[0]=" << eigM.get_evectors()[0] << endl; 00146 @endcode 00147 00148 Further, provided the eigenvalues are nonnegative, the square root of 00149 a matrix and its inverse can also be obtained, 00150 @code 00151 // print the square root of the matrix. 00152 cout << "R=sqrtm(A)=" << eigM.get_sqrtm() << endl; 00153 // print the square root of the matrix squared. 00154 cout << "(should equal A), R^T*R=" 00155 << eigM.get_sqrtm().T() * eigM.get_sqrtm() << endl; 00156 // print the inverse of the matrix. 00157 cout << "A^-1=" << eigM.get_pinv() << endl; 00158 // print the inverse square root of the matrix. 00159 cout << "C=isqrtm(A)=" << eigM.get_isqrtm() << endl; 00160 // print the inverse square root of the matrix squared. 00161 cout << "(should equal A^-1), C^T*C=" 00162 << eigM.get_isqrtm().T() * eigM.get_isqrtm() << endl; 00163 @endcode 00164 00165 This decomposition is very similar to the SVD (q.v.), and can be used to solve 00166 equations using backsub() or get_pinv(), with the same treatment of condition numbers. 00167 00168 SymEigen<> (= SymEigen<-1>) can be used to create an eigen decomposition whose size is determined at run-time. 00169 @ingroup gDecomps 00170 **/ 00171 template <int Size=Dynamic, typename Precision = double> 00172 class SymEigen { 00173 public: 00174 inline SymEigen(){} 00175 00176 /// Initialise this eigen decomposition but do no immediately 00177 /// perform a decomposition. 00178 /// 00179 /// @param m The size of the matrix to perform the eigen decomposition on. 00180 inline SymEigen(int m) : my_evectors(m,m), my_evalues(m) {} 00181 00182 /// Construct the eigen decomposition of a matrix. This initialises the class, and 00183 /// performs the decomposition immediately. 00184 template<int R, int C, typename B> 00185 inline SymEigen(const Matrix<R, C, Precision, B>& m) : my_evectors(m.num_rows(), m.num_cols()), my_evalues(m.num_rows()) { 00186 compute(m); 00187 } 00188 00189 /// Perform the eigen decomposition of a matrix. 00190 template<int R, int C, typename B> 00191 inline void compute(const Matrix<R,C,Precision,B>& m){ 00192 SizeMismatch<R, C>::test(m.num_rows(), m.num_cols()); 00193 SizeMismatch<R, Size>::test(m.num_rows(), my_evectors.num_rows()); 00194 Internal::ComputeSymEigen<Size>::compute(m, my_evectors, my_evalues); 00195 } 00196 00197 /// Calculate result of multiplying the (pseudo-)inverse of M by a vector. 00198 /// For a vector \f$b\f$, this calculates \f$M^{\dagger}b\f$ by back substitution 00199 /// (i.e. without explictly calculating the (pseudo-)inverse). 00200 /// See the SVD detailed description for a description of condition variables. 00201 template <int S, typename P, typename B> 00202 Vector<Size, Precision> backsub(const Vector<S,P,B>& rhs) const { 00203 return (my_evectors.T() * diagmult(get_inv_diag(Internal::symeigen_condition_no),(my_evectors * rhs))); 00204 } 00205 00206 /// Calculate result of multiplying the (pseudo-)inverse of M by another matrix. 00207 /// For a matrix \f$A\f$, this calculates \f$M^{\dagger}A\f$ by back substitution 00208 /// (i.e. without explictly calculating the (pseudo-)inverse). 00209 /// See the SVD detailed description for a description of condition variables. 00210 template <int R, int C, typename P, typename B> 00211 Matrix<Size,C, Precision> backsub(const Matrix<R,C,P,B>& rhs) const { 00212 return (my_evectors.T() * diagmult(get_inv_diag(Internal::symeigen_condition_no),(my_evectors * rhs))); 00213 } 00214 00215 /// Calculate (pseudo-)inverse of the matrix. This is not usually needed: 00216 /// if you need the inverse just to multiply it by a matrix or a vector, use 00217 /// one of the backsub() functions, which will be faster. 00218 /// See the SVD detailed description for a description of the pseudo-inverse 00219 /// and condition variables. 00220 Matrix<Size, Size, Precision> get_pinv(const double condition=Internal::symeigen_condition_no) const { 00221 return my_evectors.T() * diagmult(get_inv_diag(condition),my_evectors); 00222 } 00223 00224 /// Calculates the reciprocals of the eigenvalues of the matrix. 00225 /// The vector <code>invdiag</code> lists the eigenvalues in order, from 00226 /// the largest (i.e. smallest reciprocal) to the smallest. 00227 /// These are also the diagonal values of the matrix \f$Lambda^{-1}\f$. 00228 /// Any eigenvalues which are too small are set to zero (see the SVD 00229 /// detailed description for a description of the and condition variables). 00230 Vector<Size, Precision> get_inv_diag(const double condition) const { 00231 Precision max_diag = -my_evalues[0] > my_evalues[my_evalues.size()-1] ? -my_evalues[0]:my_evalues[my_evalues.size()-1]; 00232 Vector<Size, Precision> invdiag(my_evalues.size()); 00233 for(int i=0; i<my_evalues.size(); i++){ 00234 if(fabs(my_evalues[i]) * condition > max_diag) { 00235 invdiag[i] = 1/my_evalues[i]; 00236 } else { 00237 invdiag[i]=0; 00238 } 00239 } 00240 return invdiag; 00241 } 00242 00243 /// Returns the eigenvectors of the matrix. 00244 /// This returns \f$U^T\f$, so that the rows of the matrix are the eigenvectors, 00245 /// which can be extracted using usual Matrix::operator[]() subscript operator. 00246 /// They are returned in order of the size of the corresponding eigenvalue, i.e. 00247 /// the vector with the largest eigenvalue is first. 00248 Matrix<Size,Size,Precision>& get_evectors() {return my_evectors;} 00249 00250 /**\overload 00251 */ 00252 const Matrix<Size,Size,Precision>& get_evectors() const {return my_evectors;} 00253 00254 00255 /// Returns the eigenvalues of the matrix. 00256 /// The eigenvalues are listed in order, from the smallest to the largest. 00257 /// These are also the diagonal values of the matrix \f$\Lambda\f$. 00258 Vector<Size, Precision>& get_evalues() {return my_evalues;} 00259 /**\overload 00260 */ 00261 const Vector<Size, Precision>& get_evalues() const {return my_evalues;} 00262 00263 /// Is the matrix positive definite? 00264 bool is_posdef() const { 00265 for (int i = 0; i < my_evalues.size(); ++i) { 00266 if (my_evalues[i] <= 0.0) 00267 return false; 00268 } 00269 return true; 00270 } 00271 00272 /// Is the matrix negative definite? 00273 bool is_negdef() const { 00274 for (int i = 0; i < my_evalues.size(); ++i) { 00275 if (my_evalues[i] >= 0.0) 00276 return false; 00277 } 00278 return true; 00279 } 00280 00281 /// Get the determinant of the matrix 00282 Precision get_determinant () const { 00283 Precision det = 1.0; 00284 for (int i = 0; i < my_evalues.size(); ++i) { 00285 det *= my_evalues[i]; 00286 } 00287 return det; 00288 } 00289 00290 /// Calculate the square root of a matrix which is a matrix M 00291 /// such that M.T*M=A. 00292 Matrix<Size, Size, Precision> get_sqrtm () const { 00293 Vector<Size, Precision> diag_sqrt(my_evalues.size()); 00294 // In the future, maybe throw an exception if an 00295 // eigenvalue is negative? 00296 for (int i = 0; i < my_evalues.size(); ++i) { 00297 diag_sqrt[i] = std::sqrt(my_evalues[i]); 00298 } 00299 return my_evectors.T() * diagmult(diag_sqrt, my_evectors); 00300 } 00301 00302 /// Calculate the inverse square root of a matrix which is a 00303 /// matrix M such that M.T*M=A^-1. 00304 /// 00305 /// Any square-rooted eigenvalues which are too small are set 00306 /// to zero (see the SVD detailed description for a 00307 /// description of the condition variables). 00308 Matrix<Size, Size, Precision> get_isqrtm (const double condition=Internal::symeigen_condition_no) const { 00309 Vector<Size, Precision> diag_isqrt(my_evalues.size()); 00310 00311 // Because sqrt is a monotonic-preserving transformation, 00312 Precision max_diag = -my_evalues[0] > my_evalues[my_evalues.size()-1] ? (-std::sqrt(my_evalues[0])):(std::sqrt(my_evalues[my_evalues.size()-1])); 00313 // In the future, maybe throw an exception if an 00314 // eigenvalue is negative? 00315 for (int i = 0; i < my_evalues.size(); ++i) { 00316 diag_isqrt[i] = std::sqrt(my_evalues[i]); 00317 if(fabs(diag_isqrt[i]) * condition > max_diag) { 00318 diag_isqrt[i] = 1/diag_isqrt[i]; 00319 } else { 00320 diag_isqrt[i] = 0; 00321 } 00322 } 00323 return my_evectors.T() * diagmult(diag_isqrt, my_evectors); 00324 } 00325 00326 private: 00327 // eigen vectors laid out row-wise so evectors[i] is the ith evector 00328 Matrix<Size,Size,Precision> my_evectors; 00329 00330 Vector<Size, Precision> my_evalues; 00331 }; 00332 00333 } 00334 00335 #endif