ExampleLeastSquares.cpp

example for LeastSquares solver

Author

MIP

/* 
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 ExampleLeastSquares.cpp
    @relates LeastSquares
    @ingroup g_examples    
    @brief  example for LeastSquares solver
    @author MIP
*/

#include <MathAlgo/LeastSquares.hh>
#include <Base/Math/Random.hh>
#include <Base/Debug/TimeMeasure.hh>

using namespace BIAS;
using namespace std;

double sigma=5.0;

void GetSystem(Matrix<double>& A, Vector<double>& b, Vector<double>& weight,
               Vector<double>& sol, unsigned solsize, unsigned numeq)
{
  Random ran;
  sol.newsize(solsize);
  for (unsigned i=0; i<solsize; i++){
    sol[i]=ran.GetUniformDistributed(-10.0, 10.0);
  }
  
  A.newsize(numeq, solsize); 
  for (unsigned i=0; i<solsize*numeq; i++){
    A[i/solsize][i%solsize]=ran.GetUniformDistributed(-100.0, 100.0);
  }

  b.newsize(numeq);
  b=A*sol;
  for (unsigned i=0; i<numeq; i++){
    b[i]+=ran.GetNormalDistributed(0.0, sigma);
  }

  weight.newsize(numeq);  
  for (unsigned i=0; i<numeq; i++){
    weight[i]=ran.GetNormalDistributed(1.0, 0.5);
  }
}

void GetSystem(Matrix<double>& A, Vector<double>& weight,
               Vector<double>& sol, unsigned solsize, unsigned numeq)
{
  Random ran;
  double sum=0.0;
  sol.newsize(solsize);
  for (unsigned i=0; i<solsize; i++){
    sol[i]=ran.GetUniformDistributed(-10.0, 10.0);
    sum+=sol[i]*sol[i];
  }
  sum=sqrt(sum);
  for (unsigned i=0; i<solsize; i++){
    sol[i]/=sum;
  }
  
  A.newsize(numeq, solsize); 
  for (unsigned r=0; r<numeq; r++){
    sum=0.0;
    for (unsigned c=0; c<solsize-1; c++){
      A[r][c]=ran.GetUniformDistributed(-100.0, 100.0);
      sum+=A[r][c]*sol[c];
    }
    A[r][solsize-1]=-sum/sol[solsize-1];
  }

  //cerr << "Ax: "<<A*sol<<endl;

  for (unsigned r=0; r<numeq; r++){
    for (unsigned c=0; c<solsize-1; c++){
      A[r][c]+=ran.GetNormalDistributed(0.0, sigma);
    }
  }

  //cerr << "Ax: "<<A*sol<<endl;

  weight.newsize(numeq);  
  for (unsigned i=0; i<numeq; i++){
    weight[i]=ran.GetNormalDistributed(1.0, 0.5);
  }
}

void TestNoWeightAxb(Matrix<double>& A, Vector<double>& b, 
                     Vector<double>& sol, unsigned solsize,
                     LeastSquaresLapack& lsla, LeastSquaresSVD lssvd)
{
  Vector<double> x;

  cout << "solution:   ";
  for (unsigned i=0; i<solsize; i++) { cout <<setw(15)<<sol[i] << "    "; }
  cout << endl;
  
  if (lsla.Solve(A, b, x)==0){
    cout << "lapack sol: ";
    for (unsigned i=0; i<solsize; i++) { cout<<setw(15)<<x[i]<<"    "; } 
    cout <<endl;
  }
  if (lssvd.Solve(A, b, x)==0){
    cout << "svd sol:    ";
    for (unsigned i=0; i<solsize; i++) { cout <<setw(15)<< x[i] << "    "; } 
    cout <<endl;
  }
}

void TestWeightAxb(Matrix<double>& A, Vector<double>& b, 
                   Vector<double> &weight, 
                   Vector<double>& sol, unsigned solsize,
                   LeastSquaresLapack& lsla, LeastSquaresSVD lssvd)
{
  Vector<double> x;
  cout << "solution:   ";
  for (unsigned i=0; i<solsize; i++) { cout <<setw(15)<< sol[i] << "    "; }
  cout << endl;

  if (lsla.WeightedSolve(A, b, weight, x)==0){
    cout << "lapack sol: ";
    for (unsigned i=0; i<solsize; i++) { cout <<setw(15)<< x[i] << "    "; } 
    cout <<endl;
  }
  if (lssvd.WeightedSolve(A, b, weight, x)==0){ 
    cout << "svd sol:    ";
    for (unsigned i=0; i<solsize; i++) { cout <<setw(15)<< x[i] << "    "; } 
    cout <<endl;
  }

}

void TestNoWeightAx(Matrix<double>& A, 
                    Vector<double>& sol, unsigned solsize,
                    LeastSquaresLapack& lsla, LeastSquaresSVD lssvd)
{
  Vector<double> x;

  cout << "solution:   ";
  for (unsigned i=0; i<solsize; i++) { cout <<setw(15)<<sol[i] << "    "; }
  cout << endl;
  
  if (lsla.Solve(A, x)==0){
    cout << "lapack sol: ";
    for (unsigned i=0; i<solsize; i++) { cout<<setw(15)<<x[i]<<"    "; } 
    cout <<endl;
  }
  if (lssvd.Solve(A, x)==0){
    cout << "svd sol:    ";
    for (unsigned i=0; i<solsize; i++) { cout <<setw(15)<< x[i] << "    "; } 
    cout <<endl;
  }
}

void TestWeightAx(Matrix<double>& A, 
                  Vector<double> &weight, 
                  Vector<double>& sol, unsigned solsize,
                   LeastSquaresLapack& lsla, LeastSquaresSVD lssvd)
{
  Vector<double> x;
  cout << "solution:   ";
  for (unsigned i=0; i<solsize; i++) { cout <<setw(15)<< sol[i] << "    "; }
  cout << endl;

  if (lsla.WeightedSolve(A, weight, x)==0){
    cout << "lapack sol: ";
    for (unsigned i=0; i<solsize; i++) { cout <<setw(15)<< x[i] << "    "; } 
    cout <<endl;
  }
  if (lssvd.WeightedSolve(A, weight, x)==0){ 
    cout << "svd sol:    ";
    for (unsigned i=0; i<solsize; i++) { cout <<setw(15)<< x[i] << "    "; } 
    cout <<endl;
  }

}

//int main(int argc, char *argv[])
int main()
{
  Matrix<double> A;
  Vector<double> b, weight, sol, x;
  unsigned solsize=5, numeq=10;
 
  LeastSquaresLapack lsla;
  LeastSquaresSVD lssvd;

  GetSystem(A, b, weight, sol, solsize, numeq);

  lsla.Init(solsize, true);
  lssvd.Init(solsize, true);

  cout << "\nno weights, reduction to ATA, Ax=b\n";
  TestNoWeightAxb(A, b, sol, solsize, lsla, lssvd);

  cout << "\nweights, reduction to ATA, Ax=b\n";
  TestWeightAxb(A, b, weight, sol, solsize, lsla, lssvd);
 
  lsla.Init(solsize, false);
  lssvd.Init(solsize, false);

  cout << "\nno weights, no reduction to ATA, Ax=b\n";
  TestNoWeightAxb(A, b, sol, solsize, lsla, lssvd);

  cout << "\nweights, no reduction to ATA, Ax=b\n";
  TestWeightAxb(A, b, weight, sol, solsize, lsla, lssvd);



  GetSystem(A, weight, sol, solsize, numeq);

  lsla.Init(solsize, true);
  lssvd.Init(solsize, true);

  cout << "\nno weights, reduction to ATA, Ax=0\n";
  TestNoWeightAx(A, sol, solsize, lsla, lssvd);

  cout << "\nweights, reduction to ATA, Ax=0\n";
  TestWeightAx(A, weight, sol, solsize, lsla, lssvd);
 
  lsla.Init(solsize, false);
  lssvd.Init(solsize, false);

  cout << "\nno weights, no reduction to ATA, Ax=0\n";
  TestNoWeightAx(A, sol, solsize, lsla, lssvd);

  cout << "\nweights, no reduction to ATA, Ax=0\n";
  TestWeightAx(A, weight, sol, solsize, lsla, lssvd);

  
 

 return 0;
}