Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Compound Members   File Members

# LeastSquaresMinimizer.cpp

Go to the documentation of this file.
```00001 // LeastSquaresMinimizer.cpp: implementation of the LeastSquaresMinimizer class.
00002 //
00004
00005 #include "LeastSquaresMinimizer.h"
00006
00007
00009 // Construction/Destruction
00011
00012 CLeastSquaresMinimizer::CLeastSquaresMinimizer()
00013 {
00014
00015 }
00016
00017 CLeastSquaresMinimizer::~CLeastSquaresMinimizer()
00018 {
00019
00020 }
00021
00022
00024 // Compute approximate Hessian for use in Least Squares algs.
00025 // The hessian is symmetric so we only need compute 1/2 of it then
00026 // copy the rest.  This routine assumes analytic second derivatives
00027 // are not available, so finite differencing with a 3-pt central
00028 // formula:
00029 //                      f'(x) ~ (f(x+h) - f(x-h))/(2*h)
00030 // is employed to obtain derivatives.
00032
00033 void CLeastSquaresMinimizer::ComputeDerivativeInformation(double *parameters, NLLSMinimizable *nlls)
00034 {
00035         int i,j,k,nR;
00036         double transparameter = 0.0;
00037         double saveparameter = 0.0;
00038         double *fxPlus = new double[m_iNResiduals];
00039         double *fx = new double[m_iNResiduals];
00040         double *fxMinus = new double[m_iNResiduals];
00041         double **jacobian = new double *[m_iNResiduals];
00042         jacobian[0] = new double[m_iNResiduals*nParameters];
00043         for(nR = 1; nR < m_iNResiduals; nR++)
00044         {
00045                 jacobian[nR] = &(jacobian[0][nR*nParameters]);
00046         }
00047
00048         double chiSq = nlls->ObjectiveFunction(parameters);
00049         for(k = 0; k < m_iNResiduals; k++)
00050         {
00051                 fx[k] = (nlls->GetResiduals())[k];
00052         }
00053
00054         // loop over parameters
00055         for(i = 0; i < nParameters; i++)
00056         {
00057                 // transform the parameter and save its old value
00058                 saveparameter = parameters[i];
00059                 transparameter = m_pFilter->Operator(parameters[i]);
00060
00061                 double h = pow(m_dFuncAccuracy,0.333333)*fabs(transparameter);
00062
00063                 // if statement warns of tiny stepsizes
00064                 if(h < DBL_EPSILON)
00065                 {
00066                         cout << "WARNING:Parmeter step small in ImprovedLevenbergMarquardtMinimizer::ComputeDerivativeInformation !" << endl;
00067                         if(transparameter < DBL_MIN)
00068                         {
00069                                 cout << "WARNING:Zero parameter detected, truncating . . ." << endl;
00070                                 // pretend like the parameter was really small
00071                                 h = pow(m_dFuncAccuracy,0.333333)*fabs(m_pFilter->Operator(parameters[i] + DBL_EPSILON));
00072                         }
00073                 }
00074
00075
00076                 // forward step
00077                 parameters[i] = m_pFilter->OperatorInverse(transparameter + h);
00078                 double chiSqPlus = nlls->ObjectiveFunction(parameters);
00079                 for(k = 0; k < m_iNResiduals; k++)
00080                 {
00081                         fxPlus[k] = (nlls->GetResiduals())[k];
00082                 }
00083                 // backward step
00084                 parameters[i] = m_pFilter->OperatorInverse(transparameter - h);
00085                 double chiSqMinus = nlls->ObjectiveFunction(parameters);
00086                 for(k = 0; k < m_iNResiduals; k++)
00087                 {
00088                         fxMinus[k] = (nlls->GetResiduals())[k];
00089                 }
00090
00091                 parameters[i] = saveparameter;
00092
00093                 // fill in jacobian
00094                 for(j = 0; j < m_iNResiduals; j++)
00095                 {
00096                         double delta = (fxPlus[j] - fxMinus[j])/(2*h);
00097                         jacobian[j][i] = delta;
00098                 }
00099         }
00100
00101         // compute -J^T f and store the ForceMatrix in the process
00102         for(i = 0; i < nParameters; i++)
00103         {
00105                 for(int j = 0; j < m_iNResiduals; j++)
00106                 {
00108                         m_pdForceMatrix[i][j] = -2.0*jacobian[j][i]*fx[j];
00109                 }
00110         }
00111
00112         // write the force matrix
00113         Notify();
00114
00115         // compute the scaled hessian (it is symmetric)
00116         for(i = 0; i < nParameters; i++)
00117         {
00118                 for(k = i; k < nParameters; k++)
00119                 {
00120                         m_pdAlpha[i][k] = 0.0;
00121                         for(j = 0; j < m_iNResiduals; j++)
00122                         {
00123                                 m_pdAlpha[i][k] += 2.0*jacobian[j][i]*jacobian[j][k];
00124                         }
00125                         m_pdAlpha[k][i] = m_pdAlpha[i][k];
00126                 }
00127         }
00128
00129         // Save the diagonal elements
00130         this->SaveDiagonal();
00131
00132         delete [] fxPlus;
00133         delete [] fx;
00134         delete [] fxMinus;
00135         delete [] jacobian[0];
00136         delete [] jacobian;
00137
00138         cout << "Hessian calculated . . ." << endl;
00139
00140         return;
00141 }
00142
00144 // Saves the current diagonal elements of the Hessian in the member
00145 // array m_pdDiagonal
00147
00148 void CLeastSquaresMinimizer::SaveDiagonal()
00149 {
00150         for(int i = 0; i < nParameters; i++)
00151         {
00152                 m_pdDiagonal[i] = m_pdAlpha[i][i];
00153         }
00154 }
```

Generated on Mon Nov 3 09:37:56 2003 by 1.2.14 written by Dimitri van Heesch, © 1997-2002