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 }