1 |
|
#include "NLModel.hpp" |
2 |
|
|
3 |
– |
|
3 |
|
/** |
4 |
|
* calculate gradient using backward finite difference |
5 |
|
* df(Xk)/dXi = (f(Xk) - f(Xk - h*ei)) /h |
33 |
|
|
34 |
|
hi = copysign(h[i], tempX[i]); |
35 |
|
|
37 |
– |
//tempX[i] = x[i] + hi; |
36 |
|
tempX[i] -= hi; |
37 |
|
} |
38 |
|
|
171 |
|
* calculate hessian using finite difference |
172 |
|
* d2f(Xk)/dxidxj = (df(Xk+h*ej)/dXi + df(Xk - h*ej)/dXi) /2h |
173 |
|
* where h is a small positive scalar and ei is the ith unit vector (ith column of the identity Matrix) |
174 |
+ |
* |
175 |
|
*/ |
176 |
< |
SymMatrix NLModel0::FDHessian(vector<double>& h){ |
176 |
> |
SymMatrix NLModel0::FiniteHessian(vector<double>& x, double fx, vector<double>& h){ |
177 |
|
SymMatrix H(ndim); |
178 |
+ |
vector<double> tempX; |
179 |
+ |
vector<double> fi(ndim); |
180 |
+ |
vector<double> hi; |
181 |
+ |
double fij; |
182 |
+ |
|
183 |
+ |
tempX = x; |
184 |
|
|
185 |
+ |
//calculate f(X+ h*ei) |
186 |
|
for(int i = 0; i < ndim; i++){ |
187 |
|
|
188 |
+ |
#ifndef IS_MPI |
189 |
+ |
|
190 |
+ |
hi[i] = copysign(h[i],tempX[i]); |
191 |
+ |
|
192 |
+ |
tempX[i] += hi[i]; |
193 |
+ |
|
194 |
+ |
fi[i] = calcF(tempX); |
195 |
+ |
|
196 |
+ |
tempX[i] -= hi; |
197 |
+ |
|
198 |
+ |
#else |
199 |
+ |
|
200 |
+ |
if(procMappingArray[i] == myRank){ |
201 |
+ |
|
202 |
+ |
hi[i] = copysign(h[i],tempX[i]); |
203 |
+ |
|
204 |
+ |
tempX[i] += hi[i]; |
205 |
+ |
} |
206 |
+ |
|
207 |
+ |
fi[i] = calcF(tempX); |
208 |
+ |
|
209 |
+ |
if(procMappingArray[i] == myRank){ |
210 |
+ |
tempX[i] -= hi[i]; |
211 |
+ |
} |
212 |
+ |
|
213 |
+ |
#endif |
214 |
+ |
} |
215 |
+ |
|
216 |
+ |
#ifndef IS_MPI |
217 |
+ |
|
218 |
+ |
for(int i = 0; i < ndim; i++){ |
219 |
+ |
|
220 |
+ |
tempX[i] += hi[i] * 2; |
221 |
+ |
|
222 |
+ |
fii = calcF(tempX); |
223 |
+ |
|
224 |
+ |
H(i,j) = ((fx - fi[i]) + (fii - fi[i])) / (hi[i]*hi[i]); |
225 |
+ |
|
226 |
+ |
tempX[i] -= hi[i]; |
227 |
+ |
|
228 |
|
for(int j = i + 1; j < ndim; j++){ |
229 |
|
|
230 |
+ |
tempX[j] += hi[j]; |
231 |
+ |
|
232 |
+ |
fij = calcF(tempX); |
233 |
+ |
|
234 |
+ |
H(i,j) = ((f - fi[i]) + (fij - fi[j])) / (hi[i]*hi[j]); |
235 |
+ |
|
236 |
+ |
tempX[j] -= hi[j]; |
237 |
|
} |
238 |
+ |
|
239 |
+ |
tempX[i] -= hi[i]; |
240 |
|
} |
241 |
< |
|
241 |
> |
|
242 |
> |
#else |
243 |
> |
for(int i = 0; i < ndim; i++){ |
244 |
> |
|
245 |
> |
if(procMappingArray[i] == myRank){ |
246 |
> |
tempX[i] += hi[i] * 2; |
247 |
> |
} |
248 |
> |
|
249 |
> |
fii = calcF(tempX); |
250 |
> |
|
251 |
> |
H(i,j) = ((fx - fi[i]) + (fii - fi[i])) / (hi[i]*hi[i]); |
252 |
> |
|
253 |
> |
if(procMappingArray[i] == myRank){ |
254 |
> |
tempX[i] -= hi[i]; |
255 |
> |
} |
256 |
> |
|
257 |
> |
for(int j = i + 1; j < ndim; j++){ |
258 |
> |
|
259 |
> |
if(procMappingArray[j] == myRank){ |
260 |
> |
tempX[j] += hi[j]; |
261 |
> |
} |
262 |
> |
|
263 |
> |
fij = calcF(tempX); |
264 |
> |
H(i,j) = ((f - fi[i]) + (fij - fi[j])) / (hi[i]*hi[j]); |
265 |
> |
|
266 |
> |
if(procMappingArray[j] == myRank){ |
267 |
> |
tempX[j] -= hi[j]; |
268 |
> |
} |
269 |
> |
|
270 |
> |
} |
271 |
> |
|
272 |
> |
if(procMappingArray[i] == myRank){ |
273 |
> |
tempX[i] -= hi[i]; |
274 |
> |
} |
275 |
> |
|
276 |
> |
} |
277 |
> |
|
278 |
> |
#endif |
279 |
> |
|
280 |
> |
return H; |
281 |
|
} |