---

ConjugateGradientMinimizer.cpp

Go to the documentation of this file.

// ConjugateGradientMinimizer.cpp: implementation of the CConjugateGradientMinimizer class.
//

#include "ConjugateGradientMinimizer.h"

// Construction/Destruction

const double CConjugateGradientMinimizer::m_dGOLD = 1.618034;


CConjugateGradientMinimizer::CConjugateGradientMinimizer(bool searchFlag, CParameterFilter *pFilter, int nIterations, double funcAccuracy, int nRestarts, double gradTol, double funcTol)
{
        m_pFilter = pFilter;
        m_iNIterations = nIterations;
        m_iNRestarts = nRestarts;
        m_dFuncAccuracy = funcAccuracy;
        m_dGradTol = gradTol;
        m_dFuncTol = funcTol;
        m_pMO = new CMatrixOperations();
        m_dMinTol = m_dFuncAccuracy;
        m_bSearchFlag = searchFlag;
}

CConjugateGradientMinimizer::~CConjugateGradientMinimizer()
{
        delete m_pMO;
}

double CConjugateGradientMinimizer::Minimize(double *parameters, Minimizable *minimizable)
{
        // loop variables
        int i,j;
        double fbest = 0.0;
        double GammaNum = 0.0;
        double GammaDen = 0.0;
        double Gamma = 0.0;
        double xA = 0.0;
        double xB = 0.0;
        double xC = 0.0;
        double xM = 0.0;
        double fA = 0.0;
        double fB = 0.0;
        double fC = 0.0;
        double fM = 0.0;
        double supNorm = 0.0;
        // set local minimizable pointer
        m_pMinimizable = minimizable;
        // set number of parameters
        nParameters = m_pMinimizable->GetNParameters();
        // allocate
        // all parameters/gradients/etc. are stored in their transformed (filtered) form
        m_pdCurrentParameters = new double[nParameters];
        m_pdSearchDirection = new double[nParameters];
        m_pdGradient = new double[nParameters];
        double *oldGrad = new double[nParameters];
        m_pdForceVector = new double[nParameters];
        
        // copy filtered parameters into current p's
        m_pFilter->ForwardTransformation(parameters,nParameters);
        m_pMO->ElementCopy(parameters,m_pdCurrentParameters,nParameters);
        m_pFilter->BackwardTransformation(parameters,nParameters);

        
        // now get the gradient at the initial point
        this->ComputeGradient(m_pdCurrentParameters);
        
        // search direction = -g if none is provided
        if(m_bSearchFlag == false)
        {
                // set the initial search direction equal to -g
                m_pMO->ElementCopy(m_pdGradient,m_pdSearchDirection,nParameters);
                m_pMO->RescaleVector(m_pdSearchDirection,-1.0,nParameters);
        }
        else // search direction is read from a file
        {
                ParameterReader *pReader = new ParameterReader("search.par");
                for(j = 0; j < nParameters; j++)
                {
                        m_pdSearchDirection[j] = pReader->ReadParameter();
                }
                delete pReader;
        }
        
        xA = 0.0;
        // use 1/(sup norm) for the first iteration
        supNorm = m_pMO->VectorSupNorm(m_pdSearchDirection,nParameters);
        xB = 1.0/supNorm;
        
        int restartCount = 0;
        bool restartFlag = false;
        // now iterate
        for(i = 0; i < m_iNIterations; i++)
        {
                // evaluate f at p and p + 1.0*s
                fA = this->LineEvaluate(xA);
                cout << "Cost at beginning of iteration " << i << " is " << fA << endl;

                fB = this->LineEvaluate(xB);
                cout << "Cost in conjugate direction is " << fB << endl;

                // bracket the minimum
                this->BracketMinimum(&xA,&fA,&xB,&fB,&xC,&fC);
                // now line minimize in the bracket to find the minimum
                this->BrentLine(xA,xB,xC,&xM,&fM);
                // set the current set of parameters equal to p + xM*s
                m_pMO->TwoVectorLC(1.0,m_pdCurrentParameters,xM,m_pdSearchDirection,m_pdCurrentParameters,nParameters);
                fbest = fM;

                // write the accepted parameters to a file so we can save
                // incremental progress
                fstream *pftemp = new fstream("cgtemp.par",ios::out);
                *pftemp << "##################################################################" << endl;
                *pftemp << "# BEST PARAMETERS" << endl;
                *pftemp << "# COST : " << fbest << endl;
                *pftemp << "##################################################################" << endl;
                int pcount;
                for(pcount = 0; pcount < nParameters; pcount++)
                {
                        *pftemp << setprecision(10) << m_pFilter->OperatorInverse(m_pdCurrentParameters[pcount]) << endl;
                }
                delete pftemp;
                // END WRITE

                // save the old gradient and compute the new gradient
                m_pMO->ElementCopy(m_pdGradient,oldGrad,nParameters);
                this->ComputeGradient(m_pdCurrentParameters);
                // now compute gamma
                GammaNum = m_pMO->DotProduct(m_pdGradient,m_pdGradient,nParameters);
                // if the magnitude of the new gradient is small, terminate
                if(this->CheckGradTol(sqrt(GammaNum)) ) {break;}
                if(this->CheckFuncTol(fA - fM)) 
                {
                        // convergence suspected. do a restart
                        if(restartCount < m_iNRestarts)
                        {
                                cout << "Restart number " << restartCount << endl;
                                // this combined with later statements restarts along the
                                // gradient direction
                                GammaNum = 0.0;
                                GammaDen = 1.0;
                                restartCount++;
                                restartFlag = true;
                        }
                        else
                        {
                                // out of restarts
                                break;
                        }
                }
                else
                {
                        // continue computing the new search direction
                        GammaNum -= m_pMO->DotProduct(oldGrad,m_pdGradient,nParameters);
                        GammaDen = m_pMO->DotProduct(oldGrad,oldGrad,nParameters);
                }       
                Gamma = GammaNum/GammaDen;
                // set the new search direction
                m_pMO->TwoVectorLC(-1.0,m_pdGradient,Gamma,m_pdSearchDirection,m_pdSearchDirection,nParameters);
                
                // write the new search direction to a file
                fstream *pfsrch = new fstream("searchtemp.par",ios::out);
                *pfsrch << "##################################################################" << endl;
                *pfsrch << "# SEARCH DIRECTION" << endl;
                *pfsrch << "##################################################################" << endl;
                for(pcount = 0; pcount < nParameters; pcount++)
                {
                        *pfsrch << setprecision(10) << m_pdSearchDirection[pcount] << endl;
                }
                delete pfsrch;
                
                xA = 0.0;
                // for the step in the next iteration, use either the sup norm or the distance to 
                // the previous minimium, whichever is smaller
                supNorm = m_pMO->VectorSupNorm(m_pdSearchDirection,nParameters);
                if( (fabs(xM) < 1.0/supNorm) && (!restartFlag) )
                {
                        xB = fabs(xM);
                        cout << "Distance to last minimum used to set scale of next step." << endl;
                }
                else
                {
                        xB = 1.0/supNorm;
                        cout << "Sup norm used to set scale of next step." << endl;
                        restartFlag = false;
                }
        }

        // return best parameters (which are the current ones) in parameter array,
        // remembering to transform
        m_pFilter->BackwardTransformation(m_pdCurrentParameters,nParameters);
        m_pMO->ElementCopy(m_pdCurrentParameters,parameters,nParameters);

        // cleanup
        delete [] m_pdCurrentParameters;
        delete [] m_pdSearchDirection;
        delete [] m_pdGradient;
        delete [] oldGrad;

        return fbest;
}

void CConjugateGradientMinimizer::ComputeGradient(double *parameters)
{
        int i;
        double saveparameter = 0.0;
        double fPlus = 0.0;
        double fMinus = 0.0;
        
        // loop over parameters
        for(i = 0; i < nParameters; i++)
        {
                // parameters are already transformed, so just save the old value
                saveparameter = parameters[i];

                double h = pow(m_dFuncAccuracy,0.333333)*fabs(parameters[i]);
                // if statement generates warning if step is really small
                if(h < DBL_EPSILON)
                {
                        cout << "WARNING:Parmeter step small in ConjugateGradientMinimizer::ComputeGradient!" << endl;
                        if(fabs(parameters[i]) < DBL_MIN)
                        {
                                cout << "WARNING:Zero parameter detected, truncating . . ." << endl;
                                // pretend like the parameter was really small
                                double changedpar = m_pFilter->OperatorInverse(parameters[i]) + DBL_EPSILON;
                                h = pow(m_dFuncAccuracy,0.333333)*fabs(m_pFilter->Operator(changedpar));
                        }
                        
                }
                                
                // forward step
                parameters[i] += h;
                m_pFilter->BackwardTransformation(parameters,nParameters);
                fPlus = m_pMinimizable->ObjectiveFunction(parameters);
                m_pFilter->ForwardTransformation(parameters,nParameters);
                parameters[i] = saveparameter;
                // backward step
                parameters[i] -= h;
                m_pFilter->BackwardTransformation(parameters,nParameters);
                fMinus = m_pMinimizable->ObjectiveFunction(parameters);
                m_pFilter->ForwardTransformation(parameters,nParameters);

                // restore original parameter           
                parameters[i] = saveparameter;

                // fill in the appropriate piece of the gradient
                m_pdGradient[i] = (fPlus - fMinus)/(2*h);
                m_pdForceVector[i] = -1.0*m_pdGradient[i];
        }

    //print vector of forces
        Notify();

        cout << "Gradient calculated . . . " << endl;

        return;
}

void CConjugateGradientMinimizer::BracketMinimum(double *xA, double *fA, double *xB, double *fB, double *xC, double *fC)
{
        double xU,fU,r,q,xUlim;
        double glim = 100.0;

        // if f(B) > f(A), reverse the roles of A and B
        if(*fB > *fA)
        {
                double temp;
                temp = *xA;
                *xA = *xB;
                *xB = temp;
                temp = *fA;
                *fA = *fB;
                *fB = temp;
        }
        // initial guess
        *xC = *xB + m_dGOLD*(*xB - *xA);
        *fC = this->LineEvaluate(*xC);
        // now loop
        while(*fB > *fC)
        {
                r = (*xB - *xA)*(*fB - *fC);
                q = (*xB - *xC)*(*fB - *fA);
                xU = *xB - ((*xB - *xC)*q - (*xB - *xA)*r)/(2.0*(q-r));
                // largest value that we're willing to attempt
                xUlim = *xB + glim*(*xC - *xB);
                
                if((*xB - xU)*(xU - *xC) > 0.0 )  // parabolic u is in (b,c)
                {
                        // get f(U)
                        fU = this->LineEvaluate(xU);
                        if(fU < *fC)  // got a min. between b and c
                        {
                                *xA = *xB;
                                *xB = xU;
                                *fA = *fB;
                                *fB = fU;
                                return;
                        }
                        else if(fU > *fC)       // got one between a and u
                        {
                                *xC = xU;
                                *fC = fU;
                                return;
                        }

                        // parabolic fit was useless, take another step forward
                        xU = *xC + m_dGOLD*(*xB - *xC);
                        fU = this->LineEvaluate(xU);
                }
                else if((*xC - xU)*(xU - xUlim) > 0.0) // parabolic u is in (c,ulim)
                {
                        fU = this->LineEvaluate(xU);
                        if(fU < *fC)
                        {
                                *xB = *xC;
                                *xC = xU;
                                xU = *xC + m_dGOLD*(*xC - *xB);
                                *fB = *fC;
                                *fC = fU;
                                fU = this->LineEvaluate(xU);
                        }
                }
                else if((xU - xUlim)*(xUlim - *xC) >= 0.0 ) // limit parabolic u to its max
                {
                        xU = xUlim;
                        fU = this->LineEvaluate(xU);
                }
                else    // reject parabolic u entirely, use default magnification
                {
                        xU = *xC + m_dGOLD*(*xC - *xB);
                        fU = this->LineEvaluate(xU);
                }
                // dump the oldest point and continue
                *xA = *xB;
                *xB = *xC;
                *xC = xU;
                *fA = *fB;
                *fB = *fC;
                *fC = fU;
        }       // end while loop
}

void CConjugateGradientMinimizer::BrentLine(double xA, double xB, double xC, double *xM, double *fM)
{
        int i;
        double u,x,w,v,xm;
        double fu,fx,fw,fv;
        double a,b;
        double p,q,r,d,etemp;
        double tol1,tol2;
        double ZEPS = 1.0e-10;
        double RATIO = 1.0 - (1.0/m_dGOLD);

        double e = 0.0;
        
        // a and b need to be in ascending order, but the input abcissas need 
        // not be
        if(xA < xC)
        {
                a = xA;
                b = xC;
        }
        else
        {
                a = xC;
                b = xA;
        }
        x = xB;
        w = xB;
        v = xB;
        fx = this->LineEvaluate(x);
        fw = fx;
        fv = fx;

        for(i = 0; i < 100; i++)
        {
                xm = 0.5*(a + b);
                tol1 = m_dMinTol*fabs(x) + ZEPS;
                tol2 = 2.0*tol1;
                
                if(fabs(x - xm) <= (tol2 - 0.5*(b - a)) ) // test for termination
                {
                        *xM = x;
                        *fM = fx;
                        return;
                }
                if(fabs(e) > tol1) // construct a trial parabolic fit
                {
                        r = (x - w)*(fx - fv);
                        q = (x - v)*(fx - fw);
                        p = (x - v)*q - (x - w)*r;
                        q = 2.0*(q - r);
                        if(q > 0.0) {p = -1.0*p;}
                        q = fabs(q);
                        etemp = e;
                        e = d;
                        // determine the acceptability of the parabolic fit
                        if( (fabs(p) >= fabs(0.5*q*etemp)) || (p <= q*(a-x)) || (p >= q*(b-x)) )
                        {
                                if(x >= xm)
                                {
                                        e = a - x;
                                }
                                else
                                {
                                        e = b - x;
                                }
                                d = RATIO*e;
                        }
                        else  // take golden section step into the larger of the two segments
                        {
                                d = p/q;
                                u = x + d;
                                if( ((u - a) < tol2) || ((b - u) < tol2) )
                                {
                                        d = SgnProduct(tol1,xm - x);
                                }
                        }
                } // end if
                else
                {
                        if(x >= xm)
                        {
                                e = a - x;
                        }
                        else
                        {
                                e = b - x;
                        }
                        d = RATIO*e;
                }
                
                if(fabs(d) >= tol1)
                {
                        u = x + d;
                }
                else
                {
                        u = x + SgnProduct(tol1,d);
                }
                fu = this->LineEvaluate(u);

                // now decide what to do with u
                if(fu <= fx)
                {
                        if(u >= x)
                        {
                                a = x;
                        }
                        else
                        {
                                b = x;
                        }
                        v = w;
                        w = x;
                        x = u;
                        fv = fw;
                        fw = fx;
                        fx = fu;
                }
                else
                {
                        if(u < x)
                        {
                                a = u;
                        }
                        else
                        {
                                b = u;
                        }

                        if((fu <= fw) || (w == x))
                        {
                                v = w;
                                w = u;
                                fv = fw;
                                fw = fu;
                        }
                        else if((fu <= fv) || (v == x) || (v == w))
                        {
                                v = u;
                                fv = fu;
                        }
                }
        
        
        } // end for loop

        cout << "WARNING : Number of iterations in BrentLine() exceeded." << endl;

        return;
}

bool CConjugateGradientMinimizer::CheckGradTol(double norm)
{
        if(norm < m_dGradTol)
        {
                cout << "||g|| < eps_g, terminating . . . " << endl;
                return true;
        }
        else
        {
                return false;
        }
}

bool CConjugateGradientMinimizer::CheckFuncTol(double diff)
{
        if(diff < m_dFuncTol)
        {
                cout << "f_(k) - f_(k+1) < eps_f, convergence suspected . . . " << endl;
                return true;
        }
        else
        {
                return false;
        }
}

double CConjugateGradientMinimizer::LineEvaluate(double alpha)
{
        double funcVal = 0.0;
        double *z = new double[nParameters];
        // form the linear combination
        m_pMO->TwoVectorLC(1.0,m_pdCurrentParameters,alpha,m_pdSearchDirection,z,nParameters);
        // back transform
        m_pFilter->BackwardTransformation(z,nParameters);
        // evaluate
        funcVal = m_pMinimizable->ObjectiveFunction(z);
        // cleanup
        delete [] z;

        return funcVal;
}

double CConjugateGradientMinimizer::SgnProduct(double a, double b)
{
        int sign;
        if(b >= 0.0) {sign = 1;}
        else {sign = -1;}
        return sign*fabs(a);
}

---