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

svd.cpp

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  *      Singular Value Decomposition of a rectangular matrix
00007  *                           A = U * Sig * V'
00008  *
00009  * where matrices U and V are orthogonal and Sig is a digonal matrix.
00010  *
00011  * The singular value decomposition is performed by constructing an SVD
00012  * object from an M*N matrix A with M>=N (that is, at least as many rows
00013  * as columns). Note, in case M > N, matrix Sig has to be a M*N diagonal
00014  * matrix. However, it has only N diag elements, which we store in a 1:N
00015  * Vector sig.
00016  *
00017  * Algorithm
00018  *      Bidiagonalization with Householder reflections followed by a
00019  * modification of a QR-algorithm. For more details, see
00020  *   G.E. Forsythe, M.A. Malcolm, C.B. Moler
00021  *   Computer methods for mathematical computations. - Prentice-Hall, 1977
00022  * However, in the persent implementation, matrices U and V are computed
00023  * right away rather than delayed until after all Householder reflections.
00024  *
00025  * This code is based for the most part on a Algol68 code I wrote
00026  * ca. 1987
00027  *
00028  * $Id: svd.cpp,v 1.2 2004/10/02 22:23:32 maxx Exp $
00029  *
00030  ************************************************************************
00031  */
00032 
00033 #if defined(__GNUC__)
00034 #pragma implementation
00035 #endif
00036 
00037 #include <math.h>
00038 #include "svd.h"
00039 #include <float.h>
00040 #include "LAStreams.h"
00041 
00042 namespace linalg 
00043 {
00044     using namespace linalg;
00045     
00046 /*
00047  *------------------------------------------------------------------------
00048  *                              Bidiagonalization
00049  */
00050 
00051  /*
00052  *                      Left Householder Transformations
00053  *
00054  * Zero out an entire subdiagonal of the i-th column of A and compute the
00055  * modified A[i,i] by multiplication (UP' * A) with a matrix UP'
00056  *   (1)  UP' = E - UPi * UPi' / beta
00057  *
00058  * where a column-vector UPi is as follows
00059  *   (2)  UPi = [ (i-1) zeros, A[i,i] + Norm, vector A[i+1:M,i] ]
00060  * where beta = UPi' * A[,i] and Norm is the norm of a vector A[i:M,i]
00061  * (sub-diag part of the i-th column of A). Note we assign the Norm the
00062  * same sign as that of A[i,i].
00063  * By construction, (1) does not affect the first (i-1) rows of A. Since
00064  * A[*,1:i-1] is bidiagonal (the result of the i-1 previous steps of
00065  * the bidiag algorithm), transform (1) doesn't affect these i-1 columns
00066  * either as one can easily verify.
00067  * The i-th column of A is transformed as
00068  *   (3)  UP' * A[*,i] = A[*,i] - UPi
00069  * (since UPi'*A[*,i]/beta = 1 by construction of UPi and beta)
00070  * This means effectively zeroing out A[i+1:M,i] (the entire subdiagonal
00071  * of the i-th column of A) and replacing A[i,i] with the -Norm. Thus
00072  * modified A[i,i] is returned by the present function.
00073  * The other (i+1:N) columns of A are transformed as
00074  *    (4)  UP' * A[,j] = A[,j] - UPi * ( UPi' * A[,j] / beta )
00075  * Note, due to (2), only elements of rows i+1:M actually  participate
00076  * in above transforms; the upper i-1 rows of A are not affected.
00077  * As was mentioned earlier,
00078  * (5)  beta = UPi' * A[,i] = (A[i,i] + Norm)*A[i,i] + A[i+1:M,i]*A[i+1:M,i]
00079  *      = ||A[i:M,i]||^2 + Norm*A[i,i] = Norm^2 + Norm*A[i,i]
00080  * (note the sign of the Norm is the same as A[i,i])
00081  * For extra precision, vector UPi (and so is Norm and beta) are scaled,
00082  * which would not affect (4) as easy to see.
00083  *
00084  * To satisfy the definition
00085  *   (6)  .SIG = U' A V
00086  * the result of consecutive transformations (1) over matrix A is accumulated
00087  * in matrix U' (which is initialized to be a unit matrix). At each step,
00088  * U' is left-multiplied by UP' = UP (UP is symmetric by construction,
00089  * see (1)). That is, U is right-multiplied by UP, that is, rows of U are
00090  * transformed similarly to columns of A, see eq. (4). We also keep in mind
00091  * that multiplication by UP at the i-th step does not affect the first i-1
00092  * columns of U.
00093  * Note that the vector UPi doesn't have to be allocated explicitly: its
00094  * first i-1 components are zeros (which we can always imply in computations),
00095  * and the rest of the components (but the UPi[i]) are the same as those
00096  * of A[i:M,i], the subdiagonal of A[,i]. This column, A[,i] is affected only
00097  * trivially as explained above, that is, we don't need to carry this
00098  * transformation explicitly (only A[i,i] is going to be non-trivially
00099  * affected, that is, replaced by -Norm, but we will use sig[i] to store
00100  * the result).
00101  *
00102  */
00103  
00104  inline double SVD::left_householder(Matrix& A, const int i)
00105  {                                      // Note that only UPi[i:M] matter
00106    IRange range = IRange::from(i - A.q_col_lwb());
00107    MatrixColumn x(A,i);
00108    LAStreamOut UPi_str(x,range);
00109    register int j;
00110    REAL scale = 0;                      // Compute the scaling factor
00111    while( !UPi_str.eof() )
00112      scale += fabs(UPi_str.get());
00113    if( scale == 0 )                     // If A[,i] is a null vector, no
00114      return 0;                          // transform is required
00115 
00116    UPi_str.rewind();
00117    double Norm_sqr = 0;                 // Scale UPi (that is, A[,i])
00118    while( !UPi_str.eof() )              // and compute its norm, Norm^2
00119      Norm_sqr += sqr(UPi_str.get() /= scale);
00120    UPi_str.rewind();
00121    double new_Aii = sqrt(Norm_sqr);     // new_Aii = -Norm, Norm has the
00122    if( UPi_str.peek() > 0 )             // same sign as Aii (that is, UPi[i])
00123      new_Aii = -new_Aii;
00124    const float beta = - UPi_str.peek()*new_Aii + Norm_sqr;
00125    UPi_str.peek() -= new_Aii;           // UPi[i] = A[i,i] - (-Norm)
00126    
00127    for(j=i+1; j<=N; j++)        // Transform i+1:N columns of A
00128    {
00129      MatrixColumn xx(A,j);
00130      LAStreamOut Aj_str(xx,range);
00131      REAL factor = 0;
00132      while( !UPi_str.eof() )
00133        factor += UPi_str.get() * Aj_str.get();  // Compute UPi' * A[,j]
00134      factor /= beta;
00135      for(UPi_str.rewind(), Aj_str.rewind(); !UPi_str.eof(); )
00136        Aj_str.get() -= UPi_str.get() * factor;
00137      UPi_str.rewind();
00138    }
00139 
00140    for(j=1; j<=M; j++)  // Accumulate the transform in U
00141    {
00142      MatrixRow y(U,j);
00143      LAStrideStreamOut Uj_str(y,range);
00144      REAL factor = 0;
00145      while( !UPi_str.eof() )
00146        factor += UPi_str.get() * Uj_str.get();  // Compute  U[j,] * UPi
00147      factor /= beta;
00148      for(UPi_str.rewind(),Uj_str.rewind(); !UPi_str.eof(); )
00149         Uj_str.get() -= UPi_str.get() * factor;
00150      UPi_str.rewind();
00151    }
00152    return new_Aii * scale;              // Scale new Aii back (our new Sig[i])
00153  }
00154  
00155 /*
00156  *                      Right Householder Transformations
00157  *
00158  * Zero out i+2:N elements of a row A[i,] of matrix A by right
00159  * multiplication (A * VP) with a matrix VP
00160  *   (1)  VP = E - VPi * VPi' / beta
00161  *
00162  * where a vector-column .VPi is as follows
00163  *   (2)  VPi = [ i zeros, A[i,i+1] + Norm, vector A[i,i+2:N] ]
00164  * where beta = A[i,] * VPi and Norm is the norm of a vector A[i,i+1:N]
00165  * (right-diag part of the i-th row of A). Note we assign the Norm the
00166  * same sign as that of A[i,i+1].
00167  * By construction, (1) does not affect the first i columns of A. Since
00168  * A[1:i-1,] is bidiagonal (the result of the previous steps of
00169  * the bidiag algorithm), transform (1) doesn't affect these i-1 rows
00170  * either as one can easily verify.
00171  * The i-th row of A is transformed as
00172  *  (3)  A[i,*] * VP = A[i,*] - VPi'
00173  * (since A[i,*]*VPi/beta = 1 by construction of VPi and beta)
00174  * This means effectively zeroing out A[i,i+2:N] (the entire right super-
00175  * diagonal of the i-th row of A, but ONE superdiag element) and replacing
00176  * A[i,i+1] with - Norm. Thus modified A[i,i+1] is returned as the result of
00177  * the present function.
00178  * The other (i+1:M) rows of A are transformed as
00179  *    (4)  A[j,] * VP = A[j,] - VPi' * ( A[j,] * VPi / beta )
00180  * Note, due to (2), only elements of columns i+1:N actually  participate
00181  * in above transforms; the left i columns of A are not affected.
00182  * As was mentioned earlier,
00183  * (5)  beta = A[i,] * VPi = (A[i,i+1] + Norm)*A[i,i+1]
00184  *                         + A[i,i+2:N]*A[i,i+2:N]
00185  *      = ||A[i,i+1:N]||^2 + Norm*A[i,i+1] = Norm^2 + Norm*A[i,i+1]
00186  * (note the sign of the Norm is the same as A[i,i+1])
00187  * For extra precision, vector VPi (and so is Norm and beta) are scaled,
00188  * which would not affect (4) as easy to see.
00189  *
00190  * The result of consecutive transformations (1) over matrix A is accumulated
00191  * in matrix V (which is initialized to be a unit matrix). At each step,
00192  * V is right-multiplied by VP. That is, rows of V are transformed similarly
00193  * to rows of A, see eq. (4). We also keep in mind that multiplication by
00194  * VP at the i-th step does not affect the first i rows of V.
00195  * Note that vector VPi doesn't have to be allocated explicitly: its
00196  * first i components are zeros (which we can always imply in computations),
00197  * and the rest of the components (but the VPi[i+1]) are the same as those
00198  * of A[i,i+1:N], the superdiagonal of A[i,]. This row, A[i,] is affected
00199  * only trivially as explained above, that is, we don't need to carry this
00200  * transformation explicitly (only A[i,i+1] is going to be non-trivially
00201  * affected, that is, replaced by -Norm, but we will use super_diag[i+1] to
00202  * store the result).
00203  *
00204  */
00205  
00206 inline double SVD::right_householder(Matrix& A, const int i)
00207 {
00208    IRange range = IRange::from(i+1 - A.q_row_lwb());// Note only VPi[i+1:N] matter
00209    MatrixRow y(A,i);
00210    LAStrideStreamOut VPi_str(y,range);
00211    register int j;
00212    REAL scale = 0;                      // Compute the scaling factor
00213    while( !VPi_str.eof() )
00214      scale += fabs(VPi_str.get());
00215    if( scale == 0 )                     // If A[i,] is a null vector, no
00216      return 0;                          // transform is required
00217  
00218    VPi_str.rewind();
00219    double Norm_sqr = 0;                 // Scale VPi (that is, A[i,])
00220    while( !VPi_str.eof() )              // and compute its norm, Norm^2
00221      Norm_sqr += sqr(VPi_str.get() /= scale);
00222    VPi_str.rewind();
00223    double new_Aii1 = sqrt(Norm_sqr);    // new_Aii1 = -Norm, Norm has the
00224    if( VPi_str.peek() > 0 )             // same sign as
00225      new_Aii1 = -new_Aii1;              // Aii1 (that is, VPi[i+1])
00226    const float beta = - VPi_str.peek()*new_Aii1 + Norm_sqr;
00227    VPi_str.peek() -= new_Aii1;          // VPi[i+1] = A[i,i+1] - (-Norm)
00228    
00229    for(j=i+1; j<=M; j++)        // Transform i+1:M rows of A
00230    {
00231      MatrixRow yy(A,j);
00232      LAStrideStreamOut Aj_str( yy,range );
00233      REAL factor = 0;
00234      while( !VPi_str.eof() )
00235        factor += VPi_str.get() * Aj_str.get();  // Compute A[j,] * VPi
00236      factor /= beta;
00237      for(VPi_str.rewind(), Aj_str.rewind(); !VPi_str.eof(); )
00238        Aj_str.get() -= VPi_str.get() * factor;
00239      VPi_str.rewind();
00240    }
00241 
00242    for(j=1; j<=N; j++)  // Accumulate the transform in V
00243    {
00244      MatrixRow yy(V,j);
00245      LAStrideStreamOut Vj_str( yy, range );
00246      REAL factor = 0;
00247      while( !VPi_str.eof() )
00248        factor += VPi_str.get() * Vj_str.get();  // Compute  V[j,] * VPi
00249      factor /= beta;
00250      for(VPi_str.rewind(), Vj_str.rewind(); !VPi_str.eof(); )
00251        Vj_str.get() -= VPi_str.get() * factor;
00252      VPi_str.rewind();
00253    }
00254    return new_Aii1 * scale;             // Scale new Aii1 back
00255 }
00256 
00257 /*
00258  *------------------------------------------------------------------------
00259  *                        Bidiagonalization
00260  * This nethod turns matrix A into a bidiagonal one. Its N diagonal elements
00261  * are stored in a vector sig, while N-1 superdiagonal elements are stored
00262  * in a vector super_diag(2:N) (with super_diag(1) being always 0).
00263  * Matrices U and V store the record of orthogonal Householder
00264  * reflections that were used to convert A to this form. The method
00265  * returns the norm of the resulting bidiagonal matrix, that is, the
00266  * maximal column sum.
00267  */
00268 
00269 double SVD::bidiagonalize(Vector& super_diag, const Matrix& __A)
00270 {
00271   double norm_acc = 0;
00272   super_diag(1) = 0;                    // No superdiag elem above A(1,1)
00273   Matrix A = __A;                       // A being transformed
00274   A.resize_to(__A.q_nrows(),__A.q_ncols()); // Indexing from 1
00275   for(register int i=1; i<=N; i++)
00276   {
00277     const REAL& diagi = sig(i) = left_householder(A,i);
00278     if( i < N )
00279       super_diag(i+1) = right_householder(A,i);
00280     norm_acc = max(norm_acc,(double)fabs(diagi)+fabs(super_diag(i)));
00281   }
00282   return norm_acc;
00283 }
00284 
00285 /*
00286  *------------------------------------------------------------------------
00287  *              QR-diagonalization of a bidiagonal matrix
00288  *
00289  * After bidiagonalization we get a bidiagonal matrix J:
00290  *    (1)  J = U' * A * V
00291  * The present method turns J into a matrix JJ by applying a set of
00292  * orthogonal transforms
00293  *    (2)  JJ = S' * J * T
00294  * Orthogonal matrices S and T are chosen so that JJ were also a
00295  * bidiagonal matrix, but with superdiag elements smaller than those of J.
00296  * We repeat (2) until non-diag elements of JJ become smaller than EPS
00297  * and can be disregarded.
00298  * Matrices S and T are constructed as
00299  *    (3)  S = S1 * S2 * S3 ... Sn, and similarly T
00300  * where Sk and Tk are matrices of simple rotations
00301  *    (4)  Sk[i,j] = i==j ? 1 : 0 for all i>k or i<k-1
00302  *         Sk[k-1,k-1] = cos(Phk),  Sk[k-1,k] = -sin(Phk),
00303  *         SK[k,k-1] = sin(Phk),    Sk[k,k] = cos(Phk), k=2..N
00304  * Matrix Tk is constructed similarly, only with angle Thk rather than Phk.
00305  *
00306  * Thus left multiplication of J by SK' can be spelled out as
00307  *    (5)  (Sk' * J)[i,j] = J[i,j] when i>k or i<k-1,
00308  *                  [k-1,j] = cos(Phk)*J[k-1,j] + sin(Phk)*J[k,j]
00309  *                  [k,j] =  -sin(Phk)*J[k-1,j] + cos(Phk)*J[k,j]
00310  * That is, k-1 and k rows of J are replaced by their linear combinations;
00311  * the rest of J is unaffected. Right multiplication of J by Tk similarly
00312  * changes only k-1 and k columns of J.
00313  * Matrix T2 is chosen the way that T2'J'JT2 were a QR-transform with a
00314  * shift. Note that multiplying J by T2 gives rise to a J[2,1] element of
00315  * the product J (which is below the main diagonal). Angle Ph2 is then
00316  * chosen so that multiplication by S2' (which combines 1 and 2 rows of J)
00317  * gets rid of that elemnent. But this will create a [1,3] non-zero element.
00318  * T3 is made to make it disappear, but this leads to [3,2], etc.
00319  * In the end, Sn removes a [N,N-1] element of J and matrix S'JT becomes
00320  * bidiagonal again. However, because of a special choice
00321  * of T2 (QR-algorithm), its non-diag elements are smaller than those of J.
00322  *
00323  * All this process in more detail is described in
00324  *    J.H. Wilkinson, C. Reinsch. Linear algebra - Springer-Verlag, 1971
00325  *
00326  * If during transforms (1), JJ[N-1,N] turns 0, then JJ[N,N] is a singular
00327  * number (possibly with a wrong (that is, negative) sign). This is a
00328  * consequence of Frantsis' Theorem, see the book above. In that case, we can
00329  * eliminate the N-th row and column of JJ and carry out further transforms
00330  * with a smaller matrix. If any other superdiag element of JJ turns 0,
00331  * then JJ effectively falls into two independent matrices. We will process
00332  * them independently (the bottom one first).
00333  *
00334  * Since matrix J is a bidiagonal, it can be stored efficiently. As a matter
00335  * of fact, its N diagonal elements are in array Sig, and superdiag elements
00336  * are stored in array super_diag.
00337  */
00338  
00339                                 // Carry out U * S with a rotation matrix
00340                                 // S (which combines i-th and j-th columns
00341                                 // of U, i>j)
00342 inline void SVD::rotate(Matrix& U, const int i, const int j,
00343                    const double cos_ph, const double sin_ph)
00344 {
00345   MatrixColumn x(U,i);
00346   LAStreamOut Ui(x);
00347   MatrixColumn y(U,j);
00348   LAStreamOut Uj(y);
00349   while( !Ui.eof() )
00350   {
00351     REAL& Uil = Ui.get(); REAL& Ujl = Uj.get();
00352     const REAL Ujl_was = Ujl;
00353     Ujl =  cos_ph*Ujl_was + sin_ph*Uil;
00354     Uil = -sin_ph*Ujl_was + cos_ph*Uil;
00355   }
00356 }
00357 
00358 /*
00359  * A diagonal element J[l-1,l-1] turns out 0 at the k-th step of the
00360  * algorithm. That means that one of the original matrix' singular numbers
00361  * shall be zero. In that case, we multiply J by specially selected
00362  * matrices S' of simple rotations to eliminate a superdiag element J[l-1,l].
00363  * After that, matrix J falls into two pieces, which can be dealt with
00364  * in a regular way (the bottom piece first).
00365  * 
00366  * These special S transformations are accumulated into matrix U: since J
00367  * is left-multiplied by S', U would be right-multiplied by S. Transform
00368  * formulas for doing these rotations are similar to (5) above. See the
00369  * book cited above for more details.
00370  */
00371 inline void SVD::rip_through(
00372         Vector& super_diag, const int k, const int l, const double eps)
00373 {
00374   double cos_ph = 0, sin_ph = 1;        // Accumulate cos,sin of Ph
00375                                 // The first step of the loop below
00376                                 // (when i==l) would eliminate J[l-1,l],
00377                                 // which is stored in super_diag(l)
00378                                 // However, it gives rise to J[l-1,l+1]
00379                                 // and J[l,l+2]
00380                                 // The following steps eliminate these
00381                                 // until they fall below
00382                                 // significance
00383   for(register int i=l; i<=k; i++)
00384   {
00385     const double f = sin_ph * super_diag(i);
00386     super_diag(i) *= cos_ph;
00387     if( fabs(f) <= eps )
00388       break;                    // Current J[l-1,l] became unsignificant
00389     cos_ph = sig(i); sin_ph = -f;       // unnormalized sin/cos
00390     const double norm = (sig(i) = hypot(cos_ph,sin_ph)); // sqrt(sin^2+cos^2)
00391     cos_ph /= norm, sin_ph /= norm;     // Normalize sin/cos
00392     rotate(U,i,l-1,cos_ph,sin_ph);
00393   }
00394 }
00395 
00396                         // We're about to proceed doing QR-transforms
00397                         // on a (bidiag) matrix J[1:k,1:k]. It may happen
00398                         // though that the matrix splits (or can be
00399                         // split) into two independent pieces. This function
00400                         // checks for splitting and returns the lowerbound
00401                         // index l of the bottom piece, J[l:k,l:k]
00402 inline int SVD::get_submatrix_to_work_on(
00403         Vector& super_diag, const int k, const double eps)
00404 {
00405   for(register int l=k; l>1; l--)
00406     if( fabs(super_diag(l)) <= eps )
00407       return l;                         // The breaking point: zero J[l-1,l]
00408     else if( fabs(sig(l-1)) <= eps )    // Diagonal J[l,l] turns out 0
00409     {                                   // meaning J[l-1,l] _can_ be made
00410       rip_through(super_diag,k,l,eps);  // zero after some rotations
00411       return l;
00412     }
00413   return 1;                     // Deal with J[1:k,1:k] as a whole
00414 }
00415 
00416                 // Diagonalize root module
00417 void SVD::diagonalize(Vector& super_diag, const double eps)
00418 {
00419   for(register int k=N; k>=1; k--)      // QR-iterate upon J[l:k,l:k]
00420   {
00421     register int l;
00422     while(l=get_submatrix_to_work_on(super_diag,k,eps),
00423           fabs(super_diag(k)) > eps)    // until superdiag J[k-1,k] becomes 0
00424     {
00425       double shift;                     // Compute a QR-shift from a bottom
00426       {                                 // corner minor of J[l:k,l:k] order 2
00427         REAL Jk2k1 = super_diag(k-1),   // J[k-2,k-1]
00428              Jk1k  = super_diag(k),
00429              Jk1k1 = sig(k-1),          // J[k-1,k-1]
00430              Jkk   = sig(k),
00431              Jll   = sig(l);            // J[l,l]
00432         shift = (Jk1k1-Jkk)*(Jk1k1+Jkk) + (Jk2k1-Jk1k)*(Jk2k1+Jk1k);
00433         shift /= 2*Jk1k*Jk1k1;
00434         shift += (shift < 0 ? -1 : 1) * sqrt(shift*shift+1);
00435         shift = ( (Jll-Jkk)*(Jll+Jkk) + Jk1k*(Jk1k1/shift-Jk1k) )/Jll;
00436       }
00437                                 // Carry on multiplications by T2, S2, T3...
00438       double cos_th = 1, sin_th = 1;
00439       REAL Ji1i1 = sig(l);      // J[i-1,i-1] at i=l+1...k
00440       for(register int i=l+1; i<=k; i++)
00441       {
00442         REAL Ji1i = super_diag(i), Jii = sig(i);  // J[i-1,i] and J[i,i]
00443         sin_th *= Ji1i; Ji1i *= cos_th; cos_th = shift;
00444         double norm_f = (super_diag(i-1) = hypot(cos_th,sin_th));
00445         cos_th /= norm_f, sin_th /= norm_f;
00446                                         // Rotate J[i-1:i,i-1:i] by Ti
00447         shift = cos_th*Ji1i1 + sin_th*Ji1i;     // new J[i-1,i-1]
00448         Ji1i = -sin_th*Ji1i1 + cos_th*Ji1i;     // J[i-1,i] after rotation
00449         const double Jii1 = Jii*sin_th;         // Emerged J[i,i-1]
00450         Jii *= cos_th;                          // new J[i,i]
00451         rotate(V,i,i-1,cos_th,sin_th); // Accumulate T rotations in V
00452         
00453         double cos_ph = shift, sin_ph = Jii1;// Make Si to get rid of J[i,i-1]
00454         sig(i-1) = (norm_f = hypot(cos_ph,sin_ph));     // New J[i-1,i-1]
00455         if( norm_f == 0 )               // If norm =0, rotation angle
00456           cos_ph = cos_th, sin_ph = sin_th; // can be anything now
00457         else
00458           cos_ph /= norm_f, sin_ph /= norm_f;
00459                                         // Rotate J[i-1:i,i-1:i] by Si
00460         shift = cos_ph * Ji1i + sin_ph*Jii;     // New J[i-1,i]
00461         Ji1i1 = -sin_ph*Ji1i + cos_ph*Jii;      // New Jii, would carry over
00462                                                 // as J[i-1,i-1] for next i
00463         rotate(U,i,i-1,cos_ph,sin_ph);  // Accumulate S rotations in U
00464                                         // Jii1 disappears, sin_th would
00465         cos_th = cos_ph, sin_th = sin_ph; // carry over a (scaled) J[i-1,i+1]
00466                                         // to eliminate on the next i, cos_th
00467                                         // would carry over a scaled J[i,i+1]
00468       }
00469       super_diag(l) = 0;                // Supposed to be eliminated by now
00470       super_diag(k) = shift;
00471       sig(k) = Ji1i1;
00472     }           // --- end-of-QR-iterations
00473     if( sig(k) < 0 )            // Correct the sign of the sing number
00474     {
00475       sig(k) = -sig(k);
00476       MatrixColumn y(V,k);
00477       for(LAStreamOut Vk(y); !Vk.eof(); )
00478         { REAL& vk = Vk.get(); vk = -vk; }
00479     }
00480   }
00481     
00482 } 
00483 
00484 
00485 /*
00486  *------------------------------------------------------------------------
00487  *                              The root Module
00488  */
00489 
00490 SVD::SVD(const Matrix& A)
00491    : M(A.q_nrows()), N(A.q_ncols()),
00492      U(A.q_nrows(),A.q_nrows()),
00493      V(A.q_ncols(),A.q_ncols()),
00494      sig(A.q_ncols())
00495 {
00496   if( M < N )
00497     A.info(),
00498     _error("Matrix A should have at least as many rows as it has columns");
00499      
00500   U.unit_matrix(); V.unit_matrix();
00501 
00502   Vector super_diag(N);
00503   const double bidiag_norm = bidiagonalize(super_diag,A);
00504   const double eps = FLT_EPSILON * bidiag_norm; // Significance threshold
00505   diagonalize(super_diag,eps);
00506 }
00507 
00508 /*
00509  *------------------------------------------------------------------------
00510  *              Print some info about the SVD that has been built
00511  */
00512 
00513                                 // Print the info about the SVD
00514 void SVD::info(void) const
00515 {
00516   U.is_valid();
00517   message("\nSVD of an %dx%d matrix",M,N);
00518 }
00519 
00520                                 // Return min and max singular values
00521 SVD::operator MinMax(void) const
00522 {
00523   LAStreamIn sigs(sig);
00524   MinMax mm(sigs.get());
00525   while( !sigs.eof() )
00526     mm << sigs.get();
00527   return mm;
00528 }
00529 
00530                                         // sig_max/sig_min
00531 double SVD::q_cond_number(void) const
00532 {
00533   return ((MinMax)(*this)).ratio();
00534 }
00535 
00536 /*
00537  *------------------------------------------------------------------------
00538  *                      class SVD_inv_mult
00539  *      Solving an (overspecified) set of linear equations A*X=B
00540  *              using a least squares method via SVD
00541  *
00542  * Matrix A is a rectangular M*N matrix with M>=N, B is a M*K matrix
00543  * with K either =1 (that is, B is merely a column-vector) or K>1.
00544  * If B is a unit matrix of the size that of A, the present LazyMatrix is
00545  * a regularized (pseudo)inverse of A.
00546  *
00547  * Algorithm
00548  *   Matrix A is decomposed first using SVD:
00549  *      (1) A = U*Sig*V'
00550  * where matrices U and V are orthogonal and Sig is diagonal.
00551  * The set of simultaneous linear equations AX=B can be written then as
00552  *      (2) Sig*V'*X = U'*B
00553  * (where we have used the fact that U'=inv(U)), or
00554  *      (3) Sig*Vx = Ub
00555  * where we introduced
00556  *      (4) Vx = V'*X and Ub = U'*b
00557  * Since Sig is a diag matrix, eq. (3) is solved trivially:
00558  *      (5) Vx[i] = Ub[i]/sig[i], if sig[i] > tau
00559  *                = 0   otherwise
00560  * In the latter case matrix Sig has an incomplete rank, or close to that.
00561  * That is, one or more singular values are _too_ small. This means that
00562  * there is linear dependence among the equations of the set. In that case,
00563  * we will print the corresponding "null coefficients", the corresponding
00564  * columns of V. Adding them to the solution X won't change A*X.
00565  *
00566  * Having computed Vx, X is simply recovered as
00567  *      (6) X = V*Vx
00568  * since V*V' = E.
00569  *
00570  * Threshold tau in (5) is chosen as N*FLT_EPSILON*max(Sig[i,i]) unless
00571  * specified otherwise by the user.
00572  */
00573 
00574 SVD_inv_mult::SVD_inv_mult
00575         (const SVD& _svd, const Matrix& __B,const double _tau)
00576         : LazyMatrix(_svd.q_V().q_ncols(),__B.q_ncols()),
00577           svd(_svd), B(__B), tau(_tau)
00578 {
00579   if( svd.q_U().q_nrows() != B.q_nrows() )
00580     svd.info(), B.info(),
00581     _error("Unsuitable matrices for SVD*X=B set");
00582   MinMax sig_minmax = (MinMax)svd;
00583   if( tau == 0 )
00584    tau = svd.q_V().q_nrows() * sig_minmax.max() * FLT_EPSILON;
00585   are_zero_coeff = sig_minmax.min() < tau;
00586 }
00587 
00588 #if 0
00589                                 // Computing X in a special case where
00590                                 // B (and X) is a vector
00591 inline void SVD_inv_mult::fill_in_vector(Matrix& X) const
00592 {
00593   X.clear();
00594   bool x_is_vector = x.q_ncols() == 1;
00595   if( are_zero_coeff )
00596     message("\nSVD solver of AX=B detected a linear dependency among X"
00597             "\n  #  \tsingular value\tnull coefficients\n");
00598   for(register int i=1; i<=svd.q_V().q_nrows(); i++)
00599   {
00600      MatrixCol Ui(U,i);
00601      const double sigi = svd.q_sig()(i);
00602      if( sigi > tau )
00603        if( x_is_vector )
00604          X(i) = (Ui * B)/sigi;
00605        else
00606          X(i) =0, print null coeffs...
00607   }
00608   X *= V;                       // ???
00609 }
00610 #endif
00611                                 // Computing X in a general case where
00612                                 // B (and X) is a rectangular matrix
00613 void SVD_inv_mult::fill_in(Matrix& X) const
00614 {
00615   if( are_zero_coeff )
00616     message("\nSVD solver of AX=B detected a linear dependency among X"
00617             "\n  #  singular value\tnull coefficients\n");
00618   Matrix Vx = zero(X);
00619   Vector Ui(svd.q_U().q_nrows());                       // i-th col of U
00620   Vector Bj(B.q_nrows());                               // j-th col of B
00621   const Matrix& V = svd.q_V();
00622   const Matrix& U = svd.q_U();
00623   LAStreamIn sigs(svd.q_sig());
00624   for(register int i=1; i<=V.q_ncols(); i++)
00625   {
00626     const double sigi = sigs.get();
00627     MatrixRow y(Vx,i);
00628     LAStrideStreamOut Vxi( y );
00629     if( sigi > tau )
00630     {
00631       ConstMatrixColumn x(U,i);
00632       to_every(Ui) = of_every(x);
00633       for(register int j=1; j<=X.q_ncols(); j++) {
00634         ConstMatrixColumn yy(B,j);
00635         to_every(Bj) = of_every(yy), Vxi.get() = Ui * Bj/sigi;
00636       }
00637     }
00638     else
00639     {
00640       while( !Vxi.eof() ) Vxi.get() = 0;
00641       message(" %d %12.2g \t(",i,sigi);
00642       ConstMatrixColumn x(V,i);
00643       LAStreamIn Vi(x);
00644       while( !Vi.eof() )
00645         message("%10g ",Vi.get());
00646       message(")\n");
00647     }
00648   }
00649   X.mult(V,Vx);                         // that is, set X to be V*Vx
00650 }
00651  
00652 }

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