Hessian matrix estimator
andyfaff
#pragma rtGlobals=3 // Use modern global access method and strict wave access. //See the following link for central difference problems //http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/central-differences/ //the following link says why an EPS of 1e-5 should be used. Basically, if the accuracy for the central difference algorithm is A, then the step size, H, should be A^0.3333. //Thus, if you calculate in DP, the accuracy is on the order of 1e-15. The cube root of 1e-15 is 1e-5. //http://scicomp.stackexchange.com/questions/1875/adaptive-h-for-gradient-estimation CONSTANT EPS = 1e-5 Function Hessian(fun, vec, oi) //this function calculates the second order mixed derivatives of the scalar function fun, around the locus values contained in the vec array // //Function fun(vec, oi) // Wave vec // struct otherInfo, &oi // //do some calculation // return val // End // // vec is an array that specifies the locus around which the derivative matrix is to be estimated // // oi is a structure containing other information to pass to the fun function. You can alter this structure to make sense for your function. // // One of it main uses is in estimating the Hessian matrix for a least squares problem (the inverse of which gives the covariance matrix). // //If you use this to calculate the covariance matrix then you have to follow these steps: // 1) divide Hessian by 2. // 2) Take the matrix inverse // 3) Multiply by chi2/(number of fitted points - number of fitted parameters) Funcref energy, fun Wave vec Struct otherinfo &oi variable ii, jj, nvec, t0, t1, t2, t3 nvec = numpnts(vec) make/n=(nvec, nvec)/d/o M_hessian duplicate/free vec, scale, freevec, normal redimension/d scale, freevec, normal scale = 1/vec normal = 1 for(ii = 0 ; ii < nvec ; ii += 1) for(jj = ii ; jj < nvec ; jj+=1) normal = 1 if(ii == jj) normal[ii] = 1 - 2 * EPS freevec = normal/scale t0 = fun(freevec, oi) normal[ii] = 1 - EPS freevec = normal /scale t1 = fun(freevec, oi) normal[ii] = 1 freevec = normal /scale t2 = fun(freevec, oi) normal[ii] = 1 + EPS freevec = normal /scale t3 = fun(freevec, oi) normal[ii] = 1 + 2 * EPS freevec = normal /scale t4 = fun(freevec, oi) M_Hessian[ii][jj] = (-t0 + 16 * t1 - 30 * t2 + 16 * t3 - t4) / 12/EPS/EPS*scale[ii]*scale[ii] else //f -1,-1 normal[ii] = 1 * (1 - EPS) normal[jj] = 1 * (1 - EPS) freevec = normal / scale t0 = fun(freevec, oi) //f +1,+1 normal[ii] = 1 * (1 + EPS) normal[jj] = 1 * (1 + EPS) freevec = normal / scale t1 = fun(freevec, oi) //f +1,-1 normal[ii] = 1 * (1 + EPS) normal[jj] = 1 * (1 - EPS) freevec = normal / scale t2 = fun(freevec, oi) //f -1,+1 normal[ii] = 1 * (1 - EPS) normal[jj] = 1 * (1 + EPS) freevec = normal / scale t3 = fun(freevec, oi) M_hessian[ii][jj] = (t0 + t1 - t2 - t3)/4/EPS/EPS*scale[ii]*scale[jj] M_hessian[jj][ii] = M_hessian[ii][jj] // the following is a higher order calculation, which is more expensive, but may be worth while. // //f 1,-2 // normal[ii] = 1 + EPS // normal[jj] = 1 - 2 * EPS // freevec = normal / scale // t0 = fun(freevec, oi) // // //f +2, -1 // normal[ii] = 1 + 2 * EPS // normal[jj] = 1- EPS // freevec = normal / scale // t1 = fun(freevec, oi) // // //f -2,+1 // normal[ii] = 1 - 2 * EPS // normal[jj] = 1 + EPS // freevec = normal / scale // t2 = fun(freevec, oi) // // //f -1,+1 // normal[ii] = 1 - EPS // normal[jj] = 1 + 2 * EPS // freevec = normal / scale // t3 = fun(freevec, oi) // // M_hessian[ii][jj] = -63 * (t0 + t1 + t2 + t3) // // //f -1,-2 // normal[ii] = 1 - EPS // normal[jj] = 1 - 2 * EPS // freevec = normal / scale // t0 = fun(freevec, oi) // // //f -2, -1 // normal[ii] = 1 - 2 * EPS // normal[jj] = 1- EPS // freevec = normal / scale // t1 = fun(freevec, oi) // // //f +1,+2 // normal[ii] = 1 + EPS // normal[jj] = 1 + 2* EPS // freevec = normal / scale // t2 = fun(freevec, oi) // // //f+2,+1 // normal[ii] = 1 + 2 * EPS // normal[jj] = 1 + EPS // freevec = normal / scale // t3 = fun(freevec, oi) // M_hessian[ii][jj] += 63 * (t0 + t1 + t2 + t3) // // //f 2,-2 // normal[ii] = 1 + 2 * EPS // normal[jj] = 1 - 2 * EPS // freevec = normal / scale // t0 = fun(freevec, oi) // // //f -2, 2 // normal[ii] = 1 - 2 * EPS // normal[jj] = 1 + 2 * EPS // freevec = normal / scale // t1 = fun(freevec, oi) // // //f -2,-2 // normal[ii] = 1 - 2* EPS // normal[jj] = 1 - 2* EPS // freevec = normal / scale // t2 = fun(freevec, oi) // // //f2,2 // normal[ii] = 1 + 2 * EPS // normal[jj] = 1 +2* EPS // freevec = normal / scale // t3 = fun(freevec, oi) // M_hessian[ii][jj] += 44 * (t0 + t1 - t2 - t3) // // //f -1,-1 // normal[ii] = 1 -1 * EPS // normal[jj] = 1 - 1 * EPS // freevec = normal / scale // t0 = fun(freevec, oi) // // //f 1,1 // normal[ii] = 1 + EPS // normal[jj] = 1 + EPS // freevec = normal / scale // t1 = fun(freevec, oi) // // //f 1,-1 // normal[ii] = 1 + EPS // normal[jj] = 1 - EPS // freevec = normal / scale // t2 = fun(freevec, oi) // // //f-1,1 // normal[ii] = 1 -1 * EPS // normal[jj] = 1 +1* EPS // freevec = normal / scale // t3 = fun(freevec, oi) // M_hessian[ii][jj] += 74 * (t0 + t1 - t2 - t3) // // M_hessian[ii][jj] = M_hessian[ii][jj]/600/EPS/EPS*scale[ii]*scale[jj] // M_hessian[jj][ii] = M_hessian[ii][jj] endif endfor endfor End Function energy(vec, oi) Wave vec Struct otherinfo &oi End Structure otherinfo Endstructure Function rosen(w, oi) Wave w Struct otherinfo &oi return (1-w[0])^2 + 105 * (w[1] - w[0]^2)^2 End Function test(w) Wave w Struct otherinfo pop Hessian(rosen, w, pop) //should give //M_hessian[0][0]= {842,-420} //M_hessian[0][1]= {-420,210} End
Forum
Support
Gallery
Igor Pro 9
Learn More
Igor XOP Toolkit
Learn More
Igor NIDAQ Tools MX
Learn More
Rather than hard-commenting out the higher order, perhaps this change would be useful ...
--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAHuntsville
May 6, 2013 at 05:04 pm - Permalink