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