GPLIB++: CQuadInterLineSearch.cpp Source File

00001 #include "CQuadInterLineSearch.h"
00002 #include <iostream>
00003 
00004 using namespace std; 
00005 CQuadInterLineSearch::CQuadInterLineSearch()
00006 {
00007         Wolfe1Factor = 1e-4;
00008 }
00009 
00010 CQuadInterLineSearch::~CQuadInterLineSearch()
00011 {
00012 }
00013 
00014 double CQuadInterLineSearch::SearchStepsize(const double initial,const double CurrentMisfit,
00015                 const ublas::vector<double> &Searchdirection, const ublas::vector<double> &Model,
00016                 const ublas::vector<double> &Data, const ublas::vector<double> Gradient)
00017 {
00018         double stepsize = initial;
00019         double optstepsize = 2. * initial;
00020         const int maxsearch = 20;
00021         int searchiterations = 0; 
00022         double NewMisfit1 = 1e90;
00023         double NewMisfit2 = 1e90;
00024         double OptMisfit = 1e90;
00025         ublas::vector<double> NewModel(Model.size());
00026         ublas::vector<double> OptModel(Model.size());
00027         double Wolfe1 = 0; 
00028         const double gradsearchprodukt =  inner_prod(Gradient,Searchdirection);
00029         NewModel= Model + stepsize * Searchdirection;
00030         NewMisfit1 = MisfitCalculator->CalcPerformance(NewModel);
00031         NewModel= Model + optstepsize * Searchdirection;
00032         NewMisfit2 = MisfitCalculator->CalcPerformance(NewModel);
00033     // Quadratic interpolation taken from
00034     // http://www.mathworks.com/access/helpdesk_r13/help/toolbox/optim/tutori5b.html
00035         double beta12 =  - pow(stepsize,2);
00036         double beta23 = pow(stepsize,2) - pow(optstepsize,2);
00037         double beta31 = pow(optstepsize,2);
00038         double gamma12 = -stepsize;
00039         double gamma23 = (stepsize - optstepsize);
00040         double gamma31 = optstepsize;
00041         optstepsize = 0.5 * (beta23*CurrentMisfit + beta31 * NewMisfit1 + beta12 * NewMisfit2)
00042                                                 /(gamma23*CurrentMisfit + gamma31 * NewMisfit1 + gamma12 * NewMisfit2);
00043         OptModel= Model + optstepsize * Searchdirection;
00044         OptMisfit = MisfitCalculator->CalcPerformance(OptModel);
00045         cout << "OptStepsize: " << optstepsize << " OptMisfit: " << OptMisfit << endl;
00046         Wolfe1 = CurrentMisfit + Wolfe1Factor * stepsize * gradsearchprodukt;  
00047         while ((OptMisfit > Wolfe1) && (searchiterations < maxsearch))
00048         {
00049                 NewModel= Model + stepsize * Searchdirection;
00050                 OptMisfit = MisfitCalculator->CalcPerformance(NewModel);
00051                 Wolfe1 = CurrentMisfit + Wolfe1Factor * stepsize * gradsearchprodukt;  
00052                 cout << "NewMisfit: " << OptMisfit << " ";
00053                 cout << "Wolfe1: " << Wolfe1 << endl;
00054                 stepsize /= 2.;
00055                 searchiterations++;
00056         }
00057         cout << "Searches: " << searchiterations << " Stepsize: " << 2 * stepsize << endl;
00058         return stepsize*2;
00059 }