Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

svd.h

Go to the documentation of this file.
00001 // This may look like C code, but it is really -*- C++ -*-
00002 /*
00003  ************************************************************************
00004  *
00005  *                        Numerical Math Package
00006  *
00007  *      Singular Value Decomposition of a rectangular matrix
00008  *                         A = U * Sig * V'
00009  *                      and its applications
00010  *
00011  * In the decomposition above, matrices U and V are orthogonal; matrix
00012  * Sig is a diagonal matrix: its diagonal elements, which are all
00013  * non-negative, are singular values (numbers) of the original matrix A.
00014  * In another interpretation, the singular values are eigenvalues
00015  * of matrix A'A.
00016  *
00017  * $Id: svd.h,v 1.1 2004/05/21 21:02:52 maxx Exp $
00018  *
00019  ************************************************************************
00020  */
00021 
00022 #ifndef __GNUC__
00023 #pragma once
00024 #else
00025 #pragma interface
00026 #endif
00027 #ifndef _svd_h
00028 #define _svd_h 1
00029 
00030 
00031 #include "LinAlg.h"
00032 
00033                                 // A class that holds U,V,Sig - the singular
00034                                 // value decomposition of a matrix
00035 namespace linalg 
00036 {
00037     using namespace linalg;
00038     
00039 class SVD
00040 {
00041   const int M,N;                        // Dimensions of the problem (M>=N)
00042   Matrix U;                             // M*M orthogonal matrix U
00043   Matrix V;                             // N*N orthogonal matrix V
00044   Vector sig;                           // Vector(1:N) of N onordered singular
00045                                         // values
00046   
00047                                         // Internal procedures used in SVD
00048  inline double left_householder(Matrix& A, const int i);
00049  inline double right_householder(Matrix& A, const int i);
00050  double bidiagonalize(Vector& super_diag, const Matrix& __A);
00051 
00052  inline void rotate(Matrix& U, const int i, const int j,
00053                     const double cos_ph, const double sin_ph);
00054  inline void rip_through(
00055         Vector& super_diag, const int k, const int l, const double eps);
00056  inline int get_submatrix_to_work_on(
00057         Vector& super_diag, const int k, const double eps);
00058  void diagonalize(Vector& super_diag, const double eps);
00059  
00060 public:
00061   SVD(const Matrix& A);                 // Decompose Matrix A, of M rows and
00062                                         // N columns, M>=N
00063   
00064                                         // Regularization: make all sig(i)
00065                                         // that are smaller than min_sig
00066                                         // exactly zeros
00067   void cut_off(const double min_sig);
00068   
00069                                 // Inquiries
00070   const Matrix& q_U(void) const         { return U; }
00071   const Matrix& q_V(void) const         { return V; }
00072   const Vector& q_sig(void) const       { return sig; }
00073 
00074   operator MinMax(void) const;          // Return min and max singular values
00075   double q_cond_number(void) const;     // sig_max/sig_min
00076   void info(void) const;                // Print the info about the SVD
00077 };
00078 
00079 /*
00080  *------------------------------------------------------------------------
00081  * An application of SVD: a regularized solution of a set of simultaneous
00082  * linear equations Ax=B
00083  * B can be either a vector (Mx1-matrix), or a full-blown matrix.
00084  * Note, if B=Unit(A), SVD_inv_mult class gives a (pseudo)inverse matrix;
00085  * btw, A doesn't even have to be a square matrix.
00086  * 
00087  * In the case of a rectangular MxN matrix A with M>N, the set
00088  * Ax=b is obviously overspecified. The solution x produced by a
00089  * SVD_inv_mult below is the least-norm solution, which is the solution
00090  * returned by a least-squares method.
00091  * 
00092  * tau is a regularization criterion: only singular values bigger than tau
00093  * would participate. If tau is not given, it's assumed to be
00094  * dim(sig)*max(sig)*FLT_EPSILON
00095  *
00096  * Use the SVD_inv_mult object as follows:
00097  *      SVD svd(A);
00098  *      cout << "condition number of matrix A " << svd.q_cond_number();
00099  *      Vector x = SVD_inv_mult(svd,b);         // Solution of Ax=b
00100  *
00101  */
00102 
00103 class SVD_inv_mult : public LazyMatrix
00104 {
00105   const SVD& svd;
00106   const Matrix& B;
00107   double tau;                   // smallness threshold for sig[i]
00108   bool are_zero_coeff;          // true if A had an incomplete rank
00109   void fill_in(Matrix& m) const;
00110 public: 
00111   SVD_inv_mult(const SVD& _svd, const Matrix& __B,const double tau=0);
00112 };
00113 
00114 #endif
00115 }

Generated on Wed Dec 15 21:20:32 2004 for vuVolume by  doxygen 1.3.9.1