---

leastsquare.h

// Emacs will be in -*- Mode: c++ -*-
//
// ************ DO NOT REMOVE THIS BANNER ****************
//
//  Nicolas Di Cesare <Nicolas.Dicesare@free.fr>
//  http://acm.emath.fr/~dicesare
//
//********************************************************
//
//  Least Square Method using conjugate gradient method
//
//********************************************************
#ifndef __leastsquare_h__
#define __leastsquare_h__

#include <cg.h>

template <class V, class M>
class LeastSquareProblem {
public:
  virtual int size()  = 0;
  virtual void targetAndValue(const V& x, V& target, V& fct2fit) = 0;
  virtual void targetValueAndfirstDerivative(const V& x, M& grad_fct2fit, V& target, V& fct2fit) = 0;
};

template <class V, class M>
class LeastSquareFunction : public CostFunction<V> {
public:
  typedef double value_type;
protected:
  LeastSquareProblem<V,M> & lsp_;
public:
  LeastSquareFunction(LeastSquareProblem<V,M> & lsp) : lsp_(lsp) {}
  virtual ~LeastSquareFunction() {}
  
  virtual value_type value(const V& x);
  virtual void firstDerivative(V& grad_f, const V& x);
  virtual value_type valueAndFirstDerivative(V& grad_f, const V& x);
  
  virtual void Update() {}
  virtual void Save(OptimizationMethod<V>&) {}
};


template <class V, class M>
typename LeastSquareFunction<V,M>::value_type 
LeastSquareFunction<V,M>::value(const V& x) 
{
  // size of target and function to fit vectors
  V target(lsp_.size()), fct2fit(lsp_.size()); 
  // compute its values
  lsp_.targetAndValue(x,target,fct2fit);
  // do the difference 
  V diff = target - fct2fit;
  // and compute the scalar product (square of the norm)
  return DotProduct(diff,diff);
}

template <class V, class M>
void LeastSquareFunction<V,M>::firstDerivative(V& grad_f, const V& x)
{
  // size of target and function to fit vectors
  V target(lsp_.size()), fct2fit(lsp_.size()); 
  // size of gradient matrix
  M grad_fct2fit(lsp_.size(),x.size());
  // compute its values
  lsp_.targetValueAndfirstDerivative(x,grad_fct2fit,target,fct2fit);
  // do the difference 
  V diff = target - fct2fit;
  // compute derivative
  grad_f = -2.*(transpose(grad_fct2fit)*diff);
}

template <class V, class M>
typename LeastSquareFunction<V,M>::value_type 
LeastSquareFunction<V,M>::valueAndFirstDerivative(V& grad_f, const V& x) 
{
  // size of target and function to fit vectors
  V target(lsp_.size()), fct2fit(lsp_.size()); 
  // size of gradient matrix
  M grad_fct2fit(lsp_.size(),x.size());
  // compute its values
  lsp_.targetValueAndfirstDerivative(x,grad_fct2fit,target,fct2fit);
  // do the difference 
  V diff = target - fct2fit;
  // compute derivative
  grad_f = -2.*(transpose(grad_fct2fit)*diff);
  // and compute the scalar product (square of the norm)
  return DotProduct(diff,diff);
}



template <class V>
class NonLinearLeastSquare {
  V results_, initialValue_;
  double resnorm_;
  int exitFlag_;
  double accuracy_, bestAccuracy_;
  unsigned int maxIterations_, nbIterations_;
  QuantLib::Handle< OptimizationMethod<V> > om_;
public:
  inline NonLinearLeastSquare(double accuracy = 1e-4,
                  int maxiter = 100);
  inline NonLinearLeastSquare(double accuracy,
                  int maxiter,
                  QuantLib::Handle< OptimizationMethod<V> > om);  
  inline ~NonLinearLeastSquare() {}

  template <class M> inline 
  V& Perform(LeastSquareProblem<V,M> & lsProblem)
  {
    double eps = accuracy_;

    // set initial value of the optimization method
    om_->setInitialValue(initialValue_);
    // set end criteria with a given maximum number of iteration and a given error eps
    om_->setEndCriteria(OptimizationEndCriteria(maxIterations_, eps));
    om_->endCriteria().setPositiveOptimization();

    // wrap the least square problem in an optimization function
    LeastSquareFunction<V,M> lsf(lsProblem);

    // define optimization problem
    OptimizationProblem<V> P(lsf,*om_);

    // minimize
    P.Minimize();

    // summarize results of minimization
    exitFlag_ = om_->endCriteria().criteria();
    nbIterations_ = om_->iterationNumber();

    results_ = om_->x();
    resnorm_ = om_->functionValue();
    bestAccuracy_ = om_->functionValue();
  
    return results_;
  }

  inline void setInitialValue(const V& initialValue)
  { initialValue_ = initialValue;}

  inline V& results() { return results_;}
  inline double residualNorm() { return resnorm_;}
  inline double lastValue() { return bestAccuracy_;}
  inline int exitFlag() { return exitFlag_;}
  inline int iterationsNumber() { return nbIterations_;}
};


template <class V> inline 
NonLinearLeastSquare<V>::NonLinearLeastSquare(double accuracy,
                          int maxiter)
  : exitFlag_(-1), accuracy_(accuracy), maxIterations_(maxiter), 
    om_(QuantLib::Handle< OptimizationMethod<V> >(new ConjugateGradient<V>()))
{}

template <class V> inline 
NonLinearLeastSquare<V>::NonLinearLeastSquare(double accuracy,
                          int maxiter,
                          QuantLib::Handle< OptimizationMethod<V> > om)
  : exitFlag_(-1), accuracy_(accuracy), maxIterations_(maxiter), om_(om)
{}


#endif

---