Curve Fitting and Optimization





Doxygen documentation


Table of contents

I- Least Square : Curve Fitting example
II- Optimization
III- Optimization Example

I- Least Square : Curve Fitting example

(see CurveFitting.cpp)

We present here an example of curve fitting. Curve fitting is a particular case of least square problems. We solve least square problem using conjugate gradient algorithm. The most interesting part is the design of optimization problem (see next section). Provided classes are just a way to turn a least square problem into an optimization one.

The example is an inverse problem : we know the solution and we try to recover it by an optimization process. In our example, we want to fit the Brace volatility term structure function ( similar to Nelson-Siegel function for zero rate) to given data. These data are previously computed, using the function, to ensure the existence of a solution. The Brace function is defined with 4 parameters a, b, c and d.

The function (operator() (double t)) and its partial derivatives (da(double t), db(double t), dc(double t) and dd(double t)) values are defined int the class BraceVolatilityFunction.

 class BraceVolatilityFunction {
   double a_, b_, c_, d_;
 public :
   inline BraceVolatilityFunction(double a, double b, double c, double d)
     : a_(a), b_(b), c_(c), d_(d) {}
   inline ~BraceVolatilityFunction() {}

   inline double operator() (double t) 
   { return (a_+b_*t)*exp(-c_*t) + d_; }
 
   inline double da(double t) 
   { return exp(-c_*t); }
 
   inline double db(double t) 
   { return t*exp(-c_*t); }
 
   inline double dc(double t) 
   { return -t*(a_+b_*t)*exp(-c_*t); }
 
   inline double dd(double t) 
   { return 1.; }
 };


I first show how to proceed to solve the problem and later how to define it.

Let targetValue be the vector of values to fit and let timeToMaturity the the vector of time to maturity values. First, we define the accuracy of our optimization process and the maximum number of iteration maxiter. The optimization process ends when the asked accuracy is obtained or if the number of iteration is larger than maxiter. Finally, we define a vector of initial values initialValue set to constant value 0.1. These quantities are the inputs of the least square problem.

Then we define the least square solver NonLinearLeastSquare<Array > with given accuracy and maximum number of iterations. We also define the least square problem CurveFittingProblem with given time to maturity timeToMaturity and target values targetValue vectors.

We set the initial value in the solver lsqnonlin.setInitialValue(initialValue) and we finally perform the optimization by calling the member Perform on the CurveFittingProblem cfp. The member returns the solution vector.

 // Accuracy of the optimization method
 double accuracy = 1e-10;// It is the square of the accuracy
 // Maximum number of iterations
 int maxiter = 10000;
 
 Array initialValue(4, 0.1);// Initial value
 
 // Least square optimizer
 NonLinearLeastSquare<Array > lsqnonlin(accuracy,maxiter);
 
 // Define the least square problem
 CurveFittingProblem cfp(timeToMaturity, targetValue);
 
 // Set initial values
 lsqnonlin.setInitialValue(initialValue);
 // perform fitting
 Array solution = lsqnonlin.Perform(cfp);
 
 cout << endl 
      << "Optimal values  : " << abcd
      << "Solution values : " << solution
      << "Solution Error  : " << (solution - abcd ) << endl;


Thus, to perform the curve fitting task, we have to define the underlying least-square problem. It is done by the definition of class CurveFittingProblem deriving from the abstract class LeastSquareProblem<Array,Matrix >.

 class CurveFittingProblem 
     : public LeastSquareProblem<Array,Matrix > {
     const Array &ttm_;
     const Array &target_;
 public :
   CurveFittingProblem(const Array& ttm, const Array& target)
       : ttm_(ttm), target_(target) {}
   virtual ~CurveFittingProblem() {}

   virtual int size() { return ttm_.size(); }

   virtual void targetAndValue(const Array& abcd, 
                               Array& target, 
                               Array& fct2fit);
 
   virtual void targetValueAndfirstDerivative(const Array& abcd, 
 };


To fully satisfied the requirements of the abstract class LeastSquareProblem, we have to define the following members :
  • int size() : return the size of the least square problem. It is the size of the vector that contains data to fit (target values).
  • void targetAndValue(const Array& abcd, Array& target, Array& fct2fit) : return function and target values.
  • void targetValueAndfirstDerivative(const Array& abcd, Matrix& grad_fct2fit, Array& target, Array& fct2fit) : return function, target and first derivatives values.


 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;
 };

II- Optimization

It is an introduction to the design of multi-dimensional optimization classes.

The goal is to provide a maximum of flexibilty with a common way to proceed. By flexibility, I mean that we can define new optimization problem, new optimization method and new line search by simply respecting some constraints (abstract class interfaces).

Here are the idea that lead to the abstract classes :
  • The Cost Function is the function to minimize. We have to define the value of the function at a given point. The first derivative (firstDerivative) can be provided or computed by finite difference.
  • The Optimization Method has to minimize an Optimization Problem.
  • The Optimization Problem is the minimization of a Cost Function with a given Optimization Method.
The first idea is the Cost Function, thus we have written an abstract class called CostFunction with four virtual members. The member value(const V& x) is "pure virtual", thus you have to overload it to define your cost function. It is not required to overload firstDerivative(V& grad_f, const V& x) that can be computed by finite difference but we highly encourage people to overload it with the analytic form of the first derivative to avoid truncation errors. For finite difference approximation of the first derivative, it's possible to specify the small quantity used for computation (epsilon). It's done by overloading finiteDifferenceEpsilon() which default value is 1e-8. Some time, it is usefull (faster) to compute the value and the first derivative of the function at the same time, thus we provide the member valueAndFirstDerivative(V& grad_f, const V& x) to do that.

 template <class V>
 class CostFunction {
 public:
   typedef double value_type;
 
   virtual value_type value(const V& x) = 0;
 
   virtual void firstDerivative(V& grad_f, const V& x)
   { /* default is finite difference computation */ }
 
   virtual value_type valueAndFirstDerivative(V& grad_f, const V& x) 
   { firstDerivative(grad_f, x); return value(x);}// default case
 
   virtual void Update() {}

   virtual double finiteDifferenceEpsilon() { return 1e-8;}
 }; 


The Optimization Method idea leads to OptimizationMethod abstract class. It only requires to provide (to overload) a member to minimize an Optimization Problem . We also store several useful informations : the iteration number, the number of evaluation of the function and the first derivative, the last function value, the last norm of the first derivative, the current and initial values of the unknown, the search direction and the optimization end criteria.

 template <class V>
 class OptimizationMethod {
 public:
   typedef double value_type;
 protected:
   V initialValue_;
   int iterationNumber_;
   OptimizationEndCriteria endCriteria_;
   int functionEvaluation_, gradientEvaluation_;
   value_type functionValue_, squaredNorm_;
   V x_, searchDirection_;
 public:
   explicit OptimizationMethod(); 
   virtual ~OptimizationMethod() {}
 
   inline void setInitialValue(const V& initialValue); 
   inline void setEndCriteria(const OptimizationEndCriteria& endCriteria); 

   inline int& iterationNumber() { return iterationNumber_;}
   inline OptimizationEndCriteria& endCriteria() { return endCriteria_;}
   inline int & functionEvaluation() { return functionEvaluation_; }
   inline int & gradientEvaluation() { return gradientEvaluation_; }
   inline value_type& functionValue() { return functionValue_; }
   inline value_type& gradientNormValue() { return squaredNorm_; }
   V& x() { return x_;}
   V& searchDirection() { return searchDirection_;}
   
   virtual void Minimize(OptimizationProblem<V> & P) = 0;
 };


Steepest Descent and Conjugate Gradient optimization methods are implemented using Armijo line search algorithm.

The class OptimizationProblem is the implementation of the last idea. It is the core of the optimization design. Given references to CostFunction and OptimizationMethod classes, it performs the optimization task.

 template <class V>
 class OptimizationProblem {
 public:
   typedef double value_type;
 protected:
   CostFunction<V> & costFunction_;
   OptimizationMethod<V> & method_;
   QuantLib::Handle<OptimizationProblemOutput<V> > opo_;
 public:
   OptimizationProblem(CostFunction<V> & f,
               OptimizationMethod<V> & meth);  
   OptimizationProblem(CostFunction<V> & f,
               OptimizationMethod<V> & meth,
               QuantLib::Handle<OptimizationProblemOutput<V> >& opo);
 
   ~OptimizationProblem() {}
   value_type value(const V& x);
   void firstDerivative(V& grad_f, const V& x);
   value_type valueAndFirstDerivative(V& grad_f, const V& x);

   OptimizationMethod<V> & optimisationMethod() { return method_;}
 
   void Minimize() { method_.Minimize( *this ); }
 
   V& minimumValue() { return method_.x();}
 
   void Save(int iterationNumber, 
         value_type function, 
         value_type normGradient,
         value_type lineSearchStep,
         OptimizationMethod<V>& method);
 };

III- Optimization Example

It's a one dimensional example, just increase the Array size to N and you will obtain a N dimensional example.

Definition of a new cost function

 template <class V>
 class xSquare : public CostFunction<V> {
 public:
   xSquare() {}
   virtual ~xSquare() {}
 
   virtual double value(const V& x)
   {
     double x0 = x[0];
     return x0*x0;
   }
 
   virtual void firstDerivative(V& grad_f, const V& x)
   { 
     double x0 = x[0];
     grad_f[0] = 2.*x0; 
   }
 
 };


Minimization example

 double eps = 1e-8;
 int maxIter = 10000;
 
 Array x0(1,4.);

 xSquare<Array> Square;
 
 double alpha = 0.5, beta = 0.65;
 Handle<LineSearch<Array> > armijoLineSearch = 
            new ArmijoLineSearch<Array>(eps,alpha,beta);
 ConjugateGradient<Array> conjugateGradient(armijoLineSearch);

 conjugateGradient.setInitialValue(x0);
 conjugateGradient.setEndCriteria(OptimizationEndCriteria(maxIter, eps) );
     
 OptimizationProblem<Array> P(Square,conjugateGradient);
 
 P.Minimize();
 
 std::cout << std::endl 
           << "criteria = " 
           << conjugateGradient.endCriteria().criteria() << std::endl
           << "Minimum value : " 
           << P.minimumValue();