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 }