home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / LinAlg 3.1 / LinAlg / svd.h < prev    next >
Encoding:
Text File  |  1995-12-22  |  3.4 KB  |  108 lines  |  [TEXT/ALFA]
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *              Numerical Math Package
  6.  *
  7.  *    Singular Value Decomposition of a rectangular matrix
  8.  *               A = U * Sig * V'
  9.  *            and its applications
  10.  *
  11.  * In the decomposition above, matrices U and V are orthogonal; matrix
  12.  * Sig is a diagonal matrix: its diagonal elements, which are all
  13.  * non-negative, are singular values (numbers) of the original matrix A.
  14.  * In another interpretation, the singular values are eigenvalues
  15.  * of matrix A'A.
  16.  *
  17.  * $Id$
  18.  *
  19.  ************************************************************************
  20.  */
  21.  
  22. #ifndef __GNUC__
  23. #pragma once
  24. #endif
  25. #ifndef _svd_h
  26. #define _svd_h 1
  27.  
  28. #pragma interface
  29.  
  30. #include "LinAlg.h"
  31.  
  32.                 // A class that holds U,V,Sig - the singular
  33.                 // value decomposition of a matrix
  34. class SVD
  35. {
  36.   const int M,N;            // Dimensions of the problem (M>=N)
  37.   Matrix U;                // M*M orthogonal matrix U
  38.   Matrix V;                // N*N orthogonal matrix V
  39.   Vector sig;                // Vector(1:N) of N onordered singular
  40.                       // values
  41.   
  42.                       // Internal procedures used in SVD
  43.  inline double left_hausholder(Matrix& A, const int i);
  44.  inline double right_hausholder(Matrix& A, const int i);
  45.  double bidiagonalize(Vector& super_diag, const Matrix& _A);
  46.  
  47.  inline void rotate(Matrix& U, const int i, const int j,
  48.             const double cos_ph, const double sin_ph);
  49.  inline void rip_through(
  50.     Vector& super_diag, const int k, const int l, const double eps);
  51.  inline int get_submatrix_to_work_on(
  52.     Vector& super_diag, const int k, const double eps);
  53.  void diagonalize(Vector& super_diag, const double eps);
  54.  
  55. public:
  56.   SVD(const Matrix& A);            // Decompose Matrix A, of M rows and
  57.                       // N columns, M>=N
  58.   
  59.                       // Regularization: make all sig(i)
  60.                       // that are smaller than min_sig
  61.                       // exactly zeros
  62.   void cut_off(const double min_sig);
  63.   
  64.                   // Inquiries
  65.   const Matrix& q_U(void) const        { return U; }
  66.   const Matrix& q_V(void) const        { return V; }
  67.   const Vector& q_sig(void) const    { return sig; }
  68.  
  69.   operator MinMax(void) const;        // Return min and max singular values
  70.   double q_cond_number(void) const;    // sig_max/sig_min
  71.   void info(void) const;        // Print the info about the SVD
  72. };
  73.  
  74. /*
  75.  *------------------------------------------------------------------------
  76.  * Application of SVD to a regularized solution of a set of simultaneous
  77.  * linear equations Ax=B
  78.  * B can be either a vector (Mx1-matrix), or a full-blown matrix.
  79.  * Note, if B=Unit(A), SVD_inv_mult class gives a (pseudo)inverse matrix;
  80.  * btw, A doesn't even have to be a square matrix.
  81.  * In the case of a rectangular matrix MxN matrix A with M>N, the set
  82.  * Ax=b is obviously overspecified. The solution x produced by SVD_inv_mult
  83.  * is the least-norm solution: the solution returned by the least-squares
  84.  * method.
  85.  * tau is a regularization criterion: only singular values bigger than tau
  86.  * would participate. If tau is not given, it's assumed to be
  87.  * dim(sig)*max(sig)*FLT_EPSILON
  88.  *
  89.  * Use the SVD_inv_mult object as follows:
  90.  *     SVD svd(A);
  91.  *    cout << "condition number of matrix A " << svd.q_cond_number();
  92.  *    Vector x = SVD_inv_mult(svd,b);        // Solution of Ax=b
  93.  *
  94.  */
  95.  
  96. class SVD_inv_mult : public LazyMatrix
  97. {
  98.   const SVD& svd;
  99.   const Matrix& B;
  100.   double tau;            // smallness threshold for sig[i]
  101.   bool are_zero_coeff;        // true if A had an incomplete rank
  102.   void fill_in(Matrix& m) const;
  103. public:    
  104.   SVD_inv_mult(const SVD& _svd, const Matrix& _B,const double tau=0);
  105. };
  106.  
  107. #endif
  108.