ExampleSVD.cpp

/* This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

/** @example ExampleSVD.cpp
    @relates SVD3x3, SVD
    @ingroup g_examples    
    @brief   Example for singular value decomposition
    @author  woelk 07/2004
*/

#include <bias_config.h>
#include <MathAlgo/SVD.hh>
#include <MathAlgo/SVD3x3.hh>
#include <Base/Math/Random.hh>
#include <Base/Math/Matrix3x3.hh>
#include <Base/Debug/TimeMeasure.hh>

using namespace BIAS;
using namespace std;

#define MAX_PRINT_SIZE 8

int main(int argc, char **argv)
{
  int size = 3;
  if (argc > 1) {
    size = atoi(argv[1]);
  }
  bool useOpenCV = false;
  if (argc > 2)
    useOpenCV = (atoi(argv[2]) != 0);
  if (size < 1) {
    cout << "Usage : ExampleSVD [<matrix size>] [<use OpenCV?>]" << endl << endl;
    return 0;
  }

  if (useOpenCV) {
#ifdef BIAS_HAVE_OPENCV
    cout << endl << "Using OpenCV singular value decomposition" << endl;
#else
    cout << "BIAS does not support OpenCV, using LAPACK SVD instead" << endl;
    useOpenCV = false;
#endif
  } else {
    cout << endl << "Using LAPACK singular value decomposition" << endl;
  }
  double errorThreshold = 1e-6 * size * size;

  // 1.) Compute SVD for random NxN matrix

  cout << endl << "1.) Compute SVD for random " << size << "x" << size << " matrix" << endl;

  Random rand;
  Matrix<double> M(size, size, MatrixZero);
  SVD svd;
  SVD3x3 svd3x3;
  TimeMeasure timer, timer3x3;
 
  for (int i = 0; i < size; i++)
    for (int j = 0; j < size; j++)
      M[i][j] = rand.GetUniformDistributed(-10.0, 10.0);
  if (size <= MAX_PRINT_SIZE)
    M.PrintPretty(cout, "M");

  int numRepeats = 100;
  if (size >= 200)
    numRepeats = 1;
  else if (size >= 100)
    numRepeats = 10;
  if (size == 3) {
    timer3x3.Start();
    for (int k = 0; k < numRepeats; k++)
      svd3x3.Compute(M);
    timer3x3.Stop();
  }
  timer.Start();
  for (int k = 0; k < numRepeats; k++)
#ifdef BIAS_HAVE_OPENCV
    if (useOpenCV)
      svd.ComputeOpenCV(M);
    else
      svd.Compute(M);
#else
    svd.Compute(M);
#endif
  timer.Stop();

  if (size <= MAX_PRINT_SIZE)
    svd.GetU().PrintPretty(cout, "SVD    U");
  if (size == 3) {
    svd3x3.GetU().PrintPretty(cout, "SVD3x3 U");
    Matrix<double> dU(svd.GetU() - svd3x3.GetU());
    dU.PrintPretty(cout, "diff.  U ");
  }

  if (size <= MAX_PRINT_SIZE)
    cout << endl << "SVD    S  = " << svd.GetS() << endl;
  if (size == 3) {
    Vector<double> dS(svd.GetS() - svd3x3.GetS());
    cout << endl << "SVD3x3 S  = " << svd3x3.GetS() << endl
         << endl << "diff.  S  = " << dS << endl;
  }

  if (size <= MAX_PRINT_SIZE)
    svd.GetVT().PrintPretty(cout, "SVD    VT");
  if (size == 3) {
    svd3x3.GetVT().PrintPretty(cout, "SVD3x3 VT");
    Matrix<double> dVT(svd.GetVT() - svd3x3.GetVT());
    dVT.PrintPretty(cout, "diff.  VT");
  }

  Matrix<double> S(size, size, MatrixZero);
  for (int i = 0; i < size; i++)
    S[i][i] = svd.GetS()[i];
  Matrix<double> dM = M - svd.GetU() * S * svd.GetVT();
  cout << endl << "SVD    |M - U*S*VT| = " << dM.NormL2() << endl;
  if (dM.NormL2() > errorThreshold) {
    BIASERR("SVD failed!");
    dM.PrintPretty(cout, "M - U*S*VT");
  }
  if (size == 3) {
    for (int i = 0; i < size; i++)
      S[i][i] = svd3x3.GetS()[i];
    dM = M - svd3x3.GetU() * S * svd3x3.GetVT();
    cout << endl << "SVD3x3 |M - U*S*VT| = " << dM.NormL2() << endl;
    if (dM.NormL2() > errorThreshold) {
      BIASERR("SVD3x3 failed!");
      if (size <= MAX_PRINT_SIZE)
        dM.PrintPretty(cout, "M - U*S*VT");
    }
  }

  cout << endl << "SVD    time = " << (0.001 * timer.GetRealTime()) / numRepeats << " ms";
  if (size == 3) {
    cout << endl << "SVD3x3 time = " << (0.001 * timer3x3.GetRealTime()) / numRepeats << " ms";
  }

  // 2.) Invert matrix with SVD

  cout << endl << endl << "2.) Invert matrix with SVD" << endl;

  Matrix<double> I(size, size, MatrixIdentity);
  Matrix<double> invM;
#ifdef BIAS_HAVE_OPENCV
    if (useOpenCV)
      invM = svd.InvertOpenCV(M);
    else
      invM = svd.Invert(M);
#else
    invM = svd.Invert(M);
#endif
  dM = M * invM - I;
  cout << endl << "SVD    |M*inv(M) - I| = " << dM.NormL2() << endl;
  if (dM.NormL2() > errorThreshold) {
    BIASERR("Matrix inversion with SVD failed!");
    if (size <= MAX_PRINT_SIZE)
      dM.PrintPretty(cout, "M*inv(M) - I");
  }
  if (size == 3) {
    invM = svd3x3.Invert();
    dM = M * invM - I;
    cout << endl << "SVD3x3 |M*inv(M) - I| = " << dM.NormL2() << endl;
    if (dM.NormL2() > errorThreshold) {
      BIASERR("Matrix inversion with SVD3x3 failed!");
      dM.PrintPretty(cout, "M*inv(M) - I");
    }
  }

  // 3.) Compute matrix square root with SVD

  cout << endl << endl << "3.) Compute matrix square root with SVD" << endl;

  Matrix<double> MTM = M.Transpose() * M;
  if (size <= MAX_PRINT_SIZE)
    MTM.PrintPretty(cout, "MT*M");

  Matrix<double> sqrtMTM = svd.Sqrt(MTM);
  if (size <= MAX_PRINT_SIZE)
    sqrtMTM.PrintPretty(cout, "SVD    sqrt(MT*M)");
  dM = MTM - (sqrtMTM * sqrtMTM);
  cout << endl << "SVD    |MT*M - sqrt(MT*M)^2| = " << dM.NormL2() << endl;
  if (dM.NormL2() > errorThreshold) {
    BIASERR("Matrix square root failed!");
    if (size <= MAX_PRINT_SIZE)
      dM.PrintPretty(cout, "MT*M - sqrt(MT*M)^2");
  }
  if (size == 3) {
    sqrtMTM = svd3x3.Sqrt(MTM);
    sqrtMTM.PrintPretty(cout, "SVD3x3 sqrt(MT*M)");
    dM = MTM - (sqrtMTM * sqrtMTM);
    cout << endl << "SVD3x3 |MT*M - sqrt(MT*M)^2| = " << dM.NormL2() << endl;
    if (dM.NormL2() > errorThreshold) {
      BIASERR("Matrix square root failed!");
      dM.PrintPretty(cout, "MT*M - sqrt(MT*M)^2");
    }
  }
  cout << endl;

  return 0;
}