---

LeastSquaresMinimizer.cpp

Go to the documentation of this file.

// LeastSquaresMinimizer.cpp: implementation of the LeastSquaresMinimizer class.
//

#include "LeastSquaresMinimizer.h"


// Construction/Destruction

CLeastSquaresMinimizer::CLeastSquaresMinimizer()
{

}

CLeastSquaresMinimizer::~CLeastSquaresMinimizer()
{

}


// Compute approximate Hessian for use in Least Squares algs.
// The hessian is symmetric so we only need compute 1/2 of it then
// copy the rest.  This routine assumes analytic second derivatives
// are not available, so finite differencing with a 3-pt central
// formula:
//                      f'(x) ~ (f(x+h) - f(x-h))/(2*h)
// is employed to obtain derivatives.

void CLeastSquaresMinimizer::ComputeDerivativeInformation(double *parameters, NLLSMinimizable *nlls)
{
        int i,j,k,nR;
        double transparameter = 0.0;
        double saveparameter = 0.0;
        double *fxPlus = new double[m_iNResiduals];
        double *fx = new double[m_iNResiduals];
        double *fxMinus = new double[m_iNResiduals];
        double **jacobian = new double *[m_iNResiduals];
        jacobian[0] = new double[m_iNResiduals*nParameters];
        for(nR = 1; nR < m_iNResiduals; nR++)
        {
                jacobian[nR] = &(jacobian[0][nR*nParameters]);
        }

        double chiSq = nlls->ObjectiveFunction(parameters);
        for(k = 0; k < m_iNResiduals; k++)
        {
                fx[k] = (nlls->GetResiduals())[k];
        }

        // loop over parameters
        for(i = 0; i < nParameters; i++)
        {
                // transform the parameter and save its old value
                saveparameter = parameters[i];
                transparameter = m_pFilter->Operator(parameters[i]);

                double h = pow(m_dFuncAccuracy,0.333333)*fabs(transparameter);
                
                // if statement warns of tiny stepsizes
                if(h < DBL_EPSILON)
                {
                        cout << "WARNING:Parmeter step small in ImprovedLevenbergMarquardtMinimizer::ComputeDerivativeInformation !" << endl;
                        if(transparameter < DBL_MIN)
                        {
                                cout << "WARNING:Zero parameter detected, truncating . . ." << endl;
                                // pretend like the parameter was really small
                                h = pow(m_dFuncAccuracy,0.333333)*fabs(m_pFilter->Operator(parameters[i] + DBL_EPSILON));
                        }
                }

                        
                // forward step
                parameters[i] = m_pFilter->OperatorInverse(transparameter + h);
                double chiSqPlus = nlls->ObjectiveFunction(parameters);
                for(k = 0; k < m_iNResiduals; k++)
                {
                        fxPlus[k] = (nlls->GetResiduals())[k];
                }
                // backward step
                parameters[i] = m_pFilter->OperatorInverse(transparameter - h);
                double chiSqMinus = nlls->ObjectiveFunction(parameters);
                for(k = 0; k < m_iNResiduals; k++)
                {
                        fxMinus[k] = (nlls->GetResiduals())[k];
                }
                
                parameters[i] = saveparameter;

                // fill in jacobian
                for(j = 0; j < m_iNResiduals; j++)
                {
                        double delta = (fxPlus[j] - fxMinus[j])/(2*h);
                        jacobian[j][i] = delta;
                }
        }
        
        // compute -J^T f and store the ForceMatrix in the process
        for(i = 0; i < nParameters; i++)
        {
                m_pdGrad[i] = 0.0;
                for(int j = 0; j < m_iNResiduals; j++)
                {
                        m_pdGrad[i] += -2.0*jacobian[j][i]*fx[j]; 
                        m_pdForceMatrix[i][j] = -2.0*jacobian[j][i]*fx[j];
                }
        }

        // write the force matrix
        Notify();

        // compute the scaled hessian (it is symmetric)
        for(i = 0; i < nParameters; i++)
        {
                for(k = i; k < nParameters; k++)
                {
                        m_pdAlpha[i][k] = 0.0;
                        for(j = 0; j < m_iNResiduals; j++)
                        {
                                m_pdAlpha[i][k] += 2.0*jacobian[j][i]*jacobian[j][k];
                        }
                        m_pdAlpha[k][i] = m_pdAlpha[i][k];
                }
        }

        // Save the diagonal elements
        this->SaveDiagonal();

        delete [] fxPlus;
        delete [] fx;
        delete [] fxMinus;
        delete [] jacobian[0];
        delete [] jacobian;

        cout << "Hessian calculated . . ." << endl;

        return;
}

// Saves the current diagonal elements of the Hessian in the member
// array m_pdDiagonal

void CLeastSquaresMinimizer::SaveDiagonal()
{
        for(int i = 0; i < nParameters; i++)
        {
                m_pdDiagonal[i] = m_pdAlpha[i][i];
        }
}

---