Coverage for /usr/lib/python3/dist-packages/scipy/optimize/_lsq/least_squares.py: 10%
259 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""Generic interface for least-squares minimization."""
2from warnings import warn
4import numpy as np
5from numpy.linalg import norm
7from scipy.sparse import issparse
8from scipy.sparse.linalg import LinearOperator
9from scipy.optimize import _minpack, OptimizeResult
10from scipy.optimize._numdiff import approx_derivative, group_columns
11from scipy.optimize._minimize import Bounds
13from .trf import trf
14from .dogbox import dogbox
15from .common import EPS, in_bounds, make_strictly_feasible
18TERMINATION_MESSAGES = {
19 -1: "Improper input parameters status returned from `leastsq`",
20 0: "The maximum number of function evaluations is exceeded.",
21 1: "`gtol` termination condition is satisfied.",
22 2: "`ftol` termination condition is satisfied.",
23 3: "`xtol` termination condition is satisfied.",
24 4: "Both `ftol` and `xtol` termination conditions are satisfied."
25}
28FROM_MINPACK_TO_COMMON = {
29 0: -1, # Improper input parameters from MINPACK.
30 1: 2,
31 2: 3,
32 3: 4,
33 4: 1,
34 5: 0
35 # There are 6, 7, 8 for too small tolerance parameters,
36 # but we guard against it by checking ftol, xtol, gtol beforehand.
37}
40def call_minpack(fun, x0, jac, ftol, xtol, gtol, max_nfev, x_scale, diff_step):
41 n = x0.size
43 if diff_step is None:
44 epsfcn = EPS
45 else:
46 epsfcn = diff_step**2
48 # Compute MINPACK's `diag`, which is inverse of our `x_scale` and
49 # ``x_scale='jac'`` corresponds to ``diag=None``.
50 if isinstance(x_scale, str) and x_scale == 'jac':
51 diag = None
52 else:
53 diag = 1 / x_scale
55 full_output = True
56 col_deriv = False
57 factor = 100.0
59 if jac is None:
60 if max_nfev is None:
61 # n squared to account for Jacobian evaluations.
62 max_nfev = 100 * n * (n + 1)
63 x, info, status = _minpack._lmdif(
64 fun, x0, (), full_output, ftol, xtol, gtol,
65 max_nfev, epsfcn, factor, diag)
66 else:
67 if max_nfev is None:
68 max_nfev = 100 * n
69 x, info, status = _minpack._lmder(
70 fun, jac, x0, (), full_output, col_deriv,
71 ftol, xtol, gtol, max_nfev, factor, diag)
73 f = info['fvec']
75 if callable(jac):
76 J = jac(x)
77 else:
78 J = np.atleast_2d(approx_derivative(fun, x))
80 cost = 0.5 * np.dot(f, f)
81 g = J.T.dot(f)
82 g_norm = norm(g, ord=np.inf)
84 nfev = info['nfev']
85 njev = info.get('njev', None)
87 status = FROM_MINPACK_TO_COMMON[status]
88 active_mask = np.zeros_like(x0, dtype=int)
90 return OptimizeResult(
91 x=x, cost=cost, fun=f, jac=J, grad=g, optimality=g_norm,
92 active_mask=active_mask, nfev=nfev, njev=njev, status=status)
95def prepare_bounds(bounds, n):
96 lb, ub = (np.asarray(b, dtype=float) for b in bounds)
97 if lb.ndim == 0:
98 lb = np.resize(lb, n)
100 if ub.ndim == 0:
101 ub = np.resize(ub, n)
103 return lb, ub
106def check_tolerance(ftol, xtol, gtol, method):
107 def check(tol, name):
108 if tol is None:
109 tol = 0
110 elif tol < EPS:
111 warn("Setting `{}` below the machine epsilon ({:.2e}) effectively "
112 "disables the corresponding termination condition."
113 .format(name, EPS))
114 return tol
116 ftol = check(ftol, "ftol")
117 xtol = check(xtol, "xtol")
118 gtol = check(gtol, "gtol")
120 if method == "lm" and (ftol < EPS or xtol < EPS or gtol < EPS):
121 raise ValueError("All tolerances must be higher than machine epsilon "
122 "({:.2e}) for method 'lm'.".format(EPS))
123 elif ftol < EPS and xtol < EPS and gtol < EPS:
124 raise ValueError("At least one of the tolerances must be higher than "
125 "machine epsilon ({:.2e}).".format(EPS))
127 return ftol, xtol, gtol
130def check_x_scale(x_scale, x0):
131 if isinstance(x_scale, str) and x_scale == 'jac':
132 return x_scale
134 try:
135 x_scale = np.asarray(x_scale, dtype=float)
136 valid = np.all(np.isfinite(x_scale)) and np.all(x_scale > 0)
137 except (ValueError, TypeError):
138 valid = False
140 if not valid:
141 raise ValueError("`x_scale` must be 'jac' or array_like with "
142 "positive numbers.")
144 if x_scale.ndim == 0:
145 x_scale = np.resize(x_scale, x0.shape)
147 if x_scale.shape != x0.shape:
148 raise ValueError("Inconsistent shapes between `x_scale` and `x0`.")
150 return x_scale
153def check_jac_sparsity(jac_sparsity, m, n):
154 if jac_sparsity is None:
155 return None
157 if not issparse(jac_sparsity):
158 jac_sparsity = np.atleast_2d(jac_sparsity)
160 if jac_sparsity.shape != (m, n):
161 raise ValueError("`jac_sparsity` has wrong shape.")
163 return jac_sparsity, group_columns(jac_sparsity)
166# Loss functions.
169def huber(z, rho, cost_only):
170 mask = z <= 1
171 rho[0, mask] = z[mask]
172 rho[0, ~mask] = 2 * z[~mask]**0.5 - 1
173 if cost_only:
174 return
175 rho[1, mask] = 1
176 rho[1, ~mask] = z[~mask]**-0.5
177 rho[2, mask] = 0
178 rho[2, ~mask] = -0.5 * z[~mask]**-1.5
181def soft_l1(z, rho, cost_only):
182 t = 1 + z
183 rho[0] = 2 * (t**0.5 - 1)
184 if cost_only:
185 return
186 rho[1] = t**-0.5
187 rho[2] = -0.5 * t**-1.5
190def cauchy(z, rho, cost_only):
191 rho[0] = np.log1p(z)
192 if cost_only:
193 return
194 t = 1 + z
195 rho[1] = 1 / t
196 rho[2] = -1 / t**2
199def arctan(z, rho, cost_only):
200 rho[0] = np.arctan(z)
201 if cost_only:
202 return
203 t = 1 + z**2
204 rho[1] = 1 / t
205 rho[2] = -2 * z / t**2
208IMPLEMENTED_LOSSES = dict(linear=None, huber=huber, soft_l1=soft_l1,
209 cauchy=cauchy, arctan=arctan)
212def construct_loss_function(m, loss, f_scale):
213 if loss == 'linear':
214 return None
216 if not callable(loss):
217 loss = IMPLEMENTED_LOSSES[loss]
218 rho = np.empty((3, m))
220 def loss_function(f, cost_only=False):
221 z = (f / f_scale) ** 2
222 loss(z, rho, cost_only=cost_only)
223 if cost_only:
224 return 0.5 * f_scale ** 2 * np.sum(rho[0])
225 rho[0] *= f_scale ** 2
226 rho[2] /= f_scale ** 2
227 return rho
228 else:
229 def loss_function(f, cost_only=False):
230 z = (f / f_scale) ** 2
231 rho = loss(z)
232 if cost_only:
233 return 0.5 * f_scale ** 2 * np.sum(rho[0])
234 rho[0] *= f_scale ** 2
235 rho[2] /= f_scale ** 2
236 return rho
238 return loss_function
241def least_squares(
242 fun, x0, jac='2-point', bounds=(-np.inf, np.inf), method='trf',
243 ftol=1e-8, xtol=1e-8, gtol=1e-8, x_scale=1.0, loss='linear',
244 f_scale=1.0, diff_step=None, tr_solver=None, tr_options={},
245 jac_sparsity=None, max_nfev=None, verbose=0, args=(), kwargs={}):
246 """Solve a nonlinear least-squares problem with bounds on the variables.
248 Given the residuals f(x) (an m-D real function of n real
249 variables) and the loss function rho(s) (a scalar function), `least_squares`
250 finds a local minimum of the cost function F(x)::
252 minimize F(x) = 0.5 * sum(rho(f_i(x)**2), i = 0, ..., m - 1)
253 subject to lb <= x <= ub
255 The purpose of the loss function rho(s) is to reduce the influence of
256 outliers on the solution.
258 Parameters
259 ----------
260 fun : callable
261 Function which computes the vector of residuals, with the signature
262 ``fun(x, *args, **kwargs)``, i.e., the minimization proceeds with
263 respect to its first argument. The argument ``x`` passed to this
264 function is an ndarray of shape (n,) (never a scalar, even for n=1).
265 It must allocate and return a 1-D array_like of shape (m,) or a scalar.
266 If the argument ``x`` is complex or the function ``fun`` returns
267 complex residuals, it must be wrapped in a real function of real
268 arguments, as shown at the end of the Examples section.
269 x0 : array_like with shape (n,) or float
270 Initial guess on independent variables. If float, it will be treated
271 as a 1-D array with one element. When `method` is 'trf', the initial
272 guess might be slightly adjusted to lie sufficiently within the given
273 `bounds`.
274 jac : {'2-point', '3-point', 'cs', callable}, optional
275 Method of computing the Jacobian matrix (an m-by-n matrix, where
276 element (i, j) is the partial derivative of f[i] with respect to
277 x[j]). The keywords select a finite difference scheme for numerical
278 estimation. The scheme '3-point' is more accurate, but requires
279 twice as many operations as '2-point' (default). The scheme 'cs'
280 uses complex steps, and while potentially the most accurate, it is
281 applicable only when `fun` correctly handles complex inputs and
282 can be analytically continued to the complex plane. Method 'lm'
283 always uses the '2-point' scheme. If callable, it is used as
284 ``jac(x, *args, **kwargs)`` and should return a good approximation
285 (or the exact value) for the Jacobian as an array_like (np.atleast_2d
286 is applied), a sparse matrix (csr_matrix preferred for performance) or
287 a `scipy.sparse.linalg.LinearOperator`.
288 bounds : 2-tuple of array_like or `Bounds`, optional
289 There are two ways to specify bounds:
291 1. Instance of `Bounds` class
292 2. Lower and upper bounds on independent variables. Defaults to no
293 bounds. Each array must match the size of `x0` or be a scalar,
294 in the latter case a bound will be the same for all variables.
295 Use ``np.inf`` with an appropriate sign to disable bounds on all
296 or some variables.
297 method : {'trf', 'dogbox', 'lm'}, optional
298 Algorithm to perform minimization.
300 * 'trf' : Trust Region Reflective algorithm, particularly suitable
301 for large sparse problems with bounds. Generally robust method.
302 * 'dogbox' : dogleg algorithm with rectangular trust regions,
303 typical use case is small problems with bounds. Not recommended
304 for problems with rank-deficient Jacobian.
305 * 'lm' : Levenberg-Marquardt algorithm as implemented in MINPACK.
306 Doesn't handle bounds and sparse Jacobians. Usually the most
307 efficient method for small unconstrained problems.
309 Default is 'trf'. See Notes for more information.
310 ftol : float or None, optional
311 Tolerance for termination by the change of the cost function. Default
312 is 1e-8. The optimization process is stopped when ``dF < ftol * F``,
313 and there was an adequate agreement between a local quadratic model and
314 the true model in the last step.
316 If None and 'method' is not 'lm', the termination by this condition is
317 disabled. If 'method' is 'lm', this tolerance must be higher than
318 machine epsilon.
319 xtol : float or None, optional
320 Tolerance for termination by the change of the independent variables.
321 Default is 1e-8. The exact condition depends on the `method` used:
323 * For 'trf' and 'dogbox' : ``norm(dx) < xtol * (xtol + norm(x))``.
324 * For 'lm' : ``Delta < xtol * norm(xs)``, where ``Delta`` is
325 a trust-region radius and ``xs`` is the value of ``x``
326 scaled according to `x_scale` parameter (see below).
328 If None and 'method' is not 'lm', the termination by this condition is
329 disabled. If 'method' is 'lm', this tolerance must be higher than
330 machine epsilon.
331 gtol : float or None, optional
332 Tolerance for termination by the norm of the gradient. Default is 1e-8.
333 The exact condition depends on a `method` used:
335 * For 'trf' : ``norm(g_scaled, ord=np.inf) < gtol``, where
336 ``g_scaled`` is the value of the gradient scaled to account for
337 the presence of the bounds [STIR]_.
338 * For 'dogbox' : ``norm(g_free, ord=np.inf) < gtol``, where
339 ``g_free`` is the gradient with respect to the variables which
340 are not in the optimal state on the boundary.
341 * For 'lm' : the maximum absolute value of the cosine of angles
342 between columns of the Jacobian and the residual vector is less
343 than `gtol`, or the residual vector is zero.
345 If None and 'method' is not 'lm', the termination by this condition is
346 disabled. If 'method' is 'lm', this tolerance must be higher than
347 machine epsilon.
348 x_scale : array_like or 'jac', optional
349 Characteristic scale of each variable. Setting `x_scale` is equivalent
350 to reformulating the problem in scaled variables ``xs = x / x_scale``.
351 An alternative view is that the size of a trust region along jth
352 dimension is proportional to ``x_scale[j]``. Improved convergence may
353 be achieved by setting `x_scale` such that a step of a given size
354 along any of the scaled variables has a similar effect on the cost
355 function. If set to 'jac', the scale is iteratively updated using the
356 inverse norms of the columns of the Jacobian matrix (as described in
357 [JJMore]_).
358 loss : str or callable, optional
359 Determines the loss function. The following keyword values are allowed:
361 * 'linear' (default) : ``rho(z) = z``. Gives a standard
362 least-squares problem.
363 * 'soft_l1' : ``rho(z) = 2 * ((1 + z)**0.5 - 1)``. The smooth
364 approximation of l1 (absolute value) loss. Usually a good
365 choice for robust least squares.
366 * 'huber' : ``rho(z) = z if z <= 1 else 2*z**0.5 - 1``. Works
367 similarly to 'soft_l1'.
368 * 'cauchy' : ``rho(z) = ln(1 + z)``. Severely weakens outliers
369 influence, but may cause difficulties in optimization process.
370 * 'arctan' : ``rho(z) = arctan(z)``. Limits a maximum loss on
371 a single residual, has properties similar to 'cauchy'.
373 If callable, it must take a 1-D ndarray ``z=f**2`` and return an
374 array_like with shape (3, m) where row 0 contains function values,
375 row 1 contains first derivatives and row 2 contains second
376 derivatives. Method 'lm' supports only 'linear' loss.
377 f_scale : float, optional
378 Value of soft margin between inlier and outlier residuals, default
379 is 1.0. The loss function is evaluated as follows
380 ``rho_(f**2) = C**2 * rho(f**2 / C**2)``, where ``C`` is `f_scale`,
381 and ``rho`` is determined by `loss` parameter. This parameter has
382 no effect with ``loss='linear'``, but for other `loss` values it is
383 of crucial importance.
384 max_nfev : None or int, optional
385 Maximum number of function evaluations before the termination.
386 If None (default), the value is chosen automatically:
388 * For 'trf' and 'dogbox' : 100 * n.
389 * For 'lm' : 100 * n if `jac` is callable and 100 * n * (n + 1)
390 otherwise (because 'lm' counts function calls in Jacobian
391 estimation).
393 diff_step : None or array_like, optional
394 Determines the relative step size for the finite difference
395 approximation of the Jacobian. The actual step is computed as
396 ``x * diff_step``. If None (default), then `diff_step` is taken to be
397 a conventional "optimal" power of machine epsilon for the finite
398 difference scheme used [NR]_.
399 tr_solver : {None, 'exact', 'lsmr'}, optional
400 Method for solving trust-region subproblems, relevant only for 'trf'
401 and 'dogbox' methods.
403 * 'exact' is suitable for not very large problems with dense
404 Jacobian matrices. The computational complexity per iteration is
405 comparable to a singular value decomposition of the Jacobian
406 matrix.
407 * 'lsmr' is suitable for problems with sparse and large Jacobian
408 matrices. It uses the iterative procedure
409 `scipy.sparse.linalg.lsmr` for finding a solution of a linear
410 least-squares problem and only requires matrix-vector product
411 evaluations.
413 If None (default), the solver is chosen based on the type of Jacobian
414 returned on the first iteration.
415 tr_options : dict, optional
416 Keyword options passed to trust-region solver.
418 * ``tr_solver='exact'``: `tr_options` are ignored.
419 * ``tr_solver='lsmr'``: options for `scipy.sparse.linalg.lsmr`.
420 Additionally, ``method='trf'`` supports 'regularize' option
421 (bool, default is True), which adds a regularization term to the
422 normal equation, which improves convergence if the Jacobian is
423 rank-deficient [Byrd]_ (eq. 3.4).
425 jac_sparsity : {None, array_like, sparse matrix}, optional
426 Defines the sparsity structure of the Jacobian matrix for finite
427 difference estimation, its shape must be (m, n). If the Jacobian has
428 only few non-zero elements in *each* row, providing the sparsity
429 structure will greatly speed up the computations [Curtis]_. A zero
430 entry means that a corresponding element in the Jacobian is identically
431 zero. If provided, forces the use of 'lsmr' trust-region solver.
432 If None (default), then dense differencing will be used. Has no effect
433 for 'lm' method.
434 verbose : {0, 1, 2}, optional
435 Level of algorithm's verbosity:
437 * 0 (default) : work silently.
438 * 1 : display a termination report.
439 * 2 : display progress during iterations (not supported by 'lm'
440 method).
442 args, kwargs : tuple and dict, optional
443 Additional arguments passed to `fun` and `jac`. Both empty by default.
444 The calling signature is ``fun(x, *args, **kwargs)`` and the same for
445 `jac`.
447 Returns
448 -------
449 result : OptimizeResult
450 `OptimizeResult` with the following fields defined:
452 x : ndarray, shape (n,)
453 Solution found.
454 cost : float
455 Value of the cost function at the solution.
456 fun : ndarray, shape (m,)
457 Vector of residuals at the solution.
458 jac : ndarray, sparse matrix or LinearOperator, shape (m, n)
459 Modified Jacobian matrix at the solution, in the sense that J^T J
460 is a Gauss-Newton approximation of the Hessian of the cost function.
461 The type is the same as the one used by the algorithm.
462 grad : ndarray, shape (m,)
463 Gradient of the cost function at the solution.
464 optimality : float
465 First-order optimality measure. In unconstrained problems, it is
466 always the uniform norm of the gradient. In constrained problems,
467 it is the quantity which was compared with `gtol` during iterations.
468 active_mask : ndarray of int, shape (n,)
469 Each component shows whether a corresponding constraint is active
470 (that is, whether a variable is at the bound):
472 * 0 : a constraint is not active.
473 * -1 : a lower bound is active.
474 * 1 : an upper bound is active.
476 Might be somewhat arbitrary for 'trf' method as it generates a
477 sequence of strictly feasible iterates and `active_mask` is
478 determined within a tolerance threshold.
479 nfev : int
480 Number of function evaluations done. Methods 'trf' and 'dogbox' do
481 not count function calls for numerical Jacobian approximation, as
482 opposed to 'lm' method.
483 njev : int or None
484 Number of Jacobian evaluations done. If numerical Jacobian
485 approximation is used in 'lm' method, it is set to None.
486 status : int
487 The reason for algorithm termination:
489 * -1 : improper input parameters status returned from MINPACK.
490 * 0 : the maximum number of function evaluations is exceeded.
491 * 1 : `gtol` termination condition is satisfied.
492 * 2 : `ftol` termination condition is satisfied.
493 * 3 : `xtol` termination condition is satisfied.
494 * 4 : Both `ftol` and `xtol` termination conditions are satisfied.
496 message : str
497 Verbal description of the termination reason.
498 success : bool
499 True if one of the convergence criteria is satisfied (`status` > 0).
501 See Also
502 --------
503 leastsq : A legacy wrapper for the MINPACK implementation of the
504 Levenberg-Marquadt algorithm.
505 curve_fit : Least-squares minimization applied to a curve-fitting problem.
507 Notes
508 -----
509 Method 'lm' (Levenberg-Marquardt) calls a wrapper over least-squares
510 algorithms implemented in MINPACK (lmder, lmdif). It runs the
511 Levenberg-Marquardt algorithm formulated as a trust-region type algorithm.
512 The implementation is based on paper [JJMore]_, it is very robust and
513 efficient with a lot of smart tricks. It should be your first choice
514 for unconstrained problems. Note that it doesn't support bounds. Also,
515 it doesn't work when m < n.
517 Method 'trf' (Trust Region Reflective) is motivated by the process of
518 solving a system of equations, which constitute the first-order optimality
519 condition for a bound-constrained minimization problem as formulated in
520 [STIR]_. The algorithm iteratively solves trust-region subproblems
521 augmented by a special diagonal quadratic term and with trust-region shape
522 determined by the distance from the bounds and the direction of the
523 gradient. This enhancements help to avoid making steps directly into bounds
524 and efficiently explore the whole space of variables. To further improve
525 convergence, the algorithm considers search directions reflected from the
526 bounds. To obey theoretical requirements, the algorithm keeps iterates
527 strictly feasible. With dense Jacobians trust-region subproblems are
528 solved by an exact method very similar to the one described in [JJMore]_
529 (and implemented in MINPACK). The difference from the MINPACK
530 implementation is that a singular value decomposition of a Jacobian
531 matrix is done once per iteration, instead of a QR decomposition and series
532 of Givens rotation eliminations. For large sparse Jacobians a 2-D subspace
533 approach of solving trust-region subproblems is used [STIR]_, [Byrd]_.
534 The subspace is spanned by a scaled gradient and an approximate
535 Gauss-Newton solution delivered by `scipy.sparse.linalg.lsmr`. When no
536 constraints are imposed the algorithm is very similar to MINPACK and has
537 generally comparable performance. The algorithm works quite robust in
538 unbounded and bounded problems, thus it is chosen as a default algorithm.
540 Method 'dogbox' operates in a trust-region framework, but considers
541 rectangular trust regions as opposed to conventional ellipsoids [Voglis]_.
542 The intersection of a current trust region and initial bounds is again
543 rectangular, so on each iteration a quadratic minimization problem subject
544 to bound constraints is solved approximately by Powell's dogleg method
545 [NumOpt]_. The required Gauss-Newton step can be computed exactly for
546 dense Jacobians or approximately by `scipy.sparse.linalg.lsmr` for large
547 sparse Jacobians. The algorithm is likely to exhibit slow convergence when
548 the rank of Jacobian is less than the number of variables. The algorithm
549 often outperforms 'trf' in bounded problems with a small number of
550 variables.
552 Robust loss functions are implemented as described in [BA]_. The idea
553 is to modify a residual vector and a Jacobian matrix on each iteration
554 such that computed gradient and Gauss-Newton Hessian approximation match
555 the true gradient and Hessian approximation of the cost function. Then
556 the algorithm proceeds in a normal way, i.e., robust loss functions are
557 implemented as a simple wrapper over standard least-squares algorithms.
559 .. versionadded:: 0.17.0
561 References
562 ----------
563 .. [STIR] M. A. Branch, T. F. Coleman, and Y. Li, "A Subspace, Interior,
564 and Conjugate Gradient Method for Large-Scale Bound-Constrained
565 Minimization Problems," SIAM Journal on Scientific Computing,
566 Vol. 21, Number 1, pp 1-23, 1999.
567 .. [NR] William H. Press et. al., "Numerical Recipes. The Art of Scientific
568 Computing. 3rd edition", Sec. 5.7.
569 .. [Byrd] R. H. Byrd, R. B. Schnabel and G. A. Shultz, "Approximate
570 solution of the trust region problem by minimization over
571 two-dimensional subspaces", Math. Programming, 40, pp. 247-263,
572 1988.
573 .. [Curtis] A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of
574 sparse Jacobian matrices", Journal of the Institute of
575 Mathematics and its Applications, 13, pp. 117-120, 1974.
576 .. [JJMore] J. J. More, "The Levenberg-Marquardt Algorithm: Implementation
577 and Theory," Numerical Analysis, ed. G. A. Watson, Lecture
578 Notes in Mathematics 630, Springer Verlag, pp. 105-116, 1977.
579 .. [Voglis] C. Voglis and I. E. Lagaris, "A Rectangular Trust Region
580 Dogleg Approach for Unconstrained and Bound Constrained
581 Nonlinear Optimization", WSEAS International Conference on
582 Applied Mathematics, Corfu, Greece, 2004.
583 .. [NumOpt] J. Nocedal and S. J. Wright, "Numerical optimization,
584 2nd edition", Chapter 4.
585 .. [BA] B. Triggs et. al., "Bundle Adjustment - A Modern Synthesis",
586 Proceedings of the International Workshop on Vision Algorithms:
587 Theory and Practice, pp. 298-372, 1999.
589 Examples
590 --------
591 In this example we find a minimum of the Rosenbrock function without bounds
592 on independent variables.
594 >>> import numpy as np
595 >>> def fun_rosenbrock(x):
596 ... return np.array([10 * (x[1] - x[0]**2), (1 - x[0])])
598 Notice that we only provide the vector of the residuals. The algorithm
599 constructs the cost function as a sum of squares of the residuals, which
600 gives the Rosenbrock function. The exact minimum is at ``x = [1.0, 1.0]``.
602 >>> from scipy.optimize import least_squares
603 >>> x0_rosenbrock = np.array([2, 2])
604 >>> res_1 = least_squares(fun_rosenbrock, x0_rosenbrock)
605 >>> res_1.x
606 array([ 1., 1.])
607 >>> res_1.cost
608 9.8669242910846867e-30
609 >>> res_1.optimality
610 8.8928864934219529e-14
612 We now constrain the variables, in such a way that the previous solution
613 becomes infeasible. Specifically, we require that ``x[1] >= 1.5``, and
614 ``x[0]`` left unconstrained. To this end, we specify the `bounds` parameter
615 to `least_squares` in the form ``bounds=([-np.inf, 1.5], np.inf)``.
617 We also provide the analytic Jacobian:
619 >>> def jac_rosenbrock(x):
620 ... return np.array([
621 ... [-20 * x[0], 10],
622 ... [-1, 0]])
624 Putting this all together, we see that the new solution lies on the bound:
626 >>> res_2 = least_squares(fun_rosenbrock, x0_rosenbrock, jac_rosenbrock,
627 ... bounds=([-np.inf, 1.5], np.inf))
628 >>> res_2.x
629 array([ 1.22437075, 1.5 ])
630 >>> res_2.cost
631 0.025213093946805685
632 >>> res_2.optimality
633 1.5885401433157753e-07
635 Now we solve a system of equations (i.e., the cost function should be zero
636 at a minimum) for a Broyden tridiagonal vector-valued function of 100000
637 variables:
639 >>> def fun_broyden(x):
640 ... f = (3 - x) * x + 1
641 ... f[1:] -= x[:-1]
642 ... f[:-1] -= 2 * x[1:]
643 ... return f
645 The corresponding Jacobian matrix is sparse. We tell the algorithm to
646 estimate it by finite differences and provide the sparsity structure of
647 Jacobian to significantly speed up this process.
649 >>> from scipy.sparse import lil_matrix
650 >>> def sparsity_broyden(n):
651 ... sparsity = lil_matrix((n, n), dtype=int)
652 ... i = np.arange(n)
653 ... sparsity[i, i] = 1
654 ... i = np.arange(1, n)
655 ... sparsity[i, i - 1] = 1
656 ... i = np.arange(n - 1)
657 ... sparsity[i, i + 1] = 1
658 ... return sparsity
659 ...
660 >>> n = 100000
661 >>> x0_broyden = -np.ones(n)
662 ...
663 >>> res_3 = least_squares(fun_broyden, x0_broyden,
664 ... jac_sparsity=sparsity_broyden(n))
665 >>> res_3.cost
666 4.5687069299604613e-23
667 >>> res_3.optimality
668 1.1650454296851518e-11
670 Let's also solve a curve fitting problem using robust loss function to
671 take care of outliers in the data. Define the model function as
672 ``y = a + b * exp(c * t)``, where t is a predictor variable, y is an
673 observation and a, b, c are parameters to estimate.
675 First, define the function which generates the data with noise and
676 outliers, define the model parameters, and generate data:
678 >>> from numpy.random import default_rng
679 >>> rng = default_rng()
680 >>> def gen_data(t, a, b, c, noise=0., n_outliers=0, seed=None):
681 ... rng = default_rng(seed)
682 ...
683 ... y = a + b * np.exp(t * c)
684 ...
685 ... error = noise * rng.standard_normal(t.size)
686 ... outliers = rng.integers(0, t.size, n_outliers)
687 ... error[outliers] *= 10
688 ...
689 ... return y + error
690 ...
691 >>> a = 0.5
692 >>> b = 2.0
693 >>> c = -1
694 >>> t_min = 0
695 >>> t_max = 10
696 >>> n_points = 15
697 ...
698 >>> t_train = np.linspace(t_min, t_max, n_points)
699 >>> y_train = gen_data(t_train, a, b, c, noise=0.1, n_outliers=3)
701 Define function for computing residuals and initial estimate of
702 parameters.
704 >>> def fun(x, t, y):
705 ... return x[0] + x[1] * np.exp(x[2] * t) - y
706 ...
707 >>> x0 = np.array([1.0, 1.0, 0.0])
709 Compute a standard least-squares solution:
711 >>> res_lsq = least_squares(fun, x0, args=(t_train, y_train))
713 Now compute two solutions with two different robust loss functions. The
714 parameter `f_scale` is set to 0.1, meaning that inlier residuals should
715 not significantly exceed 0.1 (the noise level used).
717 >>> res_soft_l1 = least_squares(fun, x0, loss='soft_l1', f_scale=0.1,
718 ... args=(t_train, y_train))
719 >>> res_log = least_squares(fun, x0, loss='cauchy', f_scale=0.1,
720 ... args=(t_train, y_train))
722 And, finally, plot all the curves. We see that by selecting an appropriate
723 `loss` we can get estimates close to optimal even in the presence of
724 strong outliers. But keep in mind that generally it is recommended to try
725 'soft_l1' or 'huber' losses first (if at all necessary) as the other two
726 options may cause difficulties in optimization process.
728 >>> t_test = np.linspace(t_min, t_max, n_points * 10)
729 >>> y_true = gen_data(t_test, a, b, c)
730 >>> y_lsq = gen_data(t_test, *res_lsq.x)
731 >>> y_soft_l1 = gen_data(t_test, *res_soft_l1.x)
732 >>> y_log = gen_data(t_test, *res_log.x)
733 ...
734 >>> import matplotlib.pyplot as plt
735 >>> plt.plot(t_train, y_train, 'o')
736 >>> plt.plot(t_test, y_true, 'k', linewidth=2, label='true')
737 >>> plt.plot(t_test, y_lsq, label='linear loss')
738 >>> plt.plot(t_test, y_soft_l1, label='soft_l1 loss')
739 >>> plt.plot(t_test, y_log, label='cauchy loss')
740 >>> plt.xlabel("t")
741 >>> plt.ylabel("y")
742 >>> plt.legend()
743 >>> plt.show()
745 In the next example, we show how complex-valued residual functions of
746 complex variables can be optimized with ``least_squares()``. Consider the
747 following function:
749 >>> def f(z):
750 ... return z - (0.5 + 0.5j)
752 We wrap it into a function of real variables that returns real residuals
753 by simply handling the real and imaginary parts as independent variables:
755 >>> def f_wrap(x):
756 ... fx = f(x[0] + 1j*x[1])
757 ... return np.array([fx.real, fx.imag])
759 Thus, instead of the original m-D complex function of n complex
760 variables we optimize a 2m-D real function of 2n real variables:
762 >>> from scipy.optimize import least_squares
763 >>> res_wrapped = least_squares(f_wrap, (0.1, 0.1), bounds=([0, 0], [1, 1]))
764 >>> z = res_wrapped.x[0] + res_wrapped.x[1]*1j
765 >>> z
766 (0.49999999999925893+0.49999999999925893j)
768 """
769 if method not in ['trf', 'dogbox', 'lm']:
770 raise ValueError("`method` must be 'trf', 'dogbox' or 'lm'.")
772 if jac not in ['2-point', '3-point', 'cs'] and not callable(jac):
773 raise ValueError("`jac` must be '2-point', '3-point', 'cs' or "
774 "callable.")
776 if tr_solver not in [None, 'exact', 'lsmr']:
777 raise ValueError("`tr_solver` must be None, 'exact' or 'lsmr'.")
779 if loss not in IMPLEMENTED_LOSSES and not callable(loss):
780 raise ValueError("`loss` must be one of {} or a callable."
781 .format(IMPLEMENTED_LOSSES.keys()))
783 if method == 'lm' and loss != 'linear':
784 raise ValueError("method='lm' supports only 'linear' loss function.")
786 if verbose not in [0, 1, 2]:
787 raise ValueError("`verbose` must be in [0, 1, 2].")
789 if max_nfev is not None and max_nfev <= 0:
790 raise ValueError("`max_nfev` must be None or positive integer.")
792 if np.iscomplexobj(x0):
793 raise ValueError("`x0` must be real.")
795 x0 = np.atleast_1d(x0).astype(float)
797 if x0.ndim > 1:
798 raise ValueError("`x0` must have at most 1 dimension.")
800 if isinstance(bounds, Bounds):
801 lb, ub = bounds.lb, bounds.ub
802 bounds = (lb, ub)
803 else:
804 if len(bounds) == 2:
805 lb, ub = prepare_bounds(bounds, x0.shape[0])
806 else:
807 raise ValueError("`bounds` must contain 2 elements.")
809 if method == 'lm' and not np.all((lb == -np.inf) & (ub == np.inf)):
810 raise ValueError("Method 'lm' doesn't support bounds.")
812 if lb.shape != x0.shape or ub.shape != x0.shape:
813 raise ValueError("Inconsistent shapes between bounds and `x0`.")
815 if np.any(lb >= ub):
816 raise ValueError("Each lower bound must be strictly less than each "
817 "upper bound.")
819 if not in_bounds(x0, lb, ub):
820 raise ValueError("`x0` is infeasible.")
822 x_scale = check_x_scale(x_scale, x0)
824 ftol, xtol, gtol = check_tolerance(ftol, xtol, gtol, method)
826 if method == 'trf':
827 x0 = make_strictly_feasible(x0, lb, ub)
829 def fun_wrapped(x):
830 return np.atleast_1d(fun(x, *args, **kwargs))
832 f0 = fun_wrapped(x0)
834 if f0.ndim != 1:
835 raise ValueError("`fun` must return at most 1-d array_like. "
836 "f0.shape: {}".format(f0.shape))
838 if not np.all(np.isfinite(f0)):
839 raise ValueError("Residuals are not finite in the initial point.")
841 n = x0.size
842 m = f0.size
844 if method == 'lm' and m < n:
845 raise ValueError("Method 'lm' doesn't work when the number of "
846 "residuals is less than the number of variables.")
848 loss_function = construct_loss_function(m, loss, f_scale)
849 if callable(loss):
850 rho = loss_function(f0)
851 if rho.shape != (3, m):
852 raise ValueError("The return value of `loss` callable has wrong "
853 "shape.")
854 initial_cost = 0.5 * np.sum(rho[0])
855 elif loss_function is not None:
856 initial_cost = loss_function(f0, cost_only=True)
857 else:
858 initial_cost = 0.5 * np.dot(f0, f0)
860 if callable(jac):
861 J0 = jac(x0, *args, **kwargs)
863 if issparse(J0):
864 J0 = J0.tocsr()
866 def jac_wrapped(x, _=None):
867 return jac(x, *args, **kwargs).tocsr()
869 elif isinstance(J0, LinearOperator):
870 def jac_wrapped(x, _=None):
871 return jac(x, *args, **kwargs)
873 else:
874 J0 = np.atleast_2d(J0)
876 def jac_wrapped(x, _=None):
877 return np.atleast_2d(jac(x, *args, **kwargs))
879 else: # Estimate Jacobian by finite differences.
880 if method == 'lm':
881 if jac_sparsity is not None:
882 raise ValueError("method='lm' does not support "
883 "`jac_sparsity`.")
885 if jac != '2-point':
886 warn("jac='{}' works equivalently to '2-point' "
887 "for method='lm'.".format(jac))
889 J0 = jac_wrapped = None
890 else:
891 if jac_sparsity is not None and tr_solver == 'exact':
892 raise ValueError("tr_solver='exact' is incompatible "
893 "with `jac_sparsity`.")
895 jac_sparsity = check_jac_sparsity(jac_sparsity, m, n)
897 def jac_wrapped(x, f):
898 J = approx_derivative(fun, x, rel_step=diff_step, method=jac,
899 f0=f, bounds=bounds, args=args,
900 kwargs=kwargs, sparsity=jac_sparsity)
901 if J.ndim != 2: # J is guaranteed not sparse.
902 J = np.atleast_2d(J)
904 return J
906 J0 = jac_wrapped(x0, f0)
908 if J0 is not None:
909 if J0.shape != (m, n):
910 raise ValueError(
911 "The return value of `jac` has wrong shape: expected {}, "
912 "actual {}.".format((m, n), J0.shape))
914 if not isinstance(J0, np.ndarray):
915 if method == 'lm':
916 raise ValueError("method='lm' works only with dense "
917 "Jacobian matrices.")
919 if tr_solver == 'exact':
920 raise ValueError(
921 "tr_solver='exact' works only with dense "
922 "Jacobian matrices.")
924 jac_scale = isinstance(x_scale, str) and x_scale == 'jac'
925 if isinstance(J0, LinearOperator) and jac_scale:
926 raise ValueError("x_scale='jac' can't be used when `jac` "
927 "returns LinearOperator.")
929 if tr_solver is None:
930 if isinstance(J0, np.ndarray):
931 tr_solver = 'exact'
932 else:
933 tr_solver = 'lsmr'
935 if method == 'lm':
936 result = call_minpack(fun_wrapped, x0, jac_wrapped, ftol, xtol, gtol,
937 max_nfev, x_scale, diff_step)
939 elif method == 'trf':
940 result = trf(fun_wrapped, jac_wrapped, x0, f0, J0, lb, ub, ftol, xtol,
941 gtol, max_nfev, x_scale, loss_function, tr_solver,
942 tr_options.copy(), verbose)
944 elif method == 'dogbox':
945 if tr_solver == 'lsmr' and 'regularize' in tr_options:
946 warn("The keyword 'regularize' in `tr_options` is not relevant "
947 "for 'dogbox' method.")
948 tr_options = tr_options.copy()
949 del tr_options['regularize']
951 result = dogbox(fun_wrapped, jac_wrapped, x0, f0, J0, lb, ub, ftol,
952 xtol, gtol, max_nfev, x_scale, loss_function,
953 tr_solver, tr_options, verbose)
955 result.message = TERMINATION_MESSAGES[result.status]
956 result.success = result.status > 0
958 if verbose >= 1:
959 print(result.message)
960 print("Function evaluations {}, initial cost {:.4e}, final cost "
961 "{:.4e}, first-order optimality {:.2e}."
962 .format(result.nfev, initial_cost, result.cost,
963 result.optimality))
965 return result