LeastSquares.hh

/* 
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
*/


#ifndef __LeastSquares_hh__
#define __LeastSquares_hh__

#include <Base/Math/Matrix.hh>
#include <Base/Math/Vector.hh>
#include "SVD.hh"


#define D_LS_MATRIXES      0x00000001
#define D_LS_RED_MATRIXES  0x00000002
#define D_LS_SVD           0x00000004

namespace BIAS {

  /** @class LeastSquaresBase
      @ingroup g_mathalgo
      @brief Base class for linear least squares solvers.
      @author MIP
  */
  class BIASMathAlgo_EXPORT LeastSquaresBase : public Debug
  {
  public:

    LeastSquaresBase() : _ATA(), _AT(), _Weight(), _ATb() {};

    virtual ~LeastSquaresBase() {};
    
    /** @brief Initialize least squares solver.
        @param[in] SolutionSize Number of parameters
        @param[in] ReduceToATA Reduce minimization of |Ax-b| to solution of
                               (A'A)x = A'b
    */
    int Init(unsigned SolutionSize, bool ReduceToATA = true);

    /** @brief Solve linear equation system A*x = b using linear least squares,
               i.e. minimize |A*x - b|^2.
        @return Returns 0 in case of success, <0 otherwise.
    */
    virtual int Solve(Matrix<double> &A, Vector<double> &b, Vector<double> &x);

    /** @brief Solve linear equation system A*x = 0 using linear least squares,
               i.e. minimize |A*x|^2, subject to |x| = 1.
        @return Returns 0 in case of success, <0 otherwise.
    */
    virtual int Solve(Matrix<double> &A, Vector<double> &x);

    /** @brief Solve linear equation system A*x = b using weighted linear least
               squares, i.e. minimize |W*(A*x - b)|^2 with diagonal matrix W
               composed of given weights.
        @return Returns 0 in case of success, <0 otherwise.
    */
    virtual int WeightedSolve(Matrix<double> &A, Vector<double> &b, 
                              Vector<double> &weights, Vector<double> &x);

    /** @brief Solve linear equation system A*x = 0 using weighted linear least
               squares, i.e. minimize |W*(A*x - b)|^2 with diagonal matrix W
               composed of given weights, subject to |x| = 1.
        @return Returns 0 in case of success, <0 otherwise.
    */
    virtual int WeightedSolve(Matrix<double> &A, Vector<double> &weights, 
                              Vector<double> &x);

  protected:

    bool _ReduceToATA;
    unsigned _SolutionSize;
    Matrix<double> _ATA, _AT, _Weight;
    Vector<double> _ATb;

  }; // class LeastSquaresBase


  /** @class LeastSquaresLapack
      @ingroup g_mathalgo
      @brief Linear least squares solver based on Lapack routines.
      @author MIP
  */
  class BIASMathAlgo_EXPORT LeastSquaresLapack : public LeastSquaresBase
  {
  public:

    LeastSquaresLapack() {};

    virtual ~LeastSquaresLapack() {};
    
    virtual int Solve(Matrix<double> &A, Vector<double> &b, Vector<double> &x);

    virtual int Solve(Matrix<double> &A, Vector<double> &x);

    virtual int WeightedSolve(Matrix<double> &A, Vector<double> &b, 
                              Vector<double> &weights, Vector<double> &x);

    virtual int WeightedSolve(Matrix<double> &A, Vector<double> &weights, 
                              Vector<double> &x);

  }; // class LeastSquaresLapack


  /** @class LeastSquaresSVD
      @ingroup g_mathalgo
      @brief Linear least squares solver based on singular value decomposition.
      @author MIP
  */
  class BIASMathAlgo_EXPORT LeastSquaresSVD : public LeastSquaresBase
  {
  public:

    LeastSquaresSVD() : _svd(), _VT(), _S() {};

    virtual ~LeastSquaresSVD() {};
    
    virtual int Solve(Matrix<double> &A, Vector<double> &b, Vector<double> &x);

    virtual int Solve(Matrix<double> &A, Vector<double> &x);

    virtual int WeightedSolve(Matrix<double> &A, Vector<double> &b, 
                              Vector<double> &weights, Vector<double> &x);

    virtual int WeightedSolve(Matrix<double> &A, Vector<double> &weights, 
                              Vector<double> &x);

  protected:

    SVD _svd;
    Matrix<double> _VT;
    Vector<double> _S;

  }; // class LeastSquaresSVD

} // namespace

#endif // __LeastSquares_hh__