SVD.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra. Eigen itself is part of the KDE project.
00003 //
00004 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN2_SVD_H
00026 #define EIGEN2_SVD_H
00027 
00043 template<typename MatrixType> class SVD
00044 {
00045   private:
00046     typedef typename MatrixType::Scalar Scalar;
00047     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00048 
00049     enum {
00050       PacketSize = internal::packet_traits<Scalar>::size,
00051       AlignmentMask = int(PacketSize)-1,
00052       MinSize = EIGEN_SIZE_MIN_PREFER_DYNAMIC(MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime)
00053     };
00054 
00055     typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVector;
00056     typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> RowVector;
00057 
00058     typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MinSize> MatrixUType;
00059     typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixVType;
00060     typedef Matrix<Scalar, MinSize, 1> SingularValuesType;
00061 
00062   public:
00063 
00064     SVD() {} // a user who relied on compiler-generated default compiler reported problems with MSVC in 2.0.7
00065     
00066     SVD(const MatrixType& matrix)
00067       : m_matU(matrix.rows(), (std::min)(matrix.rows(), matrix.cols())),
00068         m_matV(matrix.cols(),matrix.cols()),
00069         m_sigma((std::min)(matrix.rows(),matrix.cols()))
00070     {
00071       compute(matrix);
00072     }
00073 
00074     template<typename OtherDerived, typename ResultType>
00075     bool solve(const MatrixBase<OtherDerived> &b, ResultType* result) const;
00076 
00077     const MatrixUType& matrixU() const { return m_matU; }
00078     const SingularValuesType& singularValues() const { return m_sigma; }
00079     const MatrixVType& matrixV() const { return m_matV; }
00080 
00081     void compute(const MatrixType& matrix);
00082     SVD& sort();
00083 
00084     template<typename UnitaryType, typename PositiveType>
00085     void computeUnitaryPositive(UnitaryType *unitary, PositiveType *positive) const;
00086     template<typename PositiveType, typename UnitaryType>
00087     void computePositiveUnitary(PositiveType *positive, UnitaryType *unitary) const;
00088     template<typename RotationType, typename ScalingType>
00089     void computeRotationScaling(RotationType *unitary, ScalingType *positive) const;
00090     template<typename ScalingType, typename RotationType>
00091     void computeScalingRotation(ScalingType *positive, RotationType *unitary) const;
00092 
00093   protected:
00095     MatrixUType m_matU;
00097     MatrixVType m_matV;
00099     SingularValuesType m_sigma;
00100 };
00101 
00106 template<typename MatrixType>
00107 void SVD<MatrixType>::compute(const MatrixType& matrix)
00108 {
00109   const int m = matrix.rows();
00110   const int n = matrix.cols();
00111   const int nu = (std::min)(m,n);
00112   ei_assert(m>=n && "In Eigen 2.0, SVD only works for MxN matrices with M>=N. Sorry!");
00113   ei_assert(m>1 && "In Eigen 2.0, SVD doesn't work on 1x1 matrices");
00114 
00115   m_matU.resize(m, nu);
00116   m_matU.setZero();
00117   m_sigma.resize((std::min)(m,n));
00118   m_matV.resize(n,n);
00119 
00120   RowVector e(n);
00121   ColVector work(m);
00122   MatrixType matA(matrix);
00123   const bool wantu = true;
00124   const bool wantv = true;
00125   int i=0, j=0, k=0;
00126 
00127   // Reduce A to bidiagonal form, storing the diagonal elements
00128   // in s and the super-diagonal elements in e.
00129   int nct = (std::min)(m-1,n);
00130   int nrt = (std::max)(0,(std::min)(n-2,m));
00131   for (k = 0; k < (std::max)(nct,nrt); ++k)
00132   {
00133     if (k < nct)
00134     {
00135       // Compute the transformation for the k-th column and
00136       // place the k-th diagonal in m_sigma[k].
00137       m_sigma[k] = matA.col(k).end(m-k).norm();
00138       if (m_sigma[k] != 0.0) // FIXME
00139       {
00140         if (matA(k,k) < 0.0)
00141           m_sigma[k] = -m_sigma[k];
00142         matA.col(k).end(m-k) /= m_sigma[k];
00143         matA(k,k) += 1.0;
00144       }
00145       m_sigma[k] = -m_sigma[k];
00146     }
00147 
00148     for (j = k+1; j < n; ++j)
00149     {
00150       if ((k < nct) && (m_sigma[k] != 0.0))
00151       {
00152         // Apply the transformation.
00153         Scalar t = matA.col(k).end(m-k).eigen2_dot(matA.col(j).end(m-k)); // FIXME dot product or cwise prod + .sum() ??
00154         t = -t/matA(k,k);
00155         matA.col(j).end(m-k) += t * matA.col(k).end(m-k);
00156       }
00157 
00158       // Place the k-th row of A into e for the
00159       // subsequent calculation of the row transformation.
00160       e[j] = matA(k,j);
00161     }
00162 
00163     // Place the transformation in U for subsequent back multiplication.
00164     if (wantu & (k < nct))
00165       m_matU.col(k).end(m-k) = matA.col(k).end(m-k);
00166 
00167     if (k < nrt)
00168     {
00169       // Compute the k-th row transformation and place the
00170       // k-th super-diagonal in e[k].
00171       e[k] = e.end(n-k-1).norm();
00172       if (e[k] != 0.0)
00173       {
00174           if (e[k+1] < 0.0)
00175             e[k] = -e[k];
00176           e.end(n-k-1) /= e[k];
00177           e[k+1] += 1.0;
00178       }
00179       e[k] = -e[k];
00180       if ((k+1 < m) & (e[k] != 0.0))
00181       {
00182         // Apply the transformation.
00183         work.end(m-k-1) = matA.corner(BottomRight,m-k-1,n-k-1) * e.end(n-k-1);
00184         for (j = k+1; j < n; ++j)
00185           matA.col(j).end(m-k-1) += (-e[j]/e[k+1]) * work.end(m-k-1);
00186       }
00187 
00188       // Place the transformation in V for subsequent back multiplication.
00189       if (wantv)
00190         m_matV.col(k).end(n-k-1) = e.end(n-k-1);
00191     }
00192   }
00193 
00194 
00195   // Set up the final bidiagonal matrix or order p.
00196   int p = (std::min)(n,m+1);
00197   if (nct < n)
00198     m_sigma[nct] = matA(nct,nct);
00199   if (m < p)
00200     m_sigma[p-1] = 0.0;
00201   if (nrt+1 < p)
00202     e[nrt] = matA(nrt,p-1);
00203   e[p-1] = 0.0;
00204 
00205   // If required, generate U.
00206   if (wantu)
00207   {
00208     for (j = nct; j < nu; ++j)
00209     {
00210       m_matU.col(j).setZero();
00211       m_matU(j,j) = 1.0;
00212     }
00213     for (k = nct-1; k >= 0; k--)
00214     {
00215       if (m_sigma[k] != 0.0)
00216       {
00217         for (j = k+1; j < nu; ++j)
00218         {
00219           Scalar t = m_matU.col(k).end(m-k).eigen2_dot(m_matU.col(j).end(m-k)); // FIXME is it really a dot product we want ?
00220           t = -t/m_matU(k,k);
00221           m_matU.col(j).end(m-k) += t * m_matU.col(k).end(m-k);
00222         }
00223         m_matU.col(k).end(m-k) = - m_matU.col(k).end(m-k);
00224         m_matU(k,k) = Scalar(1) + m_matU(k,k);
00225         if (k-1>0)
00226           m_matU.col(k).start(k-1).setZero();
00227       }
00228       else
00229       {
00230         m_matU.col(k).setZero();
00231         m_matU(k,k) = 1.0;
00232       }
00233     }
00234   }
00235 
00236   // If required, generate V.
00237   if (wantv)
00238   {
00239     for (k = n-1; k >= 0; k--)
00240     {
00241       if ((k < nrt) & (e[k] != 0.0))
00242       {
00243         for (j = k+1; j < nu; ++j)
00244         {
00245           Scalar t = m_matV.col(k).end(n-k-1).eigen2_dot(m_matV.col(j).end(n-k-1)); // FIXME is it really a dot product we want ?
00246           t = -t/m_matV(k+1,k);
00247           m_matV.col(j).end(n-k-1) += t * m_matV.col(k).end(n-k-1);
00248         }
00249       }
00250       m_matV.col(k).setZero();
00251       m_matV(k,k) = 1.0;
00252     }
00253   }
00254 
00255   // Main iteration loop for the singular values.
00256   int pp = p-1;
00257   int iter = 0;
00258   Scalar eps = ei_pow(Scalar(2),ei_is_same_type<Scalar,float>::ret ? Scalar(-23) : Scalar(-52));
00259   while (p > 0)
00260   {
00261     int k=0;
00262     int kase=0;
00263 
00264     // Here is where a test for too many iterations would go.
00265 
00266     // This section of the program inspects for
00267     // negligible elements in the s and e arrays.  On
00268     // completion the variables kase and k are set as follows.
00269 
00270     // kase = 1     if s(p) and e[k-1] are negligible and k<p
00271     // kase = 2     if s(k) is negligible and k<p
00272     // kase = 3     if e[k-1] is negligible, k<p, and
00273     //              s(k), ..., s(p) are not negligible (qr step).
00274     // kase = 4     if e(p-1) is negligible (convergence).
00275 
00276     for (k = p-2; k >= -1; --k)
00277     {
00278       if (k == -1)
00279           break;
00280       if (ei_abs(e[k]) <= eps*(ei_abs(m_sigma[k]) + ei_abs(m_sigma[k+1])))
00281       {
00282           e[k] = 0.0;
00283           break;
00284       }
00285     }
00286     if (k == p-2)
00287     {
00288       kase = 4;
00289     }
00290     else
00291     {
00292       int ks;
00293       for (ks = p-1; ks >= k; --ks)
00294       {
00295         if (ks == k)
00296           break;
00297         Scalar t = (ks != p ? ei_abs(e[ks]) : Scalar(0)) + (ks != k+1 ? ei_abs(e[ks-1]) : Scalar(0));
00298         if (ei_abs(m_sigma[ks]) <= eps*t)
00299         {
00300           m_sigma[ks] = 0.0;
00301           break;
00302         }
00303       }
00304       if (ks == k)
00305       {
00306         kase = 3;
00307       }
00308       else if (ks == p-1)
00309       {
00310         kase = 1;
00311       }
00312       else
00313       {
00314         kase = 2;
00315         k = ks;
00316       }
00317     }
00318     ++k;
00319 
00320     // Perform the task indicated by kase.
00321     switch (kase)
00322     {
00323 
00324       // Deflate negligible s(p).
00325       case 1:
00326       {
00327         Scalar f(e[p-2]);
00328         e[p-2] = 0.0;
00329         for (j = p-2; j >= k; --j)
00330         {
00331           Scalar t(internal::hypot(m_sigma[j],f));
00332           Scalar cs(m_sigma[j]/t);
00333           Scalar sn(f/t);
00334           m_sigma[j] = t;
00335           if (j != k)
00336           {
00337             f = -sn*e[j-1];
00338             e[j-1] = cs*e[j-1];
00339           }
00340           if (wantv)
00341           {
00342             for (i = 0; i < n; ++i)
00343             {
00344               t = cs*m_matV(i,j) + sn*m_matV(i,p-1);
00345               m_matV(i,p-1) = -sn*m_matV(i,j) + cs*m_matV(i,p-1);
00346               m_matV(i,j) = t;
00347             }
00348           }
00349         }
00350       }
00351       break;
00352 
00353       // Split at negligible s(k).
00354       case 2:
00355       {
00356         Scalar f(e[k-1]);
00357         e[k-1] = 0.0;
00358         for (j = k; j < p; ++j)
00359         {
00360           Scalar t(internal::hypot(m_sigma[j],f));
00361           Scalar cs( m_sigma[j]/t);
00362           Scalar sn(f/t);
00363           m_sigma[j] = t;
00364           f = -sn*e[j];
00365           e[j] = cs*e[j];
00366           if (wantu)
00367           {
00368             for (i = 0; i < m; ++i)
00369             {
00370               t = cs*m_matU(i,j) + sn*m_matU(i,k-1);
00371               m_matU(i,k-1) = -sn*m_matU(i,j) + cs*m_matU(i,k-1);
00372               m_matU(i,j) = t;
00373             }
00374           }
00375         }
00376       }
00377       break;
00378 
00379       // Perform one qr step.
00380       case 3:
00381       {
00382         // Calculate the shift.
00383         Scalar scale = (std::max)((std::max)((std::max)((std::max)(
00384                         ei_abs(m_sigma[p-1]),ei_abs(m_sigma[p-2])),ei_abs(e[p-2])),
00385                         ei_abs(m_sigma[k])),ei_abs(e[k]));
00386         Scalar sp = m_sigma[p-1]/scale;
00387         Scalar spm1 = m_sigma[p-2]/scale;
00388         Scalar epm1 = e[p-2]/scale;
00389         Scalar sk = m_sigma[k]/scale;
00390         Scalar ek = e[k]/scale;
00391         Scalar b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/Scalar(2);
00392         Scalar c = (sp*epm1)*(sp*epm1);
00393         Scalar shift = 0.0;
00394         if ((b != 0.0) || (c != 0.0))
00395         {
00396           shift = ei_sqrt(b*b + c);
00397           if (b < 0.0)
00398             shift = -shift;
00399           shift = c/(b + shift);
00400         }
00401         Scalar f = (sk + sp)*(sk - sp) + shift;
00402         Scalar g = sk*ek;
00403 
00404         // Chase zeros.
00405 
00406         for (j = k; j < p-1; ++j)
00407         {
00408           Scalar t = internal::hypot(f,g);
00409           Scalar cs = f/t;
00410           Scalar sn = g/t;
00411           if (j != k)
00412             e[j-1] = t;
00413           f = cs*m_sigma[j] + sn*e[j];
00414           e[j] = cs*e[j] - sn*m_sigma[j];
00415           g = sn*m_sigma[j+1];
00416           m_sigma[j+1] = cs*m_sigma[j+1];
00417           if (wantv)
00418           {
00419             for (i = 0; i < n; ++i)
00420             {
00421               t = cs*m_matV(i,j) + sn*m_matV(i,j+1);
00422               m_matV(i,j+1) = -sn*m_matV(i,j) + cs*m_matV(i,j+1);
00423               m_matV(i,j) = t;
00424             }
00425           }
00426           t = internal::hypot(f,g);
00427           cs = f/t;
00428           sn = g/t;
00429           m_sigma[j] = t;
00430           f = cs*e[j] + sn*m_sigma[j+1];
00431           m_sigma[j+1] = -sn*e[j] + cs*m_sigma[j+1];
00432           g = sn*e[j+1];
00433           e[j+1] = cs*e[j+1];
00434           if (wantu && (j < m-1))
00435           {
00436             for (i = 0; i < m; ++i)
00437             {
00438               t = cs*m_matU(i,j) + sn*m_matU(i,j+1);
00439               m_matU(i,j+1) = -sn*m_matU(i,j) + cs*m_matU(i,j+1);
00440               m_matU(i,j) = t;
00441             }
00442           }
00443         }
00444         e[p-2] = f;
00445         iter = iter + 1;
00446       }
00447       break;
00448 
00449       // Convergence.
00450       case 4:
00451       {
00452         // Make the singular values positive.
00453         if (m_sigma[k] <= 0.0)
00454         {
00455           m_sigma[k] = m_sigma[k] < Scalar(0) ? -m_sigma[k] : Scalar(0);
00456           if (wantv)
00457             m_matV.col(k).start(pp+1) = -m_matV.col(k).start(pp+1);
00458         }
00459 
00460         // Order the singular values.
00461         while (k < pp)
00462         {
00463           if (m_sigma[k] >= m_sigma[k+1])
00464             break;
00465           Scalar t = m_sigma[k];
00466           m_sigma[k] = m_sigma[k+1];
00467           m_sigma[k+1] = t;
00468           if (wantv && (k < n-1))
00469             m_matV.col(k).swap(m_matV.col(k+1));
00470           if (wantu && (k < m-1))
00471             m_matU.col(k).swap(m_matU.col(k+1));
00472           ++k;
00473         }
00474         iter = 0;
00475         p--;
00476       }
00477       break;
00478     } // end big switch
00479   } // end iterations
00480 }
00481 
00482 template<typename MatrixType>
00483 SVD<MatrixType>& SVD<MatrixType>::sort()
00484 {
00485   int mu = m_matU.rows();
00486   int mv = m_matV.rows();
00487   int n  = m_matU.cols();
00488 
00489   for (int i=0; i<n; ++i)
00490   {
00491     int  k = i;
00492     Scalar p = m_sigma.coeff(i);
00493 
00494     for (int j=i+1; j<n; ++j)
00495     {
00496       if (m_sigma.coeff(j) > p)
00497       {
00498         k = j;
00499         p = m_sigma.coeff(j);
00500       }
00501     }
00502     if (k != i)
00503     {
00504       m_sigma.coeffRef(k) = m_sigma.coeff(i);  // i.e.
00505       m_sigma.coeffRef(i) = p;                 // swaps the i-th and the k-th elements
00506 
00507       int j = mu;
00508       for(int s=0; j!=0; ++s, --j)
00509         std::swap(m_matU.coeffRef(s,i), m_matU.coeffRef(s,k));
00510 
00511       j = mv;
00512       for (int s=0; j!=0; ++s, --j)
00513         std::swap(m_matV.coeffRef(s,i), m_matV.coeffRef(s,k));
00514     }
00515   }
00516   return *this;
00517 }
00518 
00524 template<typename MatrixType>
00525 template<typename OtherDerived, typename ResultType>
00526 bool SVD<MatrixType>::solve(const MatrixBase<OtherDerived> &b, ResultType* result) const
00527 {
00528   const int rows = m_matU.rows();
00529   ei_assert(b.rows() == rows);
00530 
00531   Scalar maxVal = m_sigma.cwise().abs().maxCoeff();
00532   for (int j=0; j<b.cols(); ++j)
00533   {
00534     Matrix<Scalar,MatrixUType::RowsAtCompileTime,1> aux = m_matU.transpose() * b.col(j);
00535 
00536     for (int i = 0; i <m_matU.cols(); ++i)
00537     {
00538       Scalar si = m_sigma.coeff(i);
00539       if (ei_isMuchSmallerThan(ei_abs(si),maxVal))
00540         aux.coeffRef(i) = 0;
00541       else
00542         aux.coeffRef(i) /= si;
00543     }
00544 
00545     result->col(j) = m_matV * aux;
00546   }
00547   return true;
00548 }
00549 
00558 template<typename MatrixType>
00559 template<typename UnitaryType, typename PositiveType>
00560 void SVD<MatrixType>::computeUnitaryPositive(UnitaryType *unitary,
00561                                              PositiveType *positive) const
00562 {
00563   ei_assert(m_matU.cols() == m_matV.cols() && "Polar decomposition is only for square matrices");
00564   if(unitary) *unitary = m_matU * m_matV.adjoint();
00565   if(positive) *positive = m_matV * m_sigma.asDiagonal() * m_matV.adjoint();
00566 }
00567 
00576 template<typename MatrixType>
00577 template<typename UnitaryType, typename PositiveType>
00578 void SVD<MatrixType>::computePositiveUnitary(UnitaryType *positive,
00579                                              PositiveType *unitary) const
00580 {
00581   ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices");
00582   if(unitary) *unitary = m_matU * m_matV.adjoint();
00583   if(positive) *positive = m_matU * m_sigma.asDiagonal() * m_matU.adjoint();
00584 }
00585 
00595 template<typename MatrixType>
00596 template<typename RotationType, typename ScalingType>
00597 void SVD<MatrixType>::computeRotationScaling(RotationType *rotation, ScalingType *scaling) const
00598 {
00599   ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices");
00600   Scalar x = (m_matU * m_matV.adjoint()).determinant(); // so x has absolute value 1
00601   Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> sv(m_sigma);
00602   sv.coeffRef(0) *= x;
00603   if(scaling) scaling->lazyAssign(m_matV * sv.asDiagonal() * m_matV.adjoint());
00604   if(rotation)
00605   {
00606     MatrixType m(m_matU);
00607     m.col(0) /= x;
00608     rotation->lazyAssign(m * m_matV.adjoint());
00609   }
00610 }
00611 
00621 template<typename MatrixType>
00622 template<typename ScalingType, typename RotationType>
00623 void SVD<MatrixType>::computeScalingRotation(ScalingType *scaling, RotationType *rotation) const
00624 {
00625   ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices");
00626   Scalar x = (m_matU * m_matV.adjoint()).determinant(); // so x has absolute value 1
00627   Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> sv(m_sigma);
00628   sv.coeffRef(0) *= x;
00629   if(scaling) scaling->lazyAssign(m_matU * sv.asDiagonal() * m_matU.adjoint());
00630   if(rotation)
00631   {
00632     MatrixType m(m_matU);
00633     m.col(0) /= x;
00634     rotation->lazyAssign(m * m_matV.adjoint());
00635   }
00636 }
00637 
00638 
00642 template<typename Derived>
00643 inline SVD<typename MatrixBase<Derived>::PlainObject>
00644 MatrixBase<Derived>::svd() const
00645 {
00646   return SVD<PlainObject>(derived());
00647 }
00648 
00649 #endif // EIGEN2_SVD_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:44