core/vnl/algo/vnl_symmetric_eigensystem.h
Go to the documentation of this file.
00001 // This is core/vnl/algo/vnl_symmetric_eigensystem.h
00002 #ifndef vnl_symmetric_eigensystem_h_
00003 #define vnl_symmetric_eigensystem_h_
00004 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00005 #pragma interface
00006 #endif
00007 //:
00008 // \file
00009 // \brief Find eigenvalues of a symmetric matrix
00010 //
00011 //    vnl_symmetric_eigensystem_compute()
00012 //    solves the eigenproblem $A x = \lambda x$, with $A$ symmetric.
00013 //    The resulting eigenvectors and values are sorted in increasing order
00014 //    so <CODE> V.column(0) </CODE> is the eigenvector corresponding to the
00015 //    smallest eigenvalue.
00016 //
00017 //    As a matrix decomposition, this is $A = V D V^t$
00018 //
00019 //    Uses the EISPACK routine RS, which in turn calls TRED2 to reduce A
00020 //    to tridiagonal form, followed by TQL2, to find the eigensystem.
00021 //    This is summarized in Golub and van Loan, pgf 8.2.  The following are
00022 //    the original subroutine headers:
00023 //
00024 // \remark TRED2 is a translation of the Algol procedure tred2,
00025 //     Num. Math. 11, 181-195(1968) by Martin, Reinsch, and Wilkinson.
00026 //     Handbook for Auto. Comp., Vol.ii-Linear Algebra, 212-226(1971).
00027 //
00028 // \remark This subroutine reduces a real symmetric matrix to a
00029 //     symmetric tridiagonal matrix using and accumulating
00030 //     orthogonal similarity transformations.
00031 //
00032 // \remark TQL2 is a translation of the Algol procedure tql2,
00033 //     Num. Math. 11, 293-306(1968) by Bowdler, Martin, Reinsch, and Wilkinson.
00034 //     Handbook for Auto. Comp., Vol.ii-Linear Algebra, 227-240(1971).
00035 //
00036 // \remark This subroutine finds the eigenvalues and eigenvectors
00037 //     of a symmetric tridiagonal matrix by the QL method.
00038 //     the eigenvectors of a full symmetric matrix can also
00039 //     be found if  tred2  has been used to reduce this
00040 //     full matrix to tridiagonal form.
00041 //
00042 // \author Andrew W. Fitzgibbon, Oxford RRG
00043 // \date   29 Aug 96
00044 //
00045 // \verbatim
00046 //  Modifications
00047 //   fsm, 5 March 2000: templated
00048 //   dac (Manchester) 28/03/2001: tidied up documentation
00049 //   Feb.2002 - Peter Vanroose - brief doxygen comment placed on single line
00050 //   Jan.2003 - Peter Vanroose - added missing implementation for solve(b,x)
00051 //   Mar.2010 - Peter Vanroose - also made vnl_symmetric_eigensystem_compute()
00052 //                               & vnl_symmetric_eigensystem_compute_eigenvals() templated
00053 // \endverbatim
00054 
00055 #include <vnl/vnl_matrix.h>
00056 #include <vnl/vnl_diag_matrix.h>
00057 
00058 //: Find eigenvalues of a symmetric 3x3 matrix
00059 // Eigenvalues will be returned so that l1 <= l2 <= l3.
00060 // \verbatim
00061 // Matrix is   M11  M12  M13
00062 //             M12  M22  M23
00063 //             M13  M23  M33
00064 // \endverbatim
00065 template <class T>
00066 void vnl_symmetric_eigensystem_compute_eigenvals(
00067        T M11, T M12, T M13,
00068               T M22, T M23,
00069                      T M33,
00070        T &l1, T &l2, T &l3);
00071 
00072 //: Find eigenvalues of a symmetric matrix
00073 template <class T>
00074 bool vnl_symmetric_eigensystem_compute(vnl_matrix<T> const & A,
00075                                        vnl_matrix<T> & V,
00076                                        vnl_vector<T> & D);
00077 
00078 //: Computes and stores the eigensystem decomposition of a symmetric matrix.
00079 
00080 VCL_TEMPLATE_EXPORT template <class T>
00081 class vnl_symmetric_eigensystem
00082 {
00083  public:
00084   //: Solve real symmetric eigensystem $A x = \lambda x$
00085   vnl_symmetric_eigensystem(vnl_matrix<T> const & M);
00086 
00087  protected:
00088   // need this here to get inits in correct order, but still keep gentex
00089   // in the right order.
00090   int n_;
00091 
00092  public:
00093   //: Public eigenvectors.
00094   //  After construction, the columns of V are the eigenvectors, sorted by
00095   // increasing eigenvalue, from most negative to most positive.
00096   vnl_matrix<T> V;
00097 
00098   //: Public eigenvalues.
00099   //  After construction,  D contains the eigenvalues, sorted as described above.
00100   //  Note that D is a vnl_diag_matrix, and is therefore stored as a vcl_vector while behaving as a matrix.
00101   vnl_diag_matrix<T> D;
00102 
00103   //: Recover specified eigenvector after computation.
00104   vnl_vector<T> get_eigenvector(int i) const;
00105 
00106   //: Recover specified eigenvalue after computation.
00107   T             get_eigenvalue(int i) const;
00108 
00109   //: Convenience method to get least-squares nullvector.
00110   // It is deliberate that the signature is the same as on vnl_svd<T>.
00111   vnl_vector<T> nullvector() const { return get_eigenvector(0); }
00112 
00113   //: Return the matrix $V  D  V^\top$.
00114   //  This can be useful if you've modified $D$.  So an inverse is obtained using
00115   // \code
00116   //   vnl_symmetric_eigensystem} eig(A);
00117   //   eig.D.invert_in_place}();
00118   //   vnl_matrix<double> Ainverse = eig.recompose();
00119   // \endcode
00120 
00121   vnl_matrix<T> recompose() const { return V * D * V.transpose(); }
00122 
00123   //: return product of eigenvalues.
00124   T determinant() const;
00125 
00126   //: return the pseudoinverse.
00127   vnl_matrix<T> pinverse() const;
00128 
00129   //: return the square root, if positive semi-definite.
00130   vnl_matrix<T> square_root() const;
00131 
00132   //: return the inverse of the square root, if positive semi-definite.
00133   vnl_matrix<T> inverse_square_root() const;
00134 
00135   //: Solve LS problem M x = b
00136   vnl_vector<T> solve(vnl_vector<T> const & b);
00137 
00138   //: Solve LS problem M x = b
00139   void solve(vnl_vector<T> const & b, vnl_vector<T> * x) { *x = solve(b); }
00140 };
00141 
00142 #define VNL_SYMMETRIC_EIGENSYSTEM_INSTANTIATE(T) \
00143 extern "please include vnl/algo/vnl_symmetric_eigensystem.txx first"
00144 
00145 #endif // vnl_symmetric_eigensystem_h_