[ 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)