Basic Image AlgorithmS Library 2.8.0

SVD.cpp

00001 /* 
00002 This file is part of the BIAS library (Basic ImageAlgorithmS).
00003 
00004 Copyright (C) 2003-2009    (see file CONTACT for details)
00005   Multimediale Systeme der Informationsverarbeitung
00006   Institut fuer Informatik
00007   Christian-Albrechts-Universitaet Kiel
00008 
00009 
00010 BIAS is free software; you can redistribute it and/or modify
00011 it under the terms of the GNU Lesser General Public License as published by
00012 the Free Software Foundation; either version 2.1 of the License, or
00013 (at your option) any later version.
00014 
00015 BIAS is distributed in the hope that it will be useful,
00016 but WITHOUT ANY WARRANTY; without even the implied warranty of
00017 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018 GNU Lesser General Public License for more details.
00019 
00020 You should have received a copy of the GNU Lesser General Public License
00021 along with BIAS; if not, write to the Free Software
00022 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00023 */
00024 
00025 
00026 #include "SVD.hh"
00027 
00028 #include <Base/Debug/Exception.hh>
00029 #include <Base/Common/CompareFloatingPoint.hh>
00030 
00031 #include <algorithm>
00032 #include <cmath>
00033 #include <vector>
00034 
00035 #include "Lapack.hh"
00036 #ifdef BIAS_HAVE_OPENCV
00037 #include <cv.h>
00038 #include <cxcore.h>
00039 #endif
00040 using namespace BIAS;
00041 using namespace std;
00042 
00043 SVD::SVD(const Matrix<double> &M, double ZeroThreshold,
00044          bool AbsoluteZeroThreshold) {
00045   Decomposed_=false;
00046   AbsoluteZeroThreshold_ =  AbsoluteZeroThreshold;
00047   Compute( M, ZeroThreshold);
00048 }
00049 
00050 #ifdef BIAS_HAVE_OPENCV
00051 int SVD::ComputeOpenCV(const Matrix<double>& M,  double ZeroThreshold){
00052   ZeroThreshold_ = ZeroThreshold;
00053   A_ = M;
00054   unsigned cols = M.GetCols();
00055   unsigned rows = M.GetRows();
00056   CvMat* mat = cvCreateMat(rows,cols,CV_32FC1);
00057   cvZero(mat);
00058   //copy BIAS to OpenCV matrix
00059   for(unsigned i=0;i<rows;i++){
00060     for(unsigned j=0;j<cols;j++){
00061       CV_MAT_ELEM( *mat, float, i, j ) = (float)M[i][j];
00062     }
00063   }
00064 
00065   CvMat* U  = cvCreateMat(rows,rows,CV_32FC1); cvZero(U);
00066   CvMat* D  = cvCreateMat(rows,cols,CV_32FC1); cvZero(D);
00067   CvMat* V  = cvCreateMat(cols,cols,CV_32FC1); cvZero(V);
00068   cvSVD(mat, D, U, V, CV_SVD_V_T); // A = U D V^T
00069 
00070   U_.newsize(rows,rows);
00071   VT_.newsize(cols,cols);
00072   S_.newsize(rows);
00073   for(unsigned i=0;i<rows;i++){
00074     for(unsigned j=0;j<rows;j++){
00075       U_[i][j] = (double)CV_MAT_ELEM(*U,float,i,j);    
00076     }   
00077     if(i< rows && i < cols){
00078       //singular values are on diagonal of matrix D
00079       S_[i] = (double)CV_MAT_ELEM(*D,float,i,i);
00080     }
00081   }
00082   for(unsigned i=0;i<cols;i++){
00083     for(unsigned j=0;j<cols;j++){
00084       VT_[i][j] = (double)CV_MAT_ELEM(*V,float,i,j);
00085     }
00086   }
00087    
00088   cvReleaseMat(&U);
00089   cvReleaseMat(&D);
00090   cvReleaseMat(&V);
00091   cvReleaseMat(&mat);
00092   
00093   //need extra if -> ensure S_.size != 0
00094   if (S_[0] < numeric_limits<double>::epsilon())
00095     return -1;
00096   else 
00097     Decomposed_=true;
00098 
00099   return 0;
00100 }
00101 #endif
00102 int SVD::Compute(const Matrix<double> & M, double ZeroThreshold) 
00103 {
00104   Decomposed_=false;
00105   ZeroThreshold_ = ZeroThreshold;
00106   return General_Eigenproblem_GeneralMatrix_Lapack( M );
00107 }
00108 
00109 
00110 int SVD::General_Eigenproblem_GeneralMatrix_Lapack(const Matrix<double> & M) 
00111 {
00112   // S_ will contain the singular values of A sorted with s(i) <= s(i+1)  
00113   // U_ will contain m x m left orthogonal/unitary matrix with 
00114   //    eigenvectors in columns  
00115   // VT_ will contain n x n right _transposed_ orthogonal/unitary matrix 
00116   //     with eigenvectors in rows (because of transpose)    
00117   A_=M;
00118   int res=0;
00119   res = General_singular_value_decomposition(A_, S_, U_, VT_);
00120 
00121   if ( S_.size() != 0 && U_.num_rows() != 0 &&  U_.num_cols() != 0 &&
00122        VT_.num_rows() != 0 &&  VT_.num_cols() != 0 && res == 0) {
00123     //need extra if -> ensure S_.size != 0
00124     if (S_[0] < numeric_limits<double>::epsilon()) {
00125       res = -1;
00126     } else Decomposed_=true;
00127   } else {
00128     res = -1;
00129   }
00130   
00131   return res;
00132 }
00133 
00134 // -- Solve the matrix-vector system M x = y, returning x.
00135 Vector<double> SVD::Solve(const Vector<double>& y) const
00136 {
00137   /// @bug bug in this function (woelk)
00138   BEXCEPTION("Error, there is a bug in this function, use Solve(Matrix, "
00139              "vector, Vector instead");
00140   Vector<double> x(VT_.num_cols());      // solution matrix vector
00141   Vector<double> UTy;
00142 
00143   //cout << "U_=" << U_ << " VT_= " << VT_ << " S_=" << S_;
00144 
00145   if(Decomposed_) {
00146     if (y.size() <= U_.num_rows()){
00147       Matrix<double> UT=U_.Transpose();
00148       // at first multiply with U^T
00149       if (U_.num_rows() < U_.num_cols()) {// augment y with extra rows of
00150         Vector<double> yy(U_.num_rows(), 0.0);      // zeros, so that it match
00151         for (int k=0; k< y.size();k++) 
00152           yy[k] = y[k]; // cols of u.transpose.
00153         UTy = UT * yy;
00154       } else {
00155         UTy = UT * y;
00156       }
00157            
00158       // multiply with diagonal 1/S_
00159       int i=0;
00160       for (; i < UTy.size(); i++)      
00161         x[i]=(fabs(S_[i]) >= ZeroThreshold_) ? UTy[i]/S_[i] : 0;//UTy[i]*S_[i];
00162       
00163       for (;i<x.size();i++)  x[i]=0;
00164       
00165       return VT_.Transpose() * x;   // premultiply with v.
00166       
00167     } else {
00168       BIASERR("size of right hand side of Ax = y is too big. Length(y) = "
00169               << y.size());
00170       x.newsize(0);
00171       return x;
00172     }
00173   } else {
00174     BIASERR("SVD holds no decomposed matrix.");
00175     x.newsize(0);
00176     return x;
00177   }
00178 }
00179 
00180 #ifdef BIAS_HAVE_OPENCV
00181     /** \brief Calculates new inverse with OpenCV cvInvert*/
00182 Matrix<double> SVD::InvertOpenCV(const Matrix<double>& A)
00183 {
00184   unsigned cols = A.GetCols();
00185   unsigned rows = A.GetRows();
00186   Matrix<double> B(rows,cols,0);
00187   CvMat* a = cvCreateMat(rows,cols,CV_32FC1);
00188   cvZero(a);
00189   //copy BIAS to OpenCV matrix
00190   for(unsigned i=0;i<rows;i++){
00191     for(unsigned j=0;j<cols;j++){
00192       CV_MAT_ELEM( *a, float, i, j ) = (float)A[i][j];
00193     }
00194   }
00195   CvMat* b = cvCreateMat(rows,cols,CV_32FC1);
00196   cvZero(b);
00197   //CV_LU - Gaussian elimination with optimal pivot element 
00198   //CV_SVD - Singular decomposition method 
00199   double ret = cvInvert(a,b,CV_SVD);
00200   if(ret==0.0) BIASWARN("Could not invert Matrix");
00201   //copy OpenCV to BIAS matrix
00202   for(unsigned i=0;i<rows;i++){
00203     for(unsigned j=0;j<cols;j++){
00204       B[i][j] = (double)CV_MAT_ELEM( *a, float, i, j );
00205     }
00206   }
00207   cvReleaseMat(&a);
00208   cvReleaseMat(&b); 
00209 
00210   return B;
00211 }
00212 #endif
00213 
00214 Matrix<double> SVD::Invert()
00215 {
00216   if (!Decomposed_) {//no decomposition, return 0x0-matrix
00217     Matrix<double> res(0,0);
00218     return res;
00219   }
00220 
00221   Matrix<double> V(VT_.num_cols(), VT_.num_rows());
00222   Matrix<double> UT(U_.num_cols(), U_.num_rows());
00223   Matrix<double> SdT(VT_.num_rows(), U_.num_cols());
00224   SdT.SetZero();
00225 
00226   if (AbsoluteZeroThreshold_) {
00227     for (register int i = 0; i < S_.size(); i++){
00228       SdT[i][i] = (S_[i] < ZeroThreshold_ && S_[i] > -ZeroThreshold_)
00229                   ? (0.0) : (1.0/S_[i]);
00230     }
00231   } else {
00232     for (register int i = 0; i < S_.size(); i++){
00233       SdT[i][i]=(S_[i]/S_[0] < ZeroThreshold_ && S_[i]/S_[0] > -ZeroThreshold_)
00234                  ? (0.0) : (1.0/S_[i]);
00235     }
00236   }
00237   V = VT_.Transpose();
00238   UT = U_.Transpose();
00239 
00240   return V * SdT * UT;
00241 }
00242 
00243 Matrix<double> SVD::Invert(int rank)
00244 {
00245   if (!Decomposed_) {//no decomposition, return 0x0-matrix
00246     Matrix<double> res(0,0);
00247     return res;
00248   }
00249 
00250   Matrix<double> V(VT_.num_cols(), VT_.num_rows());
00251   Matrix<double> UT(U_.num_cols(), U_.num_rows());
00252   Matrix<double> SdT(VT_.num_rows(), U_.num_cols());
00253   SdT.SetZero();
00254 
00255   for (register int i = 0; i < rank; i++){
00256     SdT[i][i] = (1.0/S_[i]);
00257   }
00258   for (register int i = rank; i < S_.size(); i++){
00259     SdT[i][i] = 0;
00260   }
00261   V = VT_.Transpose();
00262   UT = U_.Transpose();
00263 
00264   return V * SdT * UT;
00265 }
00266 
00267 
00268 Matrix<double> SVD::Invert(Matrix<double> A)
00269 {
00270   if (Compute(A)==0)
00271     return Invert();
00272   else return Matrix<double>(0,0);
00273 } 
00274 
00275 
00276 Matrix<double> SVD::Sqrt()
00277 {
00278 #ifdef BIAS_DEBUG
00279   // check symmetry and symmetric positive definit'ness
00280   // 1. symmetry
00281   bool valid = (A_.GetCols() == A_.GetRows());
00282   if (valid){
00283     const int num = A_.GetCols();
00284     for (int r=0; r<num; r++){
00285       for (int c=r; c<num; c++){
00286         // dont use double epsilon but 1e-10 for numeric comparison
00287         valid = valid && Equal(A_[r][c], A_[c][r], 1e-10);
00288         if (!valid) {
00289           cout << setprecision(20) << A_[r][c] << "\t"<< setprecision(20) << A_[c][r]<<endl;
00290           break;
00291           break;
00292         }
00293       }
00294     }
00295   }
00296   if (!valid){
00297     //BIASERR("matrix is not symmetric "<<A_);
00298     BEXCEPTION("matrix is not symmetric "/*<<A_*/);
00299   }
00300   // 2. positive definit = all eigenvalues are positive
00301   Vector<double> eig = Upper_symmetric_eigenvalue_solve(A_);
00302   for (register int i=0; i<eig.size(); i++){
00303     valid = valid && (GreaterEqual(eig[i], 0.0, 1e-10));
00304     if (!valid){
00305       cout << "eigenvalue "<<setprecision(20)<<eig[i]<<endl;
00306       break;
00307     }
00308   }
00309   if (!valid){
00310     //BIASERR("matrix is not positive definit "<<A_);
00311     BEXCEPTION("matrix is not positive definit "/*<<A_*/);
00312   }
00313 #endif
00314   if (!Decomposed_){
00315     int res = 0;
00316     if ((res = Compute(A_))!=0) {
00317       BIASERR("SVD failed with result "<<res);
00318       return  Matrix<double>(VT_.num_rows(), U_.num_cols(), MatrixZero);
00319     }
00320   }
00321   // [U S V] = svd(A); sqrtm(A) = U * sqrt(S) * V'
00322   Matrix<double> SqrtS(VT_.num_rows(), U_.num_cols());
00323   SqrtS.SetZero();
00324   if (AbsoluteZeroThreshold_) {
00325     for (register int i=0; i<S_.size(); i++){
00326       SqrtS[i][i] = (fabs(S_[i])<ZeroThreshold_)?(0.0):(sqrt(S_[i]));
00327     }
00328   } else {
00329     for (register int i=0; i<S_.size(); i++){
00330       SqrtS[i][i] = (fabs(S_[i]/S_[0])<ZeroThreshold_)?(0.0):(sqrt(S_[i]));
00331     }
00332   }
00333 
00334   return U_ * SqrtS * VT_;
00335 }
00336 
00337 Matrix<double> SVD::Sqrt(const Matrix<double>& A)
00338 {
00339   if (Compute(A)==0)
00340     return Sqrt();
00341   else return Matrix<double>(0,0);
00342 } 
00343 
00344 Matrix<double> SVD::SqrtT()
00345 {
00346 #ifdef BIAS_DEBUG
00347   // check symmetry and symmetric positive definit'ness
00348   // 1. symmetry
00349   bool valid = (A_.GetCols() == A_.GetRows());
00350   if (valid){
00351     const int num = A_.GetCols();
00352     for (int r=0; r<num; r++){
00353       for (int c=r; c<num; c++){
00354         valid = valid && Equal(A_[r][c], A_[c][r], 1e-12);
00355         if (!valid) {
00356           cout<<" Invalid matrix ("<<num<<"x"<<num<<") element at "
00357               <<r<<" "<<c<<endl;
00358           cout << setprecision(20) << A_[r][c] << "\t"<< setprecision(20) 
00359                << A_[c][r]<<endl;
00360           break;
00361         }
00362       }
00363       if (!valid) { break; }
00364     }
00365   }
00366   if (!valid){
00367     BEXCEPTION("matrix is not symmetric "<<A_<<endl
00368                <<"Fix this by calling A_.MakeSymmetric() in CALLING class");
00369   }
00370 
00371   // 2. positive definit = all eigenvalues are positive
00372   Vector<double> eig = Upper_symmetric_eigenvalue_solve(A_);
00373   for (register int i=0; i<eig.size(); i++){
00374     valid = valid && (GreaterEqual(eig[i], 0.0, 1e-12));
00375     if (!valid){
00376       //cout << "eigenvalue "<<setprecision(20)<<eig[i]<<endl;
00377       break;
00378     }
00379   }
00380   if (!valid){
00381     //BIASERR("matrix is not positive definit "<<A_);
00382     //BEXCEPTION("matrix is not positive definit "/*<<A_*/);
00383   }
00384 #endif
00385   
00386   if (!Decomposed_){
00387     if (Compute(A_)!=0) {
00388       // error in decomposition
00389       return Matrix<double>(0,0); // empty matrix
00390     }
00391   }
00392   // [U S V] = svd(A); sqrtm(A) = U * sqrt(S)
00393   Matrix<double> SqrtS(VT_.num_rows(), U_.num_cols(), MatrixZero);  
00394   if (AbsoluteZeroThreshold_) {
00395     for (register int i=0; i<S_.size(); i++){
00396       SqrtS[i][i] = (fabs(S_[i])<ZeroThreshold_)?(0.0):(sqrt(S_[i]));
00397     }
00398   } else {
00399     for (register int i=0; i<S_.size(); i++){
00400       SqrtS[i][i] = (fabs(S_[i]/S_[0])<ZeroThreshold_)?(0.0):(sqrt(S_[i]));
00401     }
00402   }
00403 
00404   return U_ * SqrtS;
00405 }
00406 
00407 Matrix<double> SVD::SqrtT(const Matrix<double>& A)
00408 {
00409   if (Compute(A)==0)
00410     return SqrtT();
00411   else return Matrix<double>(0,0);
00412 } 
00413 
00414 
00415 const int SVD::AbsNullspaceDim() const 
00416 {
00417   int dim = 0;
00418   int index = S_.size()-1; 
00419 
00420   //cerr << "ZeroThreshold_: "<<ZeroThreshold_<<endl;
00421   while (index >= 0 && fabs(S_[index]) < ZeroThreshold_ ) {
00422     // singular values are sorted in descending order is still (around zero
00423     index--;
00424     dim++;
00425   }
00426   
00427   if (A_.num_cols()>A_.num_rows()){
00428       dim+=A_.num_cols()-A_.num_rows();
00429   }
00430   return dim;
00431 }
00432 
00433 const int SVD::RelNullspaceDim() const 
00434 { 
00435   int dim = 0;
00436   int index = S_.size()-1; 
00437 
00438   // use greatest singular value as reference
00439   const double RefSingValue(S_[0]);
00440   if (RefSingValue<=DBL_EPSILON) {
00441     // the greatest singular value is zero => all are zero
00442     return 0;
00443   } 
00444     
00445   while ((S_[index]/RefSingValue)<ZeroThreshold_){
00446     dim++;
00447     if (--index<1) break; // S[0]/S[0]==1 always => dont need to check that
00448   }
00449 
00450   // check if we have an "underdetermined equation system" by dimensionality
00451   if (A_.num_cols()>A_.num_rows()){
00452     dim += A_.num_cols()-A_.num_rows();
00453   }
00454    
00455   return dim;
00456 }
00457 
00458 
00459 const int SVD::AbsRightNullspaceDim() const { 
00460   return AbsNullspaceDim();
00461 }
00462 
00463 const int SVD::RelRightNullspaceDim() const { 
00464     return RelNullspaceDim();
00465 }
00466 
00467 const int SVD::AbsLeftNullspaceDim() const 
00468 { 
00469   int dim = 0;
00470   int index = S_.size()-1; 
00471 
00472   //cerr << "ZeroThreshold_: "<<ZeroThreshold_<<endl;
00473   while (index >= 0 && fabs(S_[index]) < ZeroThreshold_ ) {
00474     // singular values are sorted in descending order
00475     index--;
00476     dim++;
00477   }
00478   // compensate null space for rectangular matrices 
00479   if (A_.num_rows()>A_.num_cols()){
00480     dim += A_.num_rows()-A_.num_cols();
00481   }
00482   return dim;
00483 }
00484 
00485 const int SVD::RelLeftNullspaceDim() const 
00486 { 
00487   int dim = 0;
00488   int index = S_.size()-1; 
00489   
00490   if (GetSingularValue(0)==0){
00491     dim=NullspaceDim();
00492   } else {
00493     double maxs=S_[0];
00494     while (index>=0 && fabs(S_[index]/maxs)<ZeroThreshold_){
00495         index--;
00496         dim++;
00497     }
00498     // compensate null space for rectangular matrices 
00499     if (A_.num_rows()>A_.num_cols()){
00500         dim+=A_.num_rows()-A_.num_cols();
00501     }
00502   } 
00503   return dim;
00504 }
00505 
00506 unsigned int SVD::Rank()
00507 {
00508   BIASDOUT(D_SVD_RANK, "S_.Size(): "<<S_.Size()<<" NullspaceDim(): "
00509            << NullspaceDim()<<"  A_.num_cols(): "<<A_.num_cols()
00510            <<"   A_.num_rows(): "<<A_.num_rows()<<"  S_: "<<S_);
00511   unsigned int lowdim=A_.num_rows()<A_.num_cols()?A_.num_rows():A_.num_cols();
00512   return (lowdim-NullspaceDim());
00513 }
00514 
00515 
00516 int SVD::Solve(Matrix<double>& A, Vector<double>& b, Vector<double>& x)
00517 {
00518   int res = 0;
00519   // make squared symmetric matrix A^T * A
00520   Matrix<double> ATA(A.num_cols(), A.num_cols());
00521   //ATA = A.Transpose()*A;
00522   A.GetSystemMatrix(ATA);
00523 
00524   // compute ATA = U * S * V^T
00525   Compute(ATA);
00526   Vector<double> d(A.num_rows()), z(A.num_cols());
00527   Vector<double> S = GetS();
00528   BIASCDOUT(D_SVD_SOLVE, "S: "<<S<<endl); 
00529   
00530   // A * x = b  => x = (A^T * A )^-1 * A^T * b
00531   //                 = (U * S * V^T)^-1 * A^T * b 
00532   // by symmetry we know: A^T * A = (A^T * A)^T 
00533   //                 = (V * S * U^T)^-1 * A^T * b 
00534   // since U and V^T are square orthogonal matrices here
00535   //                 = U^-T * S^-1 * V^T * A^T * b
00536   //                 = U    * S^-1 *       d
00537   d = GetVT() * A.Transpose() * b;
00538   
00539   // compute z = S^-1 * d
00540   BIASCDOUT(D_SVD_SOLVE, "d: "<<d<<endl);
00541   if (AbsoluteZeroThreshold_) {
00542     for (int i=0; i<A.num_cols(); i++) {
00543       if (fabs(S[i])>ZeroThreshold_) {
00544         z[i] = d[i]/S[i];
00545       } else {
00546         z[i] = 0.0;
00547         res--;
00548       }
00549     }
00550   } else {
00551     for (int i=0; i<A.num_cols(); i++) {
00552       if (fabs(S[i]/S_[0])>ZeroThreshold_) {
00553         z[i] = d[i]/S[i];
00554       } else {
00555         z[i] = 0.0;
00556         res--;
00557       }
00558     }
00559   }
00560 
00561   if (res<0) {
00562       BIASDOUT(D_SVD_SOLVE,"Error in SVD, degenerate case in Solve ("
00563             << -res <<" dim. solution space), S=" 
00564             << S);
00565   }
00566   
00567   // get final solution x = U * z
00568   x = GetU() * z;
00569   return res;
00570 }
00571 
00572 
00573 BIAS::Vector<double> SVD::GetEigenValues()
00574 {
00575   std::vector<double> eigValues(VT_.GetRows());
00576   // now we need the eigenvalues (with sign !)
00577   // compute this by projecting the (normalized) eigenvectors onto their images
00578   Vector<double> eigenvector, eigenimage;
00579   for (unsigned int i=0; i<VT_.GetRows(); i++) {
00580     // get ith eigenvector and its image
00581     eigenvector=VT_.GetRow(i);
00582     eigenimage = A_ * eigenvector; 
00583     // compute scalar product
00584     eigValues[i] = eigenimage * eigenvector;
00585   }
00586   std::sort(eigValues.begin(),eigValues.end());
00587   BIAS::Vector<double> result(VT_.GetRows(),&eigValues[0]);
00588   
00589   return result;
00590 }
00591 
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends