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 }