Main Page   Groups   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Concepts

vnl_svd.h

Go to the documentation of this file.
00001 #ifndef vnl_svd_h_
00002 #define vnl_svd_h_
00003 
00004 //:
00005 //  \file
00006 //  \brief Holds the singular value decomposition of a vnl_matrix.
00007 //  \author Andrew W. Fitzgibbon, Oxford IERG, 15 Jul 96
00008 //
00009 //  Modifications
00010 // F. Schaffalitzky, Oxford IESRG, 26 Mar 1999
00011 //     1. The singular values are now stored as reals (not complexes) when T is
00012 //        complex.
00013 //     2. Fixed bug : for complex T, matrices have to be conjugated as well as 
00014 //        transposed.
00015 
00016 #include <vnl/vnl_numeric_traits.h>
00017 #include <vnl/vnl_vector.h>
00018 #include <vnl/vnl_matrix.h>
00019 #include <vnl/vnl_diag_matrix.h>
00020 #include <vcl_iosfwd.h>
00021 
00022 //: Holds the singular value decomposition of a vnl_matrix.
00023 //
00024 //  @{ The class holds three matrices U, W, V such that the original matrix 
00025 //  $M = U W V^\top$.The DiagMatrix W stores the singular values in decreasing
00026 //  order.  The columns of U which correspond to the nonzero singular values
00027 //  form a basis for range of M, while the columns of V corresponding to the
00028 //  zero singular values are the nullspace. @}
00029 //    
00030 //  The SVD is computed at construction time, and enquiries may then be made
00031 //  of the SVD.  In particular, this allows easy access to multiple
00032 //  right-hand-side solves without the bother of putting all the RHS's into a
00033 //  Matrix.
00034 //
00035 //  This class is supplied even though there is an existing vnl_matrix method
00036 //  for several reasons:
00037 //
00038 //  It is more convenient to use as it manages all the storage for
00039 //  the U,S,V matrices, allowing repeated queries of the same SVD
00040 //  results.
00041 //
00042 //  It avoids namespace clutter in the Matrix class.   While svd()
00043 //  is a perfectly reasonable method for a Matrix, there are many other
00044 //  decompositions that might be of interest, and adding them all would
00045 //  make for a very large Matrix class.
00046 //
00047 //  It demonstrates the holder model of compute class, implementing an
00048 //  algorithm on an object without adding a member that may not be of
00049 //  general interest.  A similar pattern can be used for other
00050 //  decompositions which are not defined as members of the library Matrix
00051 //  class.
00052 //
00053 //  It extends readily to n-ary operations, such as generalized
00054 //  eigensystems, which cannot be members of just one matrix.
00055 
00056 
00057 template <class T> class vnl_svd;
00058 // templated friends must be declared to be templated :
00059 template <class T> vcl_ostream& operator<<(vcl_ostream&, 
00060                                            vnl_svd<T> const& svd);
00061 
00062 export template <class T>
00063 class vnl_svd {
00064 public:
00065   // The singular values of a matrix of complex<T> are of type T,not complex<T>
00066   typedef typename vnl_numeric_traits<T>::abs_t singval_t;
00067 
00068 //: @{
00069 // Construct an vnl_svd<T> object from $m \times n$ matrix $M$.  The 
00070 // vnl_svd<T> object contains matrices $U, W, V$ such that $U W V^\top = M$.
00071 // \par
00072 // Uses linpack routine DSVDC to calculate an ``economy-size'' SVD
00073 // where the returned $U$ is the same size as $M$, while $W$ and $V$
00074 // are both $n \times n$.  This is efficient for
00075 // large rectangular solves where $m > n$, typical in least squares.
00076 // @}
00077 // The optional argument zero_out_tol is used to mark the zero singular
00078 // values: If nonnegative, any s.v. smaller than zero_out_tol in
00079 // absolute value is set to zero.  If zero_out_tol is negative, the
00080 // zeroing is relative to |zero_out_tol| * sigma_max();
00081 //
00082   vnl_svd(vnl_matrix<T> const &M, double zero_out_tol = 0.0);
00083  ~vnl_svd() {}
00084 
00085   // Data Access---------------------------------------------------------------
00086   void            zero_out_absolute(double tol = 1e-8); //sqrt(machine epsilon)
00087   void            zero_out_relative(double tol = 1e-8); //sqrt(machine epsilon)
00088   int             singularities () const { return W_.n() - rank(); }
00089   int             rank () const { return rank_; }
00090   singval_t       well_condition () const { return sigma_min()/sigma_max(); }
00091   singval_t       determinant_magnitude () const;
00092   singval_t       norm() const;
00093 
00094 //: Return the matrix U and the (i,j)th entry (to avoid svd.U()(i,j); ).
00095   vnl_matrix<T>      & U()       { return U_; }
00096   vnl_matrix<T> const& U() const { return U_; }
00097   T U(int i, int j) { return U_(i,j); }
00098 
00099 //: Get at DiagMatrix (q.v.) of singular values,sorted from largest to smallest
00100   vnl_diag_matrix<singval_t>       & W()             { return W_; }
00101   vnl_diag_matrix<singval_t> const & W() const       { return W_; }
00102   vnl_diag_matrix<singval_t>       & Winverse()      { return Winverse_; }
00103   vnl_diag_matrix<singval_t> const & Winverse() const { return Winverse_; }
00104   singval_t                   & W(int i, int j) { return W_(i,j); }
00105   singval_t                   & W(int i)        { return W_(i,i); }
00106   singval_t     sigma_max() const { return W_(0,0); }       // largest
00107   singval_t     sigma_min() const { return W_(n_-1,n_-1); } // smallest
00108 
00109 //: Return the matrix V.
00110   vnl_matrix<T>      & V()       { return V_; }
00111   vnl_matrix<T> const& V() const { return V_; }
00112   T V(int i, int j) { return V_(i,j); }
00113   
00114   //
00115   vnl_matrix<T> inverse () const;
00116   vnl_matrix<T> pinverse () const; // pseudo-inverse (for non-square matrix).
00117   vnl_matrix<T> tinverse () const;
00118   vnl_matrix<T> recompose () const;
00119 
00120   //
00121   vnl_matrix<T> solve (vnl_matrix<T> const& rhs) const;
00122   vnl_vector<T> solve (vnl_vector<T> const& rhs) const;
00123   void          solve (T const *rhs, T *lhs) const; // min ||A*lhs - rhs||
00124   void solve_preinverted(vnl_vector<T> const& rhs, vnl_vector<T>* out) const;
00125 
00126   //
00127   vnl_matrix<T> nullspace() const;
00128   vnl_matrix<T> left_nullspace() const;
00129 
00130   vnl_matrix<T> nullspace(int required_nullspace_dimension) const;
00131   vnl_matrix<T> left_nullspace(int required_nullspace_dimension) const;
00132 
00133   vnl_vector<T> nullvector() const;
00134   vnl_vector<T> left_nullvector() const;
00135 
00136 //  friend ostream& operator<<(ostream&, const vnl_svd<T>& svd);
00137 
00138 private:
00139 
00140   // Data Members--------------------------------------------------------------
00141 
00142   int m_, n_;              // Size of M, local cache.
00143   vnl_matrix<T> U_;        // Columns Ui are basis for range of M for Wi != 0
00144   vnl_diag_matrix<singval_t> W_;// Singular values, sorted in decreasing order
00145   vnl_diag_matrix<singval_t> Winverse_;
00146   vnl_matrix<T> V_;       // Columns Vi are basis for nullspace of M for Wi = 0
00147   unsigned rank_;
00148   bool have_max_; 
00149   singval_t max_;
00150   bool have_min_;
00151   singval_t min_;
00152   double last_tol_;
00153 
00154 protected:
00155 
00156   // Constructors/Destructors--------------------------------------------------
00157 
00158   // Disallow assignment until awf decides whether its a good idea
00159   vnl_svd(vnl_svd<T> const & that);
00160   vnl_svd<T>& operator=(vnl_svd<T> const & that);
00161 };
00162 
00163 template <class T>
00164 inline
00165 vnl_matrix<T> vnl_svd_inverse(vnl_matrix<T> const& m)
00166 {
00167   return vnl_svd<T>(m).inverse();
00168 }
00169 
00170 #endif // vnl_svd_h_

Generated at Wed Mar 12 01:13:16 2003 for ITK by doxygen 1.2.15 written by Dimitri van Heesch, © 1997-2000