Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import numpy as np 

2import scipy.sparse as sps 

3from ._numdiff import approx_derivative, group_columns 

4from ._hessian_update_strategy import HessianUpdateStrategy 

5from scipy.sparse.linalg import LinearOperator 

6 

7 

8FD_METHODS = ('2-point', '3-point', 'cs') 

9 

10 

11class ScalarFunction(object): 

12 """Scalar function and its derivatives. 

13 

14 This class defines a scalar function F: R^n->R and methods for 

15 computing or approximating its first and second derivatives. 

16 

17 Notes 

18 ----- 

19 This class implements a memoization logic. There are methods `fun`, 

20 `grad`, hess` and corresponding attributes `f`, `g` and `H`. The following 

21 things should be considered: 

22 

23 1. Use only public methods `fun`, `grad` and `hess`. 

24 2. After one of the methods is called, the corresponding attribute 

25 will be set. However, a subsequent call with a different argument 

26 of *any* of the methods may overwrite the attribute. 

27 """ 

28 def __init__(self, fun, x0, args, grad, hess, finite_diff_rel_step, 

29 finite_diff_bounds, epsilon=None): 

30 if not callable(grad) and grad not in FD_METHODS: 

31 raise ValueError("`grad` must be either callable or one of {}." 

32 .format(FD_METHODS)) 

33 

34 if not (callable(hess) or hess in FD_METHODS 

35 or isinstance(hess, HessianUpdateStrategy)): 

36 raise ValueError("`hess` must be either callable," 

37 "HessianUpdateStrategy or one of {}." 

38 .format(FD_METHODS)) 

39 

40 if grad in FD_METHODS and hess in FD_METHODS: 

41 raise ValueError("Whenever the gradient is estimated via " 

42 "finite-differences, we require the Hessian " 

43 "to be estimated using one of the " 

44 "quasi-Newton strategies.") 

45 

46 self.x = np.atleast_1d(x0).astype(float) 

47 self.n = self.x.size 

48 self.nfev = 0 

49 self.ngev = 0 

50 self.nhev = 0 

51 self.f_updated = False 

52 self.g_updated = False 

53 self.H_updated = False 

54 

55 finite_diff_options = {} 

56 if grad in FD_METHODS: 

57 finite_diff_options["method"] = grad 

58 finite_diff_options["rel_step"] = finite_diff_rel_step 

59 finite_diff_options["abs_step"] = epsilon 

60 finite_diff_options["bounds"] = finite_diff_bounds 

61 if hess in FD_METHODS: 

62 finite_diff_options["method"] = hess 

63 finite_diff_options["rel_step"] = finite_diff_rel_step 

64 finite_diff_options["abs_step"] = epsilon 

65 finite_diff_options["as_linear_operator"] = True 

66 

67 # Function evaluation 

68 def fun_wrapped(x): 

69 self.nfev += 1 

70 return fun(x, *args) 

71 

72 def update_fun(): 

73 self.f = fun_wrapped(self.x) 

74 

75 self._update_fun_impl = update_fun 

76 self._update_fun() 

77 

78 # Gradient evaluation 

79 if callable(grad): 

80 def grad_wrapped(x): 

81 self.ngev += 1 

82 return np.atleast_1d(grad(x, *args)) 

83 

84 def update_grad(): 

85 self.g = grad_wrapped(self.x) 

86 

87 elif grad in FD_METHODS: 

88 def update_grad(): 

89 self._update_fun() 

90 self.ngev += 1 

91 self.g = approx_derivative(fun_wrapped, self.x, f0=self.f, 

92 **finite_diff_options) 

93 

94 self._update_grad_impl = update_grad 

95 self._update_grad() 

96 

97 # Hessian Evaluation 

98 if callable(hess): 

99 self.H = hess(x0, *args) 

100 self.H_updated = True 

101 self.nhev += 1 

102 

103 if sps.issparse(self.H): 

104 def hess_wrapped(x): 

105 self.nhev += 1 

106 return sps.csr_matrix(hess(x, *args)) 

107 self.H = sps.csr_matrix(self.H) 

108 

109 elif isinstance(self.H, LinearOperator): 

110 def hess_wrapped(x): 

111 self.nhev += 1 

112 return hess(x, *args) 

113 

114 else: 

115 def hess_wrapped(x): 

116 self.nhev += 1 

117 return np.atleast_2d(np.asarray(hess(x, *args))) 

118 self.H = np.atleast_2d(np.asarray(self.H)) 

119 

120 def update_hess(): 

121 self.H = hess_wrapped(self.x) 

122 

123 elif hess in FD_METHODS: 

124 def update_hess(): 

125 self._update_grad() 

126 self.H = approx_derivative(grad_wrapped, self.x, f0=self.g, 

127 **finite_diff_options) 

128 return self.H 

129 

130 update_hess() 

131 self.H_updated = True 

132 elif isinstance(hess, HessianUpdateStrategy): 

133 self.H = hess 

134 self.H.initialize(self.n, 'hess') 

135 self.H_updated = True 

136 self.x_prev = None 

137 self.g_prev = None 

138 

139 def update_hess(): 

140 self._update_grad() 

141 self.H.update(self.x - self.x_prev, self.g - self.g_prev) 

142 

143 self._update_hess_impl = update_hess 

144 

145 if isinstance(hess, HessianUpdateStrategy): 

146 def update_x(x): 

147 self._update_grad() 

148 self.x_prev = self.x 

149 self.g_prev = self.g 

150 

151 self.x = np.atleast_1d(x).astype(float) 

152 self.f_updated = False 

153 self.g_updated = False 

154 self.H_updated = False 

155 self._update_hess() 

156 else: 

157 def update_x(x): 

158 self.x = np.atleast_1d(x).astype(float) 

159 self.f_updated = False 

160 self.g_updated = False 

161 self.H_updated = False 

162 self._update_x_impl = update_x 

163 

164 def _update_fun(self): 

165 if not self.f_updated: 

166 self._update_fun_impl() 

167 self.f_updated = True 

168 

169 def _update_grad(self): 

170 if not self.g_updated: 

171 self._update_grad_impl() 

172 self.g_updated = True 

173 

174 def _update_hess(self): 

175 if not self.H_updated: 

176 self._update_hess_impl() 

177 self.H_updated = True 

178 

179 def fun(self, x): 

180 if not np.array_equal(x, self.x): 

181 self._update_x_impl(x) 

182 self._update_fun() 

183 return self.f 

184 

185 def grad(self, x): 

186 if not np.array_equal(x, self.x): 

187 self._update_x_impl(x) 

188 self._update_grad() 

189 return self.g 

190 

191 def hess(self, x): 

192 if not np.array_equal(x, self.x): 

193 self._update_x_impl(x) 

194 self._update_hess() 

195 return self.H 

196 

197 def fun_and_grad(self, x): 

198 if not np.array_equal(x, self.x): 

199 self._update_x_impl(x) 

200 self._update_fun() 

201 self._update_grad() 

202 return self.f, self.g 

203 

204 

205class VectorFunction(object): 

206 """Vector function and its derivatives. 

207 

208 This class defines a vector function F: R^n->R^m and methods for 

209 computing or approximating its first and second derivatives. 

210 

211 Notes 

212 ----- 

213 This class implements a memoization logic. There are methods `fun`, 

214 `jac`, hess` and corresponding attributes `f`, `J` and `H`. The following 

215 things should be considered: 

216 

217 1. Use only public methods `fun`, `jac` and `hess`. 

218 2. After one of the methods is called, the corresponding attribute 

219 will be set. However, a subsequent call with a different argument 

220 of *any* of the methods may overwrite the attribute. 

221 """ 

222 def __init__(self, fun, x0, jac, hess, 

223 finite_diff_rel_step, finite_diff_jac_sparsity, 

224 finite_diff_bounds, sparse_jacobian): 

225 if not callable(jac) and jac not in FD_METHODS: 

226 raise ValueError("`jac` must be either callable or one of {}." 

227 .format(FD_METHODS)) 

228 

229 if not (callable(hess) or hess in FD_METHODS 

230 or isinstance(hess, HessianUpdateStrategy)): 

231 raise ValueError("`hess` must be either callable," 

232 "HessianUpdateStrategy or one of {}." 

233 .format(FD_METHODS)) 

234 

235 if jac in FD_METHODS and hess in FD_METHODS: 

236 raise ValueError("Whenever the Jacobian is estimated via " 

237 "finite-differences, we require the Hessian to " 

238 "be estimated using one of the quasi-Newton " 

239 "strategies.") 

240 

241 self.x = np.atleast_1d(x0).astype(float) 

242 self.n = self.x.size 

243 self.nfev = 0 

244 self.njev = 0 

245 self.nhev = 0 

246 self.f_updated = False 

247 self.J_updated = False 

248 self.H_updated = False 

249 

250 finite_diff_options = {} 

251 if jac in FD_METHODS: 

252 finite_diff_options["method"] = jac 

253 finite_diff_options["rel_step"] = finite_diff_rel_step 

254 if finite_diff_jac_sparsity is not None: 

255 sparsity_groups = group_columns(finite_diff_jac_sparsity) 

256 finite_diff_options["sparsity"] = (finite_diff_jac_sparsity, 

257 sparsity_groups) 

258 finite_diff_options["bounds"] = finite_diff_bounds 

259 self.x_diff = np.copy(self.x) 

260 if hess in FD_METHODS: 

261 finite_diff_options["method"] = hess 

262 finite_diff_options["rel_step"] = finite_diff_rel_step 

263 finite_diff_options["as_linear_operator"] = True 

264 self.x_diff = np.copy(self.x) 

265 if jac in FD_METHODS and hess in FD_METHODS: 

266 raise ValueError("Whenever the Jacobian is estimated via " 

267 "finite-differences, we require the Hessian to " 

268 "be estimated using one of the quasi-Newton " 

269 "strategies.") 

270 

271 # Function evaluation 

272 def fun_wrapped(x): 

273 self.nfev += 1 

274 return np.atleast_1d(fun(x)) 

275 

276 def update_fun(): 

277 self.f = fun_wrapped(self.x) 

278 

279 self._update_fun_impl = update_fun 

280 update_fun() 

281 

282 self.v = np.zeros_like(self.f) 

283 self.m = self.v.size 

284 

285 # Jacobian Evaluation 

286 if callable(jac): 

287 self.J = jac(self.x) 

288 self.J_updated = True 

289 self.njev += 1 

290 

291 if (sparse_jacobian or 

292 sparse_jacobian is None and sps.issparse(self.J)): 

293 def jac_wrapped(x): 

294 self.njev += 1 

295 return sps.csr_matrix(jac(x)) 

296 self.J = sps.csr_matrix(self.J) 

297 self.sparse_jacobian = True 

298 

299 elif sps.issparse(self.J): 

300 def jac_wrapped(x): 

301 self.njev += 1 

302 return jac(x).toarray() 

303 self.J = self.J.toarray() 

304 self.sparse_jacobian = False 

305 

306 else: 

307 def jac_wrapped(x): 

308 self.njev += 1 

309 return np.atleast_2d(jac(x)) 

310 self.J = np.atleast_2d(self.J) 

311 self.sparse_jacobian = False 

312 

313 def update_jac(): 

314 self.J = jac_wrapped(self.x) 

315 

316 elif jac in FD_METHODS: 

317 self.J = approx_derivative(fun_wrapped, self.x, f0=self.f, 

318 **finite_diff_options) 

319 self.J_updated = True 

320 

321 if (sparse_jacobian or 

322 sparse_jacobian is None and sps.issparse(self.J)): 

323 def update_jac(): 

324 self._update_fun() 

325 self.J = sps.csr_matrix( 

326 approx_derivative(fun_wrapped, self.x, f0=self.f, 

327 **finite_diff_options)) 

328 self.J = sps.csr_matrix(self.J) 

329 self.sparse_jacobian = True 

330 

331 elif sps.issparse(self.J): 

332 def update_jac(): 

333 self._update_fun() 

334 self.J = approx_derivative(fun_wrapped, self.x, f0=self.f, 

335 **finite_diff_options).toarray() 

336 self.J = self.J.toarray() 

337 self.sparse_jacobian = False 

338 

339 else: 

340 def update_jac(): 

341 self._update_fun() 

342 self.J = np.atleast_2d( 

343 approx_derivative(fun_wrapped, self.x, f0=self.f, 

344 **finite_diff_options)) 

345 self.J = np.atleast_2d(self.J) 

346 self.sparse_jacobian = False 

347 

348 self._update_jac_impl = update_jac 

349 

350 # Define Hessian 

351 if callable(hess): 

352 self.H = hess(self.x, self.v) 

353 self.H_updated = True 

354 self.nhev += 1 

355 

356 if sps.issparse(self.H): 

357 def hess_wrapped(x, v): 

358 self.nhev += 1 

359 return sps.csr_matrix(hess(x, v)) 

360 self.H = sps.csr_matrix(self.H) 

361 

362 elif isinstance(self.H, LinearOperator): 

363 def hess_wrapped(x, v): 

364 self.nhev += 1 

365 return hess(x, v) 

366 

367 else: 

368 def hess_wrapped(x, v): 

369 self.nhev += 1 

370 return np.atleast_2d(np.asarray(hess(x, v))) 

371 self.H = np.atleast_2d(np.asarray(self.H)) 

372 

373 def update_hess(): 

374 self.H = hess_wrapped(self.x, self.v) 

375 elif hess in FD_METHODS: 

376 def jac_dot_v(x, v): 

377 return jac_wrapped(x).T.dot(v) 

378 

379 def update_hess(): 

380 self._update_jac() 

381 self.H = approx_derivative(jac_dot_v, self.x, 

382 f0=self.J.T.dot(self.v), 

383 args=(self.v,), 

384 **finite_diff_options) 

385 update_hess() 

386 self.H_updated = True 

387 elif isinstance(hess, HessianUpdateStrategy): 

388 self.H = hess 

389 self.H.initialize(self.n, 'hess') 

390 self.H_updated = True 

391 self.x_prev = None 

392 self.J_prev = None 

393 

394 def update_hess(): 

395 self._update_jac() 

396 # When v is updated before x was updated, then x_prev and 

397 # J_prev are None and we need this check. 

398 if self.x_prev is not None and self.J_prev is not None: 

399 delta_x = self.x - self.x_prev 

400 delta_g = self.J.T.dot(self.v) - self.J_prev.T.dot(self.v) 

401 self.H.update(delta_x, delta_g) 

402 

403 self._update_hess_impl = update_hess 

404 

405 if isinstance(hess, HessianUpdateStrategy): 

406 def update_x(x): 

407 self._update_jac() 

408 self.x_prev = self.x 

409 self.J_prev = self.J 

410 self.x = np.atleast_1d(x).astype(float) 

411 self.f_updated = False 

412 self.J_updated = False 

413 self.H_updated = False 

414 self._update_hess() 

415 else: 

416 def update_x(x): 

417 self.x = np.atleast_1d(x).astype(float) 

418 self.f_updated = False 

419 self.J_updated = False 

420 self.H_updated = False 

421 

422 self._update_x_impl = update_x 

423 

424 def _update_v(self, v): 

425 if not np.array_equal(v, self.v): 

426 self.v = v 

427 self.H_updated = False 

428 

429 def _update_x(self, x): 

430 if not np.array_equal(x, self.x): 

431 self._update_x_impl(x) 

432 

433 def _update_fun(self): 

434 if not self.f_updated: 

435 self._update_fun_impl() 

436 self.f_updated = True 

437 

438 def _update_jac(self): 

439 if not self.J_updated: 

440 self._update_jac_impl() 

441 self.J_updated = True 

442 

443 def _update_hess(self): 

444 if not self.H_updated: 

445 self._update_hess_impl() 

446 self.H_updated = True 

447 

448 def fun(self, x): 

449 self._update_x(x) 

450 self._update_fun() 

451 return self.f 

452 

453 def jac(self, x): 

454 self._update_x(x) 

455 self._update_jac() 

456 return self.J 

457 

458 def hess(self, x, v): 

459 # v should be updated before x. 

460 self._update_v(v) 

461 self._update_x(x) 

462 self._update_hess() 

463 return self.H 

464 

465 

466class LinearVectorFunction(object): 

467 """Linear vector function and its derivatives. 

468 

469 Defines a linear function F = A x, where x is N-D vector and 

470 A is m-by-n matrix. The Jacobian is constant and equals to A. The Hessian 

471 is identically zero and it is returned as a csr matrix. 

472 """ 

473 def __init__(self, A, x0, sparse_jacobian): 

474 if sparse_jacobian or sparse_jacobian is None and sps.issparse(A): 

475 self.J = sps.csr_matrix(A) 

476 self.sparse_jacobian = True 

477 elif sps.issparse(A): 

478 self.J = A.toarray() 

479 self.sparse_jacobian = False 

480 else: 

481 # np.asarray makes sure A is ndarray and not matrix 

482 self.J = np.atleast_2d(np.asarray(A)) 

483 self.sparse_jacobian = False 

484 

485 self.m, self.n = self.J.shape 

486 

487 self.x = np.atleast_1d(x0).astype(float) 

488 self.f = self.J.dot(self.x) 

489 self.f_updated = True 

490 

491 self.v = np.zeros(self.m, dtype=float) 

492 self.H = sps.csr_matrix((self.n, self.n)) 

493 

494 def _update_x(self, x): 

495 if not np.array_equal(x, self.x): 

496 self.x = np.atleast_1d(x).astype(float) 

497 self.f_updated = False 

498 

499 def fun(self, x): 

500 self._update_x(x) 

501 if not self.f_updated: 

502 self.f = self.J.dot(x) 

503 self.f_updated = True 

504 return self.f 

505 

506 def jac(self, x): 

507 self._update_x(x) 

508 return self.J 

509 

510 def hess(self, x, v): 

511 self._update_x(x) 

512 self.v = v 

513 return self.H 

514 

515 

516class IdentityVectorFunction(LinearVectorFunction): 

517 """Identity vector function and its derivatives. 

518 

519 The Jacobian is the identity matrix, returned as a dense array when 

520 `sparse_jacobian=False` and as a csr matrix otherwise. The Hessian is 

521 identically zero and it is returned as a csr matrix. 

522 """ 

523 def __init__(self, x0, sparse_jacobian): 

524 n = len(x0) 

525 if sparse_jacobian or sparse_jacobian is None: 

526 A = sps.eye(n, format='csr') 

527 sparse_jacobian = True 

528 else: 

529 A = np.eye(n) 

530 sparse_jacobian = False 

531 super(IdentityVectorFunction, self).__init__(A, x0, sparse_jacobian)