vigra - vigra/singular_value_decomposition.hxx Source File

[ 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://hci.iwr.uni-heidelberg.de/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> <vigra/singular_value_decomposition.hxx> or<br>
    00070     <b>\#include</b> <vigra/linear_algebra.hxx><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.9.0 (Tue Nov 6 2012)