[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 2007 by Ullrich Koethe */ 00004 /* */ 00005 /* This file is part of the VIGRA computer vision library. */ 00006 /* ( Version 1.6.0, Aug 13 2008 ) */ 00007 /* The VIGRA Website is */ 00008 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ 00009 /* Please direct questions, bug reports, and contributions to */ 00010 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00011 /* vigra@informatik.uni-hamburg.de */ 00012 /* */ 00013 /* Permission is hereby granted, free of charge, to any person */ 00014 /* obtaining a copy of this software and associated documentation */ 00015 /* files (the "Software"), to deal in the Software without */ 00016 /* restriction, including without limitation the rights to use, */ 00017 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00018 /* sell copies of the Software, and to permit persons to whom the */ 00019 /* Software is furnished to do so, subject to the following */ 00020 /* conditions: */ 00021 /* */ 00022 /* The above copyright notice and this permission notice shall be */ 00023 /* included in all copies or substantial portions of the */ 00024 /* Software. */ 00025 /* */ 00026 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00027 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00028 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00029 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00030 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00031 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00032 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00033 /* OTHER DEALINGS IN THE SOFTWARE. */ 00034 /* */ 00035 /************************************************************************/ 00036 00037 00038 #ifndef VIGRA_SINGULAR_VALUE_DECOMPOSITION_HXX 00039 #define VIGRA_SINGULAR_VALUE_DECOMPOSITION_HXX 00040 00041 #include "matrix.hxx" 00042 #include "array_vector.hxx" 00043 00044 00045 namespace vigra 00046 { 00047 00048 namespace linalg 00049 { 00050 00051 /** Singular Value Decomposition. 00052 \ingroup MatrixAlgebra 00053 00054 For an m-by-n matrix \a A with m >= n, the singular value decomposition is 00055 an m-by-n orthogonal matrix \a U, an n-by-n diagonal matrix S, and 00056 an n-by-n orthogonal matrix \a V so that A = U*S*V'. 00057 00058 To save memory, this functions stores the matrix \a S in a column vector of 00059 appropriate length (a diagonal matrix can be obtained by <tt>diagonalMatrix(S)</tt>). 00060 The singular values, sigma[k] = S(k, 0), are ordered so that 00061 sigma[0] >= sigma[1] >= ... >= sigma[n-1]. 00062 00063 The singular value decomposition always exists, so this function will 00064 never fail (except if the shapes of the argument matrices don't match). 00065 The effective numerical rank of A is returned. 00066 00067 (Adapted from JAMA, a Java Matrix Library, developed jointly 00068 by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama). 00069 00070 <b>\#include</b> <<a href="singular__value__decomposition_8hxx-source.html">vigra/singular_value_decomposition.hxx</a>> or<br> 00071 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00072 Namespaces: vigra and vigra::linalg 00073 */ 00074 template <class T, class C1, class C2, class C3, class C4> 00075 unsigned int 00076 singularValueDecomposition(MultiArrayView<2, T, C1> const & A, 00077 MultiArrayView<2, T, C2> &U, MultiArrayView<2, T, C3> &S, MultiArrayView<2, T, C4> &V) 00078 { 00079 typedef T Real; 00080 typedef MultiArrayShape<2>::type Shape; 00081 00082 const MultiArrayIndex rows = rowCount(A); 00083 const MultiArrayIndex cols = columnCount(A); 00084 vigra_precondition(rows >= cols, 00085 "singularValueDecomposition(): Input matrix A must be rectangular with rowCount >= columnCount."); 00086 vigra_precondition(rowCount(S) == cols && columnCount(S) == 1, 00087 "singularValueDecomposition(): Output S must be column vector with rowCount == columnCount(A)."); 00088 vigra_precondition(rowCount(U) == rows && columnCount(U) == cols, 00089 "singularValueDecomposition(): Output matrix U must have the same dimensions as input matrix A."); 00090 vigra_precondition(rowCount(V) == cols && columnCount(V) == cols, 00091 "singularValueDecomposition(): Output matrix V must be square with n = columnCount(A)."); 00092 00093 MultiArrayIndex m = rows; 00094 MultiArrayIndex n = cols; 00095 MultiArrayIndex nu = n; 00096 00097 U.init(0.0); 00098 S.init(0.0); 00099 V.init(0.0); 00100 00101 ArrayVector<Real> e((unsigned int)n); 00102 ArrayVector<Real> work((unsigned int)m); 00103 Matrix<Real> a(A); 00104 MultiArrayView<1, T, C3> s = S.bindOuter(0); 00105 00106 MultiArrayIndex i=0, j=0, k=0; 00107 00108 // Reduce a to bidiagonal form, storing the diagonal elements 00109 // in s and the super-diagonal elements in e. 00110 MultiArrayIndex nct = std::min(m-1,n); 00111 MultiArrayIndex nrt = std::max((MultiArrayIndex)0,n-2); 00112 for (k = 0; k < std::max(nct,nrt); ++k) 00113 { 00114 if (k < nct) 00115 { 00116 // Compute the transformation for the k-th column and 00117 // place the k-th diagonal in s(k). 00118 // Compute 2-norm of k-th column without under/overflow. 00119 s(k) = 0.0; 00120 for (i = k; i < m; ++i) 00121 { 00122 s(k) = hypot(s(k), a(i, k)); 00123 } 00124 if (s(k) != 0.0) 00125 { 00126 if (a(k, k) < 0.0) 00127 { 00128 s(k) = -s(k); 00129 } 00130 for (i = k; i < m; ++i) 00131 { 00132 a(i, k) /= s(k); 00133 } 00134 a(k, k) += 1.0; 00135 } 00136 s(k) = -s(k); 00137 } 00138 for (j = k+1; j < n; ++j) 00139 { 00140 if ((k < nct) && (s(k) != 0.0)) 00141 { 00142 // Apply the transformation. 00143 Real t(0.0); 00144 for (i = k; i < m; ++i) 00145 { 00146 t += a(i, k)*a(i, j); 00147 } 00148 t = -t/a(k, k); 00149 for (i = k; i < m; ++i) 00150 { 00151 a(i, j) += t*a(i, k); 00152 } 00153 } 00154 00155 // Place the k-th row of a into e for the 00156 // subsequent calculation of the row transformation. 00157 00158 e[j] = a(k, j); 00159 } 00160 if (k < nct) 00161 { 00162 // Place the transformation in U for subsequent back 00163 // multiplication. 00164 00165 for (i = k; i < m; ++i) 00166 { 00167 U(i, k) = a(i, k); 00168 } 00169 } 00170 if (k < nrt) 00171 { 00172 // Compute the k-th row transformation and place the 00173 // k-th super-diagonal in e[k]. 00174 // Compute 2-norm without under/overflow. 00175 e[k] = 0; 00176 for (i = k+1; i < n; ++i) 00177 { 00178 e[k] = hypot(e[k],e[i]); 00179 } 00180 if (e[k] != 0.0) 00181 { 00182 if (e[k+1] < 0.0) 00183 { 00184 e[k] = -e[k]; 00185 } 00186 for (i = k+1; i < n; ++i) 00187 { 00188 e[i] /= e[k]; 00189 } 00190 e[k+1] += 1.0; 00191 } 00192 e[k] = -e[k]; 00193 if ((k+1 < m) & (e[k] != 0.0)) 00194 { 00195 // Apply the transformation. 00196 for (i = k+1; i < m; ++i) 00197 { 00198 work[i] = 0.0; 00199 } 00200 for (j = k+1; j < n; ++j) 00201 { 00202 for (i = k+1; i < m; ++i) 00203 { 00204 work[i] += e[j]*a(i, j); 00205 } 00206 } 00207 for (j = k+1; j < n; ++j) 00208 { 00209 Real t(-e[j]/e[k+1]); 00210 for (i = k+1; i < m; ++i) 00211 { 00212 a(i, j) += t*work[i]; 00213 } 00214 } 00215 } 00216 // Place the transformation in V for subsequent 00217 // back multiplication. 00218 for (i = k+1; i < n; ++i) 00219 { 00220 V(i, k) = e[i]; 00221 } 00222 } 00223 } 00224 00225 // Set up the final bidiagonal matrix of order p. 00226 00227 MultiArrayIndex p = n; 00228 if (nct < n) 00229 { 00230 s(nct) = a(nct, nct); 00231 } 00232 if (m < p) 00233 { 00234 s(p-1) = 0.0; 00235 } 00236 if (nrt+1 < p) 00237 { 00238 e[nrt] = a(nrt, p-1); 00239 } 00240 e[p-1] = 0.0; 00241 00242 // Generate U. 00243 for (j = nct; j < nu; ++j) 00244 { 00245 for (i = 0; i < m; ++i) 00246 { 00247 U(i, j) = 0.0; 00248 } 00249 U(j, j) = 1.0; 00250 } 00251 for (k = nct-1; k >= 0; --k) 00252 { 00253 if (s(k) != 0.0) 00254 { 00255 for (j = k+1; j < nu; ++j) 00256 { 00257 Real t(0.0); 00258 for (i = k; i < m; ++i) 00259 { 00260 t += U(i, k)*U(i, j); 00261 } 00262 t = -t/U(k, k); 00263 for (i = k; i < m; ++i) 00264 { 00265 U(i, j) += t*U(i, k); 00266 } 00267 } 00268 for (i = k; i < m; ++i ) 00269 { 00270 U(i, k) = -U(i, k); 00271 } 00272 U(k, k) = 1.0 + U(k, k); 00273 for (i = 0; i < k-1; ++i) 00274 { 00275 U(i, k) = 0.0; 00276 } 00277 } 00278 else 00279 { 00280 for (i = 0; i < m; ++i) 00281 { 00282 U(i, k) = 0.0; 00283 } 00284 U(k, k) = 1.0; 00285 } 00286 } 00287 00288 // Generate V. 00289 for (k = n-1; k >= 0; --k) 00290 { 00291 if ((k < nrt) & (e[k] != 0.0)) 00292 { 00293 for (j = k+1; j < nu; ++j) 00294 { 00295 Real t(0.0); 00296 for (i = k+1; i < n; ++i) 00297 { 00298 t += V(i, k)*V(i, j); 00299 } 00300 t = -t/V(k+1, k); 00301 for (i = k+1; i < n; ++i) 00302 { 00303 V(i, j) += t*V(i, k); 00304 } 00305 } 00306 } 00307 for (i = 0; i < n; ++i) 00308 { 00309 V(i, k) = 0.0; 00310 } 00311 V(k, k) = 1.0; 00312 } 00313 00314 // Main iteration loop for the singular values. 00315 00316 MultiArrayIndex pp = p-1; 00317 int iter = 0; 00318 Real eps = NumericTraits<Real>::epsilon()*2.0; 00319 while (p > 0) 00320 { 00321 MultiArrayIndex k=0; 00322 int kase=0; 00323 00324 // Here is where a test for too many iterations would go. 00325 00326 // This section of the program inspects for 00327 // negligible elements in the s and e arrays. On 00328 // completion the variables kase and k are set as follows. 00329 00330 // kase = 1 if s(p) and e[k-1] are negligible and k<p 00331 // kase = 2 if s(k) is negligible and k<p 00332 // kase = 3 if e[k-1] is negligible, k<p, and 00333 // s(k), ..., s(p) are not negligible (qr step). 00334 // kase = 4 if e(p-1) is negligible (convergence). 00335 00336 for (k = p-2; k >= -1; --k) 00337 { 00338 if (k == -1) 00339 { 00340 break; 00341 } 00342 if (abs(e[k]) <= eps*(abs(s(k)) + abs(s(k+1)))) 00343 { 00344 e[k] = 0.0; 00345 break; 00346 } 00347 } 00348 if (k == p-2) 00349 { 00350 kase = 4; 00351 } 00352 else 00353 { 00354 MultiArrayIndex ks; 00355 for (ks = p-1; ks >= k; --ks) 00356 { 00357 if (ks == k) 00358 { 00359 break; 00360 } 00361 Real t( (ks != p ? abs(e[ks]) : 0.) + 00362 (ks != k+1 ? abs(e[ks-1]) : 0.)); 00363 if (abs(s(ks)) <= eps*t) 00364 { 00365 s(ks) = 0.0; 00366 break; 00367 } 00368 } 00369 if (ks == k) 00370 { 00371 kase = 3; 00372 } 00373 else if (ks == p-1) 00374 { 00375 kase = 1; 00376 } 00377 else 00378 { 00379 kase = 2; 00380 k = ks; 00381 } 00382 } 00383 ++k; 00384 00385 // Perform the task indicated by kase. 00386 00387 switch (kase) 00388 { 00389 case 1: // Deflate negligible s(p). 00390 { 00391 Real f(e[p-2]); 00392 e[p-2] = 0.0; 00393 for (j = p-2; j >= k; --j) 00394 { 00395 Real t( hypot(s(j),f)); 00396 Real cs(s(j)/t); 00397 Real sn(f/t); 00398 s(j) = t; 00399 if (j != k) 00400 { 00401 f = -sn*e[j-1]; 00402 e[j-1] = cs*e[j-1]; 00403 } 00404 for (i = 0; i < n; ++i) 00405 { 00406 t = cs*V(i, j) + sn*V(i, p-1); 00407 V(i, p-1) = -sn*V(i, j) + cs*V(i, p-1); 00408 V(i, j) = t; 00409 } 00410 } 00411 break; 00412 } 00413 case 2: // Split at negligible s(k). 00414 { 00415 Real f(e[k-1]); 00416 e[k-1] = 0.0; 00417 for (j = k; j < p; ++j) 00418 { 00419 Real t(hypot(s(j),f)); 00420 Real cs( s(j)/t); 00421 Real sn(f/t); 00422 s(j) = t; 00423 f = -sn*e[j]; 00424 e[j] = cs*e[j]; 00425 for (i = 0; i < m; ++i) 00426 { 00427 t = cs*U(i, j) + sn*U(i, k-1); 00428 U(i, k-1) = -sn*U(i, j) + cs*U(i, k-1); 00429 U(i, j) = t; 00430 } 00431 } 00432 break; 00433 } 00434 case 3: // Perform one qr step. 00435 { 00436 // Calculate the shift. 00437 Real scale = std::max(std::max(std::max(std::max( 00438 abs(s(p-1)),abs(s(p-2))),abs(e[p-2])), 00439 abs(s(k))),abs(e[k])); 00440 Real sp = s(p-1)/scale; 00441 Real spm1 = s(p-2)/scale; 00442 Real epm1 = e[p-2]/scale; 00443 Real sk = s(k)/scale; 00444 Real ek = e[k]/scale; 00445 Real b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0; 00446 Real c = (sp*epm1)*(sp*epm1); 00447 Real shift = 0.0; 00448 if ((b != 0.0) || (c != 0.0)) 00449 { 00450 shift = VIGRA_CSTD::sqrt(b*b + c); 00451 if (b < 0.0) 00452 { 00453 shift = -shift; 00454 } 00455 shift = c/(b + shift); 00456 } 00457 Real f = (sk + sp)*(sk - sp) + shift; 00458 Real g = sk*ek; 00459 00460 // Chase zeros. 00461 for (j = k; j < p-1; ++j) 00462 { 00463 Real t = hypot(f,g); 00464 Real cs = f/t; 00465 Real sn = g/t; 00466 if (j != k) 00467 { 00468 e[j-1] = t; 00469 } 00470 f = cs*s(j) + sn*e[j]; 00471 e[j] = cs*e[j] - sn*s(j); 00472 g = sn*s(j+1); 00473 s(j+1) = cs*s(j+1); 00474 for (i = 0; i < n; ++i) 00475 { 00476 t = cs*V(i, j) + sn*V(i, j+1); 00477 V(i, j+1) = -sn*V(i, j) + cs*V(i, j+1); 00478 V(i, j) = t; 00479 } 00480 t = hypot(f,g); 00481 cs = f/t; 00482 sn = g/t; 00483 s(j) = t; 00484 f = cs*e[j] + sn*s(j+1); 00485 s(j+1) = -sn*e[j] + cs*s(j+1); 00486 g = sn*e[j+1]; 00487 e[j+1] = cs*e[j+1]; 00488 if (j < m-1) 00489 { 00490 for (i = 0; i < m; ++i) 00491 { 00492 t = cs*U(i, j) + sn*U(i, j+1); 00493 U(i, j+1) = -sn*U(i, j) + cs*U(i, j+1); 00494 U(i, j) = t; 00495 } 00496 } 00497 } 00498 e[p-2] = f; 00499 iter = iter + 1; 00500 break; 00501 } 00502 case 4: // Convergence. 00503 { 00504 // Make the singular values positive. 00505 if (s(k) <= 0.0) 00506 { 00507 s(k) = (s(k) < 0.0 ? -s(k) : 0.0); 00508 for (i = 0; i <= pp; ++i) 00509 { 00510 V(i, k) = -V(i, k); 00511 } 00512 } 00513 00514 // Order the singular values. 00515 00516 while (k < pp) 00517 { 00518 if (s(k) >= s(k+1)) 00519 { 00520 break; 00521 } 00522 Real t = s(k); 00523 s(k) = s(k+1); 00524 s(k+1) = t; 00525 if (k < n-1) 00526 { 00527 for (i = 0; i < n; ++i) 00528 { 00529 t = V(i, k+1); V(i, k+1) = V(i, k); V(i, k) = t; 00530 } 00531 } 00532 if (k < m-1) 00533 { 00534 for (i = 0; i < m; ++i) 00535 { 00536 t = U(i, k+1); U(i, k+1) = U(i, k); U(i, k) = t; 00537 } 00538 } 00539 ++k; 00540 } 00541 iter = 0; 00542 --p; 00543 break; 00544 } 00545 default: 00546 vigra_fail("vigra::svd(): Internal error."); 00547 } 00548 } 00549 Real tol = std::max(m,n)*s(0)*eps; 00550 unsigned int rank = 0; 00551 for (MultiArrayIndex i = 0; i < n; ++i) 00552 { 00553 if (s(i) > tol) 00554 { 00555 ++rank; 00556 } 00557 } 00558 return rank; // effective rank 00559 } 00560 00561 } // namespace linalg 00562 00563 using linalg::singularValueDecomposition; 00564 00565 } // namespace vigra 00566 00567 #endif // VIGRA_SINGULAR_VALUE_DECOMPOSITION_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|