Coverage for /usr/lib/python3/dist-packages/scipy/optimize/_lbfgsb

1"""

2Functions

3---------

4.. autosummary::

5 :toctree: generated/

7 fmin_l_bfgs_b

9"""

11## License for the Python wrapper

12## ==============================

16## Permission is hereby granted, free of charge, to any person obtaining a

17## copy of this software and associated documentation files (the "Software"),

18## to deal in the Software without restriction, including without limitation

19## the rights to use, copy, modify, merge, publish, distribute, sublicense,

20## and/or sell copies of the Software, and to permit persons to whom the

21## Software is furnished to do so, subject to the following conditions:

23## The above copyright notice and this permission notice shall be included in

24## all copies or substantial portions of the Software.

26## THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR

27## IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

28## FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE

29## AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER

30## LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING

31## FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER

32## DEALINGS IN THE SOFTWARE.

34## Modifications by Travis Oliphant and Enthought, Inc. for inclusion in SciPy

36import numpy as np

37from numpy import array, asarray, float64, zeros

38from . import _lbfgsb

39from ._optimize import (MemoizeJac, OptimizeResult, _call_callback_maybe_halt,

40 _wrap_callback, _check_unknown_options,

41 _prepare_scalar_function)

42from ._constraints import old_bound_to_new

44from scipy.sparse.linalg import LinearOperator

46__all__ = ['fmin_l_bfgs_b', 'LbfgsInvHessProduct']

49def fmin_l_bfgs_b(func, x0, fprime=None, args=(),

50 approx_grad=0,

51 bounds=None, m=10, factr=1e7, pgtol=1e-5,

52 epsilon=1e-8,

53 iprint=-1, maxfun=15000, maxiter=15000, disp=None,

54 callback=None, maxls=20):

55 """

56 Minimize a function func using the L-BFGS-B algorithm.

58 Parameters

59 ----------

60 func : callable f(x,*args)

61 Function to minimize.

62 x0 : ndarray

63 Initial guess.

64 fprime : callable fprime(x,*args), optional

65 The gradient of `func`. If None, then `func` returns the function

66 value and the gradient (``f, g = func(x, *args)``), unless

67 `approx_grad` is True in which case `func` returns only ``f``.

68 args : sequence, optional

69 Arguments to pass to `func` and `fprime`.

70 approx_grad : bool, optional

71 Whether to approximate the gradient numerically (in which case

72 `func` returns only the function value).

73 bounds : list, optional

74 ``(min, max)`` pairs for each element in ``x``, defining

75 the bounds on that parameter. Use None or +-inf for one of ``min`` or

76 ``max`` when there is no bound in that direction.

77 m : int, optional

78 The maximum number of variable metric corrections

79 used to define the limited memory matrix. (The limited memory BFGS

80 method does not store the full hessian but uses this many terms in an

81 approximation to it.)

82 factr : float, optional

83 The iteration stops when

84 ``(f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr * eps``,

85 where ``eps`` is the machine precision, which is automatically

86 generated by the code. Typical values for `factr` are: 1e12 for

87 low accuracy; 1e7 for moderate accuracy; 10.0 for extremely

88 high accuracy. See Notes for relationship to `ftol`, which is exposed

89 (instead of `factr`) by the `scipy.optimize.minimize` interface to

90 L-BFGS-B.

91 pgtol : float, optional

92 The iteration will stop when

93 ``max{|proj g_i | i = 1, ..., n} <= pgtol``

94 where ``pg_i`` is the i-th component of the projected gradient.

95 epsilon : float, optional

96 Step size used when `approx_grad` is True, for numerically

97 calculating the gradient

98 iprint : int, optional

99 Controls the frequency of output. ``iprint < 0`` means no output;

100 ``iprint = 0`` print only one line at the last iteration;

101 ``0 < iprint < 99`` print also f and ``|proj g|`` every iprint iterations;

102 ``iprint = 99`` print details of every iteration except n-vectors;

103 ``iprint = 100`` print also the changes of active set and final x;

104 ``iprint > 100`` print details of every iteration including x and g.

105 disp : int, optional

106 If zero, then no output. If a positive number, then this over-rides

107 `iprint` (i.e., `iprint` gets the value of `disp`).

108 maxfun : int, optional

109 Maximum number of function evaluations. Note that this function

110 may violate the limit because of evaluating gradients by numerical

111 differentiation.

112 maxiter : int, optional

113 Maximum number of iterations.

114 callback : callable, optional

115 Called after each iteration, as ``callback(xk)``, where ``xk`` is the

116 current parameter vector.

117 maxls : int, optional

118 Maximum number of line search steps (per iteration). Default is 20.

119

120 Returns

121 -------

122 x : array_like

123 Estimated position of the minimum.

124 f : float

125 Value of `func` at the minimum.

126 d : dict

127 Information dictionary.

128

129 * d['warnflag'] is

130

131 - 0 if converged,

132 - 1 if too many function evaluations or too many iterations,

133 - 2 if stopped for another reason, given in d['task']

134

135 * d['grad'] is the gradient at the minimum (should be 0 ish)

136 * d['funcalls'] is the number of function calls made.

137 * d['nit'] is the number of iterations.

138

139 See also

140 --------

141 minimize: Interface to minimization algorithms for multivariate

142 functions. See the 'L-BFGS-B' `method` in particular. Note that the

143 `ftol` option is made available via that interface, while `factr` is

144 provided via this interface, where `factr` is the factor multiplying

145 the default machine floating-point precision to arrive at `ftol`:

146 ``ftol = factr * numpy.finfo(float).eps``.

147

148 Notes

149 -----

150 License of L-BFGS-B (FORTRAN code):

151

152 The version included here (in fortran code) is 3.0

153 (released April 25, 2011). It was written by Ciyou Zhu, Richard Byrd,

154 and Jorge Nocedal <nocedal@ece.nwu.edu>. It carries the following

155 condition for use:

156

157 This software is freely available, but we expect that all publications

158 describing work using this software, or all commercial products using it,

159 quote at least one of the references given below. This software is released

160 under the BSD License.

161

162 References

163 ----------

164 * R. H. Byrd, P. Lu and J. Nocedal. A Limited Memory Algorithm for Bound

165 Constrained Optimization, (1995), SIAM Journal on Scientific and

166 Statistical Computing, 16, 5, pp. 1190-1208.

167 * C. Zhu, R. H. Byrd and J. Nocedal. L-BFGS-B: Algorithm 778: L-BFGS-B,

168 FORTRAN routines for large scale bound constrained optimization (1997),

169 ACM Transactions on Mathematical Software, 23, 4, pp. 550 - 560.

170 * J.L. Morales and J. Nocedal. L-BFGS-B: Remark on Algorithm 778: L-BFGS-B,

171 FORTRAN routines for large scale bound constrained optimization (2011),

172 ACM Transactions on Mathematical Software, 38, 1.

173

174 """

175 # handle fprime/approx_grad

176 if approx_grad:

177 fun = func

178 jac = None

179 elif fprime is None:

180 fun = MemoizeJac(func)

181 jac = fun.derivative

182 else:

183 fun = func

184 jac = fprime

185

186 # build options

187 callback = _wrap_callback(callback)

188 opts = {'disp': disp,

189 'iprint': iprint,

190 'maxcor': m,

191 'ftol': factr * np.finfo(float).eps,

192 'gtol': pgtol,

193 'eps': epsilon,

194 'maxfun': maxfun,

195 'maxiter': maxiter,

196 'callback': callback,

197 'maxls': maxls}

198

199 res = _minimize_lbfgsb(fun, x0, args=args, jac=jac, bounds=bounds,

200 **opts)

201 d = {'grad': res['jac'],

202 'task': res['message'],

203 'funcalls': res['nfev'],

204 'nit': res['nit'],

205 'warnflag': res['status']}

206 f = res['fun']

207 x = res['x']

208

209 return x, f, d

210

211

212def _minimize_lbfgsb(fun, x0, args=(), jac=None, bounds=None,

213 disp=None, maxcor=10, ftol=2.2204460492503131e-09,

214 gtol=1e-5, eps=1e-8, maxfun=15000, maxiter=15000,

215 iprint=-1, callback=None, maxls=20,

216 finite_diff_rel_step=None, **unknown_options):

217 """

218 Minimize a scalar function of one or more variables using the L-BFGS-B

219 algorithm.

220

221 Options

222 -------

223 disp : None or int

224 If `disp is None` (the default), then the supplied version of `iprint`

225 is used. If `disp is not None`, then it overrides the supplied version

226 of `iprint` with the behaviour you outlined.

227 maxcor : int

228 The maximum number of variable metric corrections used to

229 define the limited memory matrix. (The limited memory BFGS

230 method does not store the full hessian but uses this many terms

231 in an approximation to it.)

232 ftol : float

233 The iteration stops when ``(f^k -

234 f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= ftol``.

235 gtol : float

236 The iteration will stop when ``max{|proj g_i | i = 1, ..., n}

237 <= gtol`` where ``pg_i`` is the i-th component of the

238 projected gradient.

239 eps : float or ndarray

240 If `jac is None` the absolute step size used for numerical

241 approximation of the jacobian via forward differences.

242 maxfun : int

243 Maximum number of function evaluations. Note that this function

244 may violate the limit because of evaluating gradients by numerical

245 differentiation.

246 maxiter : int

247 Maximum number of iterations.

248 iprint : int, optional

249 Controls the frequency of output. ``iprint < 0`` means no output;

250 ``iprint = 0`` print only one line at the last iteration;

251 ``0 < iprint < 99`` print also f and ``|proj g|`` every iprint iterations;

252 ``iprint = 99`` print details of every iteration except n-vectors;

253 ``iprint = 100`` print also the changes of active set and final x;

254 ``iprint > 100`` print details of every iteration including x and g.

255 maxls : int, optional

256 Maximum number of line search steps (per iteration). Default is 20.

257 finite_diff_rel_step : None or array_like, optional

258 If `jac in ['2-point', '3-point', 'cs']` the relative step size to

259 use for numerical approximation of the jacobian. The absolute step

260 size is computed as ``h = rel_step * sign(x) * max(1, abs(x))``,

261 possibly adjusted to fit into the bounds. For ``method='3-point'``

262 the sign of `h` is ignored. If None (default) then step is selected

263 automatically.

264

265 Notes

266 -----

267 The option `ftol` is exposed via the `scipy.optimize.minimize` interface,

268 but calling `scipy.optimize.fmin_l_bfgs_b` directly exposes `factr`. The

269 relationship between the two is ``ftol = factr * numpy.finfo(float).eps``.

270 I.e., `factr` multiplies the default machine floating-point precision to

271 arrive at `ftol`.

272

273 """

274 _check_unknown_options(unknown_options)

275 m = maxcor

276 pgtol = gtol

277 factr = ftol / np.finfo(float).eps

278

279 x0 = asarray(x0).ravel()

280 n, = x0.shape

281

282 if bounds is None:

283 bounds = [(None, None)] * n

284 if len(bounds) != n:

285 raise ValueError('length of x0 != length of bounds')

286

287 # unbounded variables must use None, not +-inf, for optimizer to work properly

288 bounds = [(None if l == -np.inf else l, None if u == np.inf else u) for l, u in bounds]

289 # LBFGSB is sent 'old-style' bounds, 'new-style' bounds are required by

290 # approx_derivative and ScalarFunction

291 new_bounds = old_bound_to_new(bounds)

292

293 # check bounds

294 if (new_bounds[0] > new_bounds[1]).any():

295 raise ValueError("LBFGSB - one of the lower bounds is greater than an upper bound.")

296

297 # initial vector must lie within the bounds. Otherwise ScalarFunction and

298 # approx_derivative will cause problems

299 x0 = np.clip(x0, new_bounds[0], new_bounds[1])

300

301 if disp is not None:

302 if disp == 0:

303 iprint = -1

304 else:

305 iprint = disp

306

307 sf = _prepare_scalar_function(fun, x0, jac=jac, args=args, epsilon=eps,

308 bounds=new_bounds,

309 finite_diff_rel_step=finite_diff_rel_step)

310

311 func_and_grad = sf.fun_and_grad

312

313 fortran_int = _lbfgsb.types.intvar.dtype

314

315 nbd = zeros(n, fortran_int)

316 low_bnd = zeros(n, float64)

317 upper_bnd = zeros(n, float64)

318 bounds_map = {(None, None): 0,

319 (1, None): 1,

320 (1, 1): 2,

321 (None, 1): 3}

322 for i in range(0, n):

323 l, u = bounds[i]

324 if l is not None:

325 low_bnd[i] = l

326 l = 1

327 if u is not None:

328 upper_bnd[i] = u

329 u = 1

330 nbd[i] = bounds_map[l, u]

331

332 if not maxls > 0:

333 raise ValueError('maxls must be positive.')

334

335 x = array(x0, float64)

336 f = array(0.0, float64)

337 g = zeros((n,), float64)

338 wa = zeros(2*m*n + 5*n + 11*m*m + 8*m, float64)

339 iwa = zeros(3*n, fortran_int)

340 task = zeros(1, 'S60')

341 csave = zeros(1, 'S60')

342 lsave = zeros(4, fortran_int)

343 isave = zeros(44, fortran_int)

344 dsave = zeros(29, float64)

345

346 task[:] = 'START'

347

348 n_iterations = 0

349

350 while 1:

351 # g may become float32 if a user provides a function that calculates

352 # the Jacobian in float32 (see gh-18730). The underlying Fortran code

353 # expects float64, so upcast it

354 g = g.astype(np.float64)

355 # x, f, g, wa, iwa, task, csave, lsave, isave, dsave = \

356 _lbfgsb.setulb(m, x, low_bnd, upper_bnd, nbd, f, g, factr,

357 pgtol, wa, iwa, task, iprint, csave, lsave,

358 isave, dsave, maxls)

359 task_str = task.tobytes()

360 if task_str.startswith(b'FG'):

361 # The minimization routine wants f and g at the current x.

362 # Note that interruptions due to maxfun are postponed

363 # until the completion of the current minimization iteration.

364 # Overwrite f and g:

365 f, g = func_and_grad(x)

366 elif task_str.startswith(b'NEW_X'):

367 # new iteration

368 n_iterations += 1

369

370 intermediate_result = OptimizeResult(x=x, fun=f)

371 if _call_callback_maybe_halt(callback, intermediate_result):

372 task[:] = 'STOP: CALLBACK REQUESTED HALT'

373 if n_iterations >= maxiter:

374 task[:] = 'STOP: TOTAL NO. of ITERATIONS REACHED LIMIT'

375 elif sf.nfev > maxfun:

376 task[:] = ('STOP: TOTAL NO. of f AND g EVALUATIONS '

377 'EXCEEDS LIMIT')

378 else:

379 break

380

381 task_str = task.tobytes().strip(b'\x00').strip()

382 if task_str.startswith(b'CONV'):

383 warnflag = 0

384 elif sf.nfev > maxfun or n_iterations >= maxiter:

385 warnflag = 1

386 else:

387 warnflag = 2

388

389 # These two portions of the workspace are described in the mainlb

390 # subroutine in lbfgsb.f. See line 363.

391 s = wa[0: m*n].reshape(m, n)

392 y = wa[m*n: 2*m*n].reshape(m, n)

393

394 # See lbfgsb.f line 160 for this portion of the workspace.

395 # isave(31) = the total number of BFGS updates prior the current iteration;

396 n_bfgs_updates = isave[30]

397

398 n_corrs = min(n_bfgs_updates, maxcor)

399 hess_inv = LbfgsInvHessProduct(s[:n_corrs], y[:n_corrs])

400

401 task_str = task_str.decode()

402 return OptimizeResult(fun=f, jac=g, nfev=sf.nfev,

403 njev=sf.ngev,

404 nit=n_iterations, status=warnflag, message=task_str,

405 x=x, success=(warnflag == 0), hess_inv=hess_inv)

406

407

408class LbfgsInvHessProduct(LinearOperator):

409 """Linear operator for the L-BFGS approximate inverse Hessian.

410

411 This operator computes the product of a vector with the approximate inverse

412 of the Hessian of the objective function, using the L-BFGS limited

413 memory approximation to the inverse Hessian, accumulated during the

414 optimization.

415

416 Objects of this class implement the ``scipy.sparse.linalg.LinearOperator``

417 interface.

418

419 Parameters

420 ----------

421 sk : array_like, shape=(n_corr, n)

422 Array of `n_corr` most recent updates to the solution vector.

423 (See [1]).

424 yk : array_like, shape=(n_corr, n)

425 Array of `n_corr` most recent updates to the gradient. (See [1]).

426

427 References

428 ----------

429 .. [1] Nocedal, Jorge. "Updating quasi-Newton matrices with limited

430 storage." Mathematics of computation 35.151 (1980): 773-782.

431

432 """

433

434 def __init__(self, sk, yk):

435 """Construct the operator."""

436 if sk.shape != yk.shape or sk.ndim != 2:

437 raise ValueError('sk and yk must have matching shape, (n_corrs, n)')

438 n_corrs, n = sk.shape

439

440 super().__init__(dtype=np.float64, shape=(n, n))

441

442 self.sk = sk

443 self.yk = yk

444 self.n_corrs = n_corrs

445 self.rho = 1 / np.einsum('ij,ij->i', sk, yk)

446

447 def _matvec(self, x):

448 """Efficient matrix-vector multiply with the BFGS matrices.

449

450 This calculation is described in Section (4) of [1].

451

452 Parameters

453 ----------

454 x : ndarray

455 An array with shape (n,) or (n,1).

456

457 Returns

458 -------

459 y : ndarray

460 The matrix-vector product

461

462 """

463 s, y, n_corrs, rho = self.sk, self.yk, self.n_corrs, self.rho

464 q = np.array(x, dtype=self.dtype, copy=True)

465 if q.ndim == 2 and q.shape[1] == 1:

466 q = q.reshape(-1)

467

468 alpha = np.empty(n_corrs)

469

470 for i in range(n_corrs-1, -1, -1):

471 alpha[i] = rho[i] * np.dot(s[i], q)

472 q = q - alpha[i]*y[i]

473

474 r = q

475 for i in range(n_corrs):

476 beta = rho[i] * np.dot(y[i], r)

477 r = r + s[i] * (alpha[i] - beta)

478

479 return r

480

481 def todense(self):

482 """Return a dense array representation of this operator.

483

484 Returns

485 -------

486 arr : ndarray, shape=(n, n)

487 An array with the same shape and containing

488 the same data represented by this `LinearOperator`.

489

490 """

491 s, y, n_corrs, rho = self.sk, self.yk, self.n_corrs, self.rho

492 I = np.eye(*self.shape, dtype=self.dtype)

493 Hk = I

494

495 for i in range(n_corrs):

496 A1 = I - s[i][:, np.newaxis] * y[i][np.newaxis, :] * rho[i]

497 A2 = I - y[i][:, np.newaxis] * s[i][np.newaxis, :] * rho[i]

498

499 Hk = np.dot(A1, np.dot(Hk, A2)) + (rho[i] * s[i][:, np.newaxis] *

500 s[i][np.newaxis, :])

501 return Hk

Coverage for /usr/lib/python3/dist-packages/scipy/optimize/_lbfgsb_py.py: 10%

135 statements