1 |
|
#include "NLModel.hpp" |
2 |
+ |
#include <math.h> |
3 |
|
|
3 |
– |
|
4 |
|
/** |
5 |
|
* calculate gradient using backward finite difference |
6 |
|
* df(Xk)/dXi = (f(Xk) - f(Xk - h*ei)) /h |
22 |
|
//tempX[i] = x[i] + hi; |
23 |
|
tempX[i] -= hi; |
24 |
|
|
25 |
< |
fminus = (*objfunc)(tempX); |
25 |
> |
fminus = calcF(tempX); |
26 |
|
|
27 |
|
partialGrad[i] = (fx - fminus) / hi; |
28 |
|
|
34 |
|
|
35 |
|
hi = copysign(h[i], tempX[i]); |
36 |
|
|
37 |
– |
//tempX[i] = x[i] + hi; |
37 |
|
tempX[i] -= hi; |
38 |
|
} |
39 |
|
|
40 |
< |
fminus = (*objfunc)(tempX); |
40 |
> |
fminus = calcF(tempX); |
41 |
|
|
42 |
|
if(procMappingArray[i] == myRank){ |
43 |
|
partialGrad[i] = (fx - fminus) / hi; |
73 |
|
//tempX[i] = x[i] + hi; |
74 |
|
tempX[i] += hi; |
75 |
|
|
76 |
< |
fplus = (*objfunc)(tempX); |
76 |
> |
fplus = calcF(tempX); |
77 |
|
|
78 |
|
partialGrad[i] = (fplus - fx) / hi; |
79 |
|
|
89 |
|
tempX[i] += hi; |
90 |
|
} |
91 |
|
|
92 |
< |
fminus = (*objfunc)(tempX); |
92 |
> |
fminus = calcF(tempX); |
93 |
|
|
94 |
|
if(procMappingArray[i] == myRank){ |
95 |
|
partialGrad[i] = (fx - fminus) / hi; |
125 |
|
//tempX[i] = x[i] + hi |
126 |
|
tempX[i] += hi; |
127 |
|
|
128 |
< |
fplus = (*objfunc)(tempX); |
128 |
> |
fplus = calcF(tempX); |
129 |
|
|
130 |
|
//tempX[i] = x[i] -hi |
131 |
|
tempX[i] -= 2*hi; |
132 |
< |
fminus = (*objfunc)(tempX); |
132 |
> |
fminus = calcF(tempX); |
133 |
|
|
134 |
|
partialGrad[i] = (fplus + fminus) / (2*hi); |
135 |
|
|
145 |
|
tempX[i] += hi; |
146 |
|
} |
147 |
|
|
148 |
< |
fplus = (*objfunc)(tempX); |
148 |
> |
fplus = calcF(tempX); |
149 |
|
|
150 |
|
if(procMappingArray[i] == myRank){ |
151 |
|
partialGrad[i] = (fx - fminus) / hi; |
172 |
|
* calculate hessian using finite difference |
173 |
|
* d2f(Xk)/dxidxj = (df(Xk+h*ej)/dXi + df(Xk - h*ej)/dXi) /2h |
174 |
|
* where h is a small positive scalar and ei is the ith unit vector (ith column of the identity Matrix) |
175 |
+ |
* |
176 |
|
*/ |
177 |
< |
SymMatrix NLModel0::FDHessian(vector<double>& h){ |
177 |
> |
SymMatrix NLModel0::FiniteHessian(vector<double>& x, double fx, vector<double>& h){ |
178 |
|
SymMatrix H(ndim); |
179 |
+ |
vector<double> tempX; |
180 |
+ |
vector<double> fi(ndim); |
181 |
+ |
vector<double> hi; |
182 |
+ |
double fii, fij; |
183 |
+ |
|
184 |
+ |
tempX = x; |
185 |
|
|
186 |
+ |
//calculate f(X+ h*ei) |
187 |
|
for(int i = 0; i < ndim; i++){ |
188 |
|
|
189 |
+ |
#ifndef IS_MPI |
190 |
+ |
|
191 |
+ |
hi[i] = copysign(h[i],tempX[i]); |
192 |
+ |
|
193 |
+ |
tempX[i] += hi[i]; |
194 |
+ |
|
195 |
+ |
fi[i] = calcF(tempX); |
196 |
+ |
|
197 |
+ |
tempX[i] -= hi[i]; |
198 |
+ |
|
199 |
+ |
#else |
200 |
+ |
|
201 |
+ |
if(procMappingArray[i] == myRank){ |
202 |
+ |
|
203 |
+ |
hi[i] = copysign(h[i],tempX[i]); |
204 |
+ |
|
205 |
+ |
tempX[i] += hi[i]; |
206 |
+ |
} |
207 |
+ |
|
208 |
+ |
fi[i] = calcF(tempX); |
209 |
+ |
|
210 |
+ |
if(procMappingArray[i] == myRank){ |
211 |
+ |
tempX[i] -= hi[i]; |
212 |
+ |
} |
213 |
+ |
|
214 |
+ |
#endif |
215 |
+ |
} |
216 |
+ |
|
217 |
+ |
#ifndef IS_MPI |
218 |
+ |
|
219 |
+ |
for(int i = 0; i < ndim; i++){ |
220 |
+ |
|
221 |
+ |
tempX[i] += hi[i] * 2; |
222 |
+ |
|
223 |
+ |
fii = calcF(tempX); |
224 |
+ |
|
225 |
+ |
H(i, i) = ((fx - fi[i]) + (fii - fi[i])) / (hi[i]*hi[i]); |
226 |
+ |
|
227 |
+ |
tempX[i] -= hi[i]; |
228 |
+ |
|
229 |
|
for(int j = i + 1; j < ndim; j++){ |
230 |
|
|
231 |
+ |
tempX[j] += hi[j]; |
232 |
+ |
|
233 |
+ |
fij = calcF(tempX); |
234 |
+ |
|
235 |
+ |
H(i,j) = ((fx - fi[i]) + (fij - fi[j])) / (hi[i]*hi[j]); |
236 |
+ |
|
237 |
+ |
tempX[j] -= hi[j]; |
238 |
|
} |
239 |
+ |
|
240 |
+ |
tempX[i] -= hi[i]; |
241 |
|
} |
242 |
< |
|
242 |
> |
|
243 |
> |
#else |
244 |
> |
for(int i = 0; i < ndim; i++){ |
245 |
> |
|
246 |
> |
if(procMappingArray[i] == myRank){ |
247 |
> |
tempX[i] += hi[i] * 2; |
248 |
> |
} |
249 |
> |
|
250 |
> |
fii = calcF(tempX); |
251 |
> |
|
252 |
> |
H(i,i) = ((fx - fi[i]) + (fii - fi[i])) / (hi[i]*hi[i]); |
253 |
> |
|
254 |
> |
if(procMappingArray[i] == myRank){ |
255 |
> |
tempX[i] -= hi[i]; |
256 |
> |
} |
257 |
> |
|
258 |
> |
for(int j = i + 1; j < ndim; j++){ |
259 |
> |
|
260 |
> |
if(procMappingArray[j] == myRank){ |
261 |
> |
tempX[j] += hi[j]; |
262 |
> |
} |
263 |
> |
|
264 |
> |
fij = calcF(tempX); |
265 |
> |
H(i,j) = ((fx - fi[i]) + (fij - fi[j])) / (hi[i]*hi[j]); |
266 |
> |
|
267 |
> |
if(procMappingArray[j] == myRank){ |
268 |
> |
tempX[j] -= hi[j]; |
269 |
> |
} |
270 |
> |
|
271 |
> |
} |
272 |
> |
|
273 |
> |
if(procMappingArray[i] == myRank){ |
274 |
> |
tempX[i] -= hi[i]; |
275 |
> |
} |
276 |
> |
|
277 |
> |
} |
278 |
> |
|
279 |
> |
#endif |
280 |
> |
|
281 |
> |
return H; |
282 |
|
} |