[ 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 /*    ( 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) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python VIGRA 1.6.0 (13 Aug 2008)