PlayaBasicLMBFGS.cpp
00001 #include "PlayaBasicLMBFGS.hpp"
00002 #include "PlayaOut.hpp"
00003 #include "PlayaTabs.hpp"
00004 #include "PlayaLinearCombinationImpl.hpp"
00005 #include "PlayaLineSearchBuilder.hpp"
00006 #include "PlayaOptConvergenceTestBuilder.hpp"
00007
00008 namespace Playa
00009 {
00010 using std::endl;
00011
00012 BasicLMBFGS::BasicLMBFGS(
00013 const ParameterList& params
00014 )
00015 : LineSearchBasedOptBase(params),
00016 memSize_(getParameter<int>(params, "Max Memory Size"))
00017 {}
00018
00019
00020 RCP<DirectionGeneratorBase>
00021 BasicLMBFGS::makeDirectionGenerator() const
00022 {
00023 return rcp(new BasicLMBFGSDirection(memSize_));
00024 }
00025
00026
00027 BasicLMBFGSDirection::BasicLMBFGSDirection(int memSize)
00028 : memSize_(memSize),
00029 xPrev_(),
00030 gradPrev_(),
00031 sMem_(),
00032 yMem_()
00033 {}
00034
00035 bool BasicLMBFGSDirection::generateDirection(
00036 const RCP<ObjectiveBase>& obj,
00037 const Vector<double>& xCur,
00038 const Vector<double>& gradCur,
00039 const double& fCur,
00040 Vector<double>& p)
00041 {
00042 Vector<double> q = gradCur.copy();
00043 int numStored = sMem_.size();
00044 Array<double> alpha(numStored);
00045 Array<double> rho(numStored);
00046
00047 for (int i=numStored-1; i>=0; i--)
00048 {
00049 rho[i] = 1.0/(sMem_[i]*yMem_[i]);
00050 alpha[i] = rho[i] * (sMem_[i]*q);
00051 q = q - alpha[i]*yMem_[i];
00052 }
00053
00054 double gamma;
00055 if (numStored > 0)
00056 {
00057 int j = numStored-1;
00058 gamma = (sMem_[j]*yMem_[j])/(yMem_[j]*yMem_[j]);
00059 }
00060 else
00061 {
00062 gamma = obj->getInvHScale();
00063 }
00064
00065 Vector<double> r = gamma*q;
00066
00067 for (int i=0; i<numStored; i++)
00068 {
00069 double beta = rho[i]*(yMem_[i]*r);
00070 r = r + (alpha[i]-beta)*sMem_[i];
00071 }
00072
00073 p = -1.0*r;
00074
00075 if (xPrev_.ptr().get() != 0)
00076 {
00077 Vector<double> s = xCur - xPrev_;
00078 Vector<double> y = gradCur - gradPrev_;
00079 sMem_.push_back(s);
00080 yMem_.push_back(y);
00081 if ((int) sMem_.size() > memSize_)
00082 {
00083 sMem_.pop_front();
00084 yMem_.pop_front();
00085 }
00086 }
00087
00088 xPrev_.acceptCopyOf(xCur);
00089 gradPrev_.acceptCopyOf(gradCur);
00090
00091 return true;
00092 }
00093
00094
00095 }