[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

vigra/singular_value_decomposition.hxx

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)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
VIGRA 1.6.0 (5 Nov 2009)