dlnpyutils package

Submodules

dlnpyutils.astro module

dlnpyutils.bindata module

dlnpyutils.coords module

dlnpyutils.db module

dlnpyutils.galaxy_model module

dlnpyutils.gaps module

dlnpyutils.gaps.gap(data, refs=None, nrefs=20, ks=range(1, 11))

Compute the Gap statistic for an nxm dataset in data. Either give a precomputed set of reference distributions in refs as an (n,m,k) scipy array, or state the number k of reference distributions in nrefs for automatic generation with a uniformed distribution within the bounding box of data. Give the list of k-values for which you want to compute the statistic in ks.

dlnpyutils.job_daemon module

dlnpyutils.ladfit module

dlnpyutils.ladfit.ladfit(xin, yin)

Copyright (c) 1994-2015, Exelis Visual Information Solutions, Inc. All rights reserved. Unauthorized reproduction is prohibited.

LADFIT

This function fits the paired data {X(i), Y(i)} to the linear model, y = A + Bx, using a “robust” least absolute deviation method. The result is a two-element vector containing the model parameters, A and B.

Result = LADFIT(X, Y)

Parameters

X: An n-element vector of type integer, float or double.

Y: An n-element vector of type integer, float or double.

EXAMPLE:

Define two n-element vectors of paired data.

x = [-3.20, 4.49, -1.66, 0.64, -2.43, -0.89, -0.12, 1.41, $

2.95, 2.18, 3.72, 5.26]

y = [-7.14, -1.30, -4.26, -1.90, -6.19, -3.98, -2.87, -1.66, $

-0.78, -2.61, 0.31, 1.74]

Compute the model parameters, A and B.

result = ladfit(x, y, absdev = absdev)

The result should be the two-element vector:

[-3.15301, 0.930440]

The keyword parameter should be returned as:

absdev = 0.636851

REFERENCE:

Numerical Recipes, The Art of Scientific Computing (Second Edition) Cambridge University Press, 2nd Edition. ISBN 0-521-43108-5

This is adapted from the routine MEDFIT described in:

Fitting a Line by Minimizing Absolute Deviation, Page 703.

MODIFICATION HISTORY:

Written by: GGS, RSI, September 1994 Modified: GGS, RSI, July 1995

Corrected an infinite loop condition that occured when the X input parameter contained mostly negative data.

Modified: GGS, RSI, October 1996

If least-absolute-deviation convergence condition is not satisfied, the algorithm switches to a chi-squared model. Modified keyword checking and use of double precision.

Modified: GGS, RSI, November 1996

Fixed an error in the computation of the median with even-length input data. See EVEN keyword to MEDIAN.

Modified: DMS, RSI, June 1997

Simplified logic, remove SIGN and MDfunc2 functions.

Modified: RJF, RSI, Jan 1999

Fixed the variance computation by adding some double conversions. This prevents the function from generating NaNs on some specific datasets (bug 11680).

Modified: CT, RSI, July 2002: Convert inputs to float or double.

Change constants to double precision if necessary.

CT, March 2004: Check for quick return if we found solution.

dlnpyutils.ladfit.ladmdfunc(b, x, y, eps=1e-07)

dlnpyutils.least_squares module

Generic interface for least-square minimization.

dlnpyutils.least_squares.arctan(z, rho, cost_only)

dlnpyutils.least_squares.call_minpack(fun, x0, jac, ftol, xtol, gtol, max_nfev, x_scale, diff_step)

dlnpyutils.least_squares.cauchy(z, rho, cost_only)

dlnpyutils.least_squares.check_jac_sparsity(jac_sparsity, m, n)

dlnpyutils.least_squares.check_tolerance(ftol, xtol, gtol)

dlnpyutils.least_squares.check_x_scale(x_scale, x0)

dlnpyutils.least_squares.construct_loss_function(m, loss, f_scale)

dlnpyutils.least_squares.huber(z, rho, cost_only)

dlnpyutils.least_squares.least_squares(fun, x0, jac='2-point', bounds=(-inf, inf), method='trf', ftol=1e-08, xtol=1e-08, gtol=1e-08, x_scale=1.0, loss='linear', f_scale=1.0, diff_step=None, tr_solver=None, tr_options={}, jac_sparsity=None, max_nfev=None, verbose=0, dx_lim=None, args=(), kwargs={})

Solve a nonlinear least-squares problem with bounds on the variables.

Given the residuals f(x) (an m-dimensional real function of n real variables) and the loss function rho(s) (a scalar function), least_squares finds a local minimum of the cost function F(x):

minimize F(x) = 0.5 * sum(rho(f_i(x)**2), i = 0, ..., m - 1)
subject to lb <= x <= ub

The purpose of the loss function rho(s) is to reduce the influence of outliers on the solution.

Parameters

funcallable

Function which computes the vector of residuals, with the signature fun(x, *args, **kwargs), i.e., the minimization proceeds with respect to its first argument. The argument x passed to this function is an ndarray of shape (n,) (never a scalar, even for n=1). It must return a 1-d array_like of shape (m,) or a scalar. If the argument x is complex or the function fun returns complex residuals, it must be wrapped in a real function of real arguments, as shown at the end of the Examples section.

x0array_like with shape (n,) or float

Initial guess on independent variables. If float, it will be treated as a 1-d array with one element.

jac{‘2-point’, ‘3-point’, ‘cs’, callable}, optional

Method of computing the Jacobian matrix (an m-by-n matrix, where element (i, j) is the partial derivative of f[i] with respect to x[j]). The keywords select a finite difference scheme for numerical estimation. The scheme ‘3-point’ is more accurate, but requires twice as many operations as ‘2-point’ (default). The scheme ‘cs’ uses complex steps, and while potentially the most accurate, it is applicable only when fun correctly handles complex inputs and can be analytically continued to the complex plane. Method ‘lm’ always uses the ‘2-point’ scheme. If callable, it is used as jac(x, *args, **kwargs) and should return a good approximation (or the exact value) for the Jacobian as an array_like (np.atleast_2d is applied), a sparse matrix or a scipy.sparse.linalg.LinearOperator.

bounds2-tuple of array_like, optional

Lower and upper bounds on independent variables. Defaults to no bounds. Each array must match the size of x0 or be a scalar, in the latter case a bound will be the same for all variables. Use np.inf with an appropriate sign to disable bounds on all or some variables.

method{‘trf’, ‘dogbox’, ‘lm’}, optional

Algorithm to perform minimization.

• ‘trf’ : Trust Region Reflective algorithm, particularly suitable for large sparse problems with bounds. Generally robust method.

• ‘dogbox’ : dogleg algorithm with rectangular trust regions, typical use case is small problems with bounds. Not recommended for problems with rank-deficient Jacobian.

• ‘lm’ : Levenberg-Marquardt algorithm as implemented in MINPACK. Doesn’t handle bounds and sparse Jacobians. Usually the most efficient method for small unconstrained problems.

Default is ‘trf’. See Notes for more information.

ftolfloat or None, optional

Tolerance for termination by the change of the cost function. Default is 1e-8. The optimization process is stopped when dF < ftol * F, and there was an adequate agreement between a local quadratic model and the true model in the last step. If None, the termination by this condition is disabled.

xtolfloat or None, optional

Tolerance for termination by the change of the independent variables. Default is 1e-8. The exact condition depends on the method used:

• For ‘trf’ and ‘dogbox’ : norm(dx) < xtol * (xtol + norm(x))

• For ‘lm’ : Delta < xtol * norm(xs), where Delta is a trust-region radius and xs is the value of x scaled according to x_scale parameter (see below).

If None, the termination by this condition is disabled.

gtolfloat or None, optional

Tolerance for termination by the norm of the gradient. Default is 1e-8. The exact condition depends on a method used:

• For ‘trf’ : norm(g_scaled, ord=np.inf) < gtol, where g_scaled is the value of the gradient scaled to account for the presence of the bounds [STIR].

• For ‘dogbox’ : norm(g_free, ord=np.inf) < gtol, where g_free is the gradient with respect to the variables which are not in the optimal state on the boundary.

• For ‘lm’ : the maximum absolute value of the cosine of angles between columns of the Jacobian and the residual vector is less than gtol, or the residual vector is zero.

If None, the termination by this condition is disabled.

x_scalearray_like or ‘jac’, optional

Characteristic scale of each variable. Setting x_scale is equivalent to reformulating the problem in scaled variables xs = x / x_scale. An alternative view is that the size of a trust region along j-th dimension is proportional to x_scale[j]. Improved convergence may be achieved by setting x_scale such that a step of a given size along any of the scaled variables has a similar effect on the cost function. If set to ‘jac’, the scale is iteratively updated using the inverse norms of the columns of the Jacobian matrix (as described in [JJMore]).

lossstr or callable, optional

Determines the loss function. The following keyword values are allowed:

• ‘linear’ (default) : rho(z) = z. Gives a standard least-squares problem.

• ‘soft_l1’ : rho(z) = 2 * ((1 + z)**0.5 - 1). The smooth approximation of l1 (absolute value) loss. Usually a good choice for robust least squares.

• ‘huber’ : rho(z) = z if z <= 1 else 2*z**0.5 - 1. Works similarly to ‘soft_l1’.

• ‘cauchy’ : rho(z) = ln(1 + z). Severely weakens outliers influence, but may cause difficulties in optimization process.

• ‘arctan’ : rho(z) = arctan(z). Limits a maximum loss on a single residual, has properties similar to ‘cauchy’.

If callable, it must take a 1-d ndarray z=f**2 and return an array_like with shape (3, m) where row 0 contains function values, row 1 contains first derivatives and row 2 contains second derivatives. Method ‘lm’ supports only ‘linear’ loss.

f_scalefloat, optional

Value of soft margin between inlier and outlier residuals, default is 1.0. The loss function is evaluated as follows rho_(f**2) = C**2 * rho(f**2 / C**2), where C is f_scale, and rho is determined by loss parameter. This parameter has no effect with loss='linear', but for other loss values it is of crucial importance.

max_nfevNone or int, optional

Maximum number of function evaluations before the termination. If None (default), the value is chosen automatically:

• For ‘trf’ and ‘dogbox’ : 100 * n.

• For ‘lm’ : 100 * n if jac is callable and 100 * n * (n + 1) otherwise (because ‘lm’ counts function calls in Jacobian estimation).

diff_stepNone or array_like, optional

Determines the relative step size for the finite difference approximation of the Jacobian. The actual step is computed as x * diff_step. If None (default), then diff_step is taken to be a conventional “optimal” power of machine epsilon for the finite difference scheme used [NR].

tr_solver{None, ‘exact’, ‘lsmr’}, optional

Method for solving trust-region subproblems, relevant only for ‘trf’ and ‘dogbox’ methods.

• ‘exact’ is suitable for not very large problems with dense Jacobian matrices. The computational complexity per iteration is comparable to a singular value decomposition of the Jacobian matrix.

• ‘lsmr’ is suitable for problems with sparse and large Jacobian matrices. It uses the iterative procedure scipy.sparse.linalg.lsmr for finding a solution of a linear least-squares problem and only requires matrix-vector product evaluations.

If None (default) the solver is chosen based on the type of Jacobian returned on the first iteration.

tr_optionsdict, optional

Keyword options passed to trust-region solver.

• tr_solver='exact': tr_options are ignored.

• tr_solver='lsmr': options for scipy.sparse.linalg.lsmr. Additionally method='trf' supports ‘regularize’ option (bool, default is True) which adds a regularization term to the normal equation, which improves convergence if the Jacobian is rank-deficient [Byrd] (eq. 3.4).

jac_sparsity{None, array_like, sparse matrix}, optional

Defines the sparsity structure of the Jacobian matrix for finite difference estimation, its shape must be (m, n). If the Jacobian has only few non-zero elements in each row, providing the sparsity structure will greatly speed up the computations [Curtis]. A zero entry means that a corresponding element in the Jacobian is identically zero. If provided, forces the use of ‘lsmr’ trust-region solver. If None (default) then dense differencing will be used. Has no effect for ‘lm’ method.

verbose{0, 1, 2}, optional

Level of algorithm’s verbosity:

• 0 (default) : work silently.

• 1 : display a termination report.

• 2 : display progress during iterations (not supported by ‘lm’ method).

args, kwargstuple and dict, optional

Additional arguments passed to fun and jac. Both empty by default. The calling signature is fun(x, *args, **kwargs) and the same for jac.

Returns

OptimizeResult with the following fields defined:

xndarray, shape (n,)

Solution found.

costfloat

Value of the cost function at the solution.

funndarray, shape (m,)

Vector of residuals at the solution.

jacndarray, sparse matrix or LinearOperator, shape (m, n)

Modified Jacobian matrix at the solution, in the sense that J^T J is a Gauss-Newton approximation of the Hessian of the cost function. The type is the same as the one used by the algorithm.

gradndarray, shape (m,)

Gradient of the cost function at the solution.

optimalityfloat

First-order optimality measure. In unconstrained problems, it is always the uniform norm of the gradient. In constrained problems, it is the quantity which was compared with gtol during iterations.

active_maskndarray of int, shape (n,)

Each component shows whether a corresponding constraint is active (that is, whether a variable is at the bound):

• 0 : a constraint is not active.

• -1 : a lower bound is active.

• 1 : an upper bound is active.

Might be somewhat arbitrary for ‘trf’ method as it generates a sequence of strictly feasible iterates and active_mask is determined within a tolerance threshold.

nfevint

Number of function evaluations done. Methods ‘trf’ and ‘dogbox’ do not count function calls for numerical Jacobian approximation, as opposed to ‘lm’ method.

njevint or None

Number of Jacobian evaluations done. If numerical Jacobian approximation is used in ‘lm’ method, it is set to None.

statusint

The reason for algorithm termination:

• -1 : improper input parameters status returned from MINPACK.

• 0 : the maximum number of function evaluations is exceeded.

• 1 : gtol termination condition is satisfied.

• 2 : ftol termination condition is satisfied.

• 3 : xtol termination condition is satisfied.

• 4 : Both ftol and xtol termination conditions are satisfied.

messagestr

Verbal description of the termination reason.

successbool

True if one of the convergence criteria is satisfied (status > 0).

See also

leastsq

A legacy wrapper for the MINPACK implementation of the Levenberg-Marquadt algorithm.

curve_fit

Least-squares minimization applied to a curve fitting problem.

Notes

Method ‘lm’ (Levenberg-Marquardt) calls a wrapper over least-squares algorithms implemented in MINPACK (lmder, lmdif). It runs the Levenberg-Marquardt algorithm formulated as a trust-region type algorithm. The implementation is based on paper [JJMore], it is very robust and efficient with a lot of smart tricks. It should be your first choice for unconstrained problems. Note that it doesn’t support bounds. Also it doesn’t work when m < n.

Method ‘trf’ (Trust Region Reflective) is motivated by the process of solving a system of equations, which constitute the first-order optimality condition for a bound-constrained minimization problem as formulated in [STIR]. The algorithm iteratively solves trust-region subproblems augmented by a special diagonal quadratic term and with trust-region shape determined by the distance from the bounds and the direction of the gradient. This enhancements help to avoid making steps directly into bounds and efficiently explore the whole space of variables. To further improve convergence, the algorithm considers search directions reflected from the bounds. To obey theoretical requirements, the algorithm keeps iterates strictly feasible. With dense Jacobians trust-region subproblems are solved by an exact method very similar to the one described in [JJMore] (and implemented in MINPACK). The difference from the MINPACK implementation is that a singular value decomposition of a Jacobian matrix is done once per iteration, instead of a QR decomposition and series of Givens rotation eliminations. For large sparse Jacobians a 2-d subspace approach of solving trust-region subproblems is used [STIR], [Byrd]. The subspace is spanned by a scaled gradient and an approximate Gauss-Newton solution delivered by scipy.sparse.linalg.lsmr. When no constraints are imposed the algorithm is very similar to MINPACK and has generally comparable performance. The algorithm works quite robust in unbounded and bounded problems, thus it is chosen as a default algorithm.

Method ‘dogbox’ operates in a trust-region framework, but considers rectangular trust regions as opposed to conventional ellipsoids [Voglis]. The intersection of a current trust region and initial bounds is again rectangular, so on each iteration a quadratic minimization problem subject to bound constraints is solved approximately by Powell’s dogleg method [NumOpt]. The required Gauss-Newton step can be computed exactly for dense Jacobians or approximately by scipy.sparse.linalg.lsmr for large sparse Jacobians. The algorithm is likely to exhibit slow convergence when the rank of Jacobian is less than the number of variables. The algorithm often outperforms ‘trf’ in bounded problems with a small number of variables.

Robust loss functions are implemented as described in [BA]. The idea is to modify a residual vector and a Jacobian matrix on each iteration such that computed gradient and Gauss-Newton Hessian approximation match the true gradient and Hessian approximation of the cost function. Then the algorithm proceeds in a normal way, i.e. robust loss functions are implemented as a simple wrapper over standard least-squares algorithms.

New in version 0.17.0.

References

STIR(1,2,3)

M. A. Branch, T. F. Coleman, and Y. Li, “A Subspace, Interior, and Conjugate Gradient Method for Large-Scale Bound-Constrained Minimization Problems,” SIAM Journal on Scientific Computing, Vol. 21, Number 1, pp 1-23, 1999.

NR

William H. Press et. al., “Numerical Recipes. The Art of Scientific Computing. 3rd edition”, Sec. 5.7.

Byrd(1,2)

R. H. Byrd, R. B. Schnabel and G. A. Shultz, “Approximate solution of the trust region problem by minimization over two-dimensional subspaces”, Math. Programming, 40, pp. 247-263, 1988.

Curtis

A. Curtis, M. J. D. Powell, and J. Reid, “On the estimation of sparse Jacobian matrices”, Journal of the Institute of Mathematics and its Applications, 13, pp. 117-120, 1974.

JJMore(1,2,3)

J. J. More, “The Levenberg-Marquardt Algorithm: Implementation and Theory,” Numerical Analysis, ed. G. A. Watson, Lecture Notes in Mathematics 630, Springer Verlag, pp. 105-116, 1977.

Voglis

C. Voglis and I. E. Lagaris, “A Rectangular Trust Region Dogleg Approach for Unconstrained and Bound Constrained Nonlinear Optimization”, WSEAS International Conference on Applied Mathematics, Corfu, Greece, 2004.

NumOpt

J. Nocedal and S. J. Wright, “Numerical optimization, 2nd edition”, Chapter 4.

BA

B. Triggs et. al., “Bundle Adjustment - A Modern Synthesis”, Proceedings of the International Workshop on Vision Algorithms: Theory and Practice, pp. 298-372, 1999.

Examples

In this example we find a minimum of the Rosenbrock function without bounds on independent variables.

>>> def fun_rosenbrock(x):
...     return np.array([10 * (x[1] - x[0]**2), (1 - x[0])])

Notice that we only provide the vector of the residuals. The algorithm constructs the cost function as a sum of squares of the residuals, which gives the Rosenbrock function. The exact minimum is at x = [1.0, 1.0].

>>> from scipy.optimize import least_squares
>>> x0_rosenbrock = np.array([2, 2])
>>> res_1 = least_squares(fun_rosenbrock, x0_rosenbrock)
>>> res_1.x
array([ 1.,  1.])
>>> res_1.cost
9.8669242910846867e-30
>>> res_1.optimality
8.8928864934219529e-14

We now constrain the variables, in such a way that the previous solution becomes infeasible. Specifically, we require that x[1] >= 1.5, and x[0] left unconstrained. To this end, we specify the bounds parameter to least_squares in the form bounds=([-np.inf, 1.5], np.inf).

We also provide the analytic Jacobian:

>>> def jac_rosenbrock(x):
...     return np.array([
...         [-20 * x[0], 10],
...         [-1, 0]])

Putting this all together, we see that the new solution lies on the bound:

>>> res_2 = least_squares(fun_rosenbrock, x0_rosenbrock, jac_rosenbrock,
...                       bounds=([-np.inf, 1.5], np.inf))
>>> res_2.x
array([ 1.22437075,  1.5       ])
>>> res_2.cost
0.025213093946805685
>>> res_2.optimality
1.5885401433157753e-07

Now we solve a system of equations (i.e., the cost function should be zero at a minimum) for a Broyden tridiagonal vector-valued function of 100000 variables:

>>> def fun_broyden(x):
...     f = (3 - x) * x + 1
...     f[1:] -= x[:-1]
...     f[:-1] -= 2 * x[1:]
...     return f

The corresponding Jacobian matrix is sparse. We tell the algorithm to estimate it by finite differences and provide the sparsity structure of Jacobian to significantly speed up this process.

>>> from scipy.sparse import lil_matrix
>>> def sparsity_broyden(n):
...     sparsity = lil_matrix((n, n), dtype=int)
...     i = np.arange(n)
...     sparsity[i, i] = 1
...     i = np.arange(1, n)
...     sparsity[i, i - 1] = 1
...     i = np.arange(n - 1)
...     sparsity[i, i + 1] = 1
...     return sparsity
...
>>> n = 100000
>>> x0_broyden = -np.ones(n)
...
>>> res_3 = least_squares(fun_broyden, x0_broyden,
...                       jac_sparsity=sparsity_broyden(n))
>>> res_3.cost
4.5687069299604613e-23
>>> res_3.optimality
1.1650454296851518e-11

Let’s also solve a curve fitting problem using robust loss function to take care of outliers in the data. Define the model function as y = a + b * exp(c * t), where t is a predictor variable, y is an observation and a, b, c are parameters to estimate.

First, define the function which generates the data with noise and outliers, define the model parameters, and generate data:

>>> def gen_data(t, a, b, c, noise=0, n_outliers=0, random_state=0):
...     y = a + b * np.exp(t * c)
...
...     rnd = np.random.RandomState(random_state)
...     error = noise * rnd.randn(t.size)
...     outliers = rnd.randint(0, t.size, n_outliers)
...     error[outliers] *= 10
...
...     return y + error
...
>>> a = 0.5
>>> b = 2.0
>>> c = -1
>>> t_min = 0
>>> t_max = 10
>>> n_points = 15
...
>>> t_train = np.linspace(t_min, t_max, n_points)
>>> y_train = gen_data(t_train, a, b, c, noise=0.1, n_outliers=3)

Define function for computing residuals and initial estimate of parameters.

>>> def fun(x, t, y):
...     return x[0] + x[1] * np.exp(x[2] * t) - y
...
>>> x0 = np.array([1.0, 1.0, 0.0])

Compute a standard least-squares solution:

>>> res_lsq = least_squares(fun, x0, args=(t_train, y_train))

Now compute two solutions with two different robust loss functions. The parameter f_scale is set to 0.1, meaning that inlier residuals should not significantly exceed 0.1 (the noise level used).

>>> res_soft_l1 = least_squares(fun, x0, loss='soft_l1', f_scale=0.1,
...                             args=(t_train, y_train))
>>> res_log = least_squares(fun, x0, loss='cauchy', f_scale=0.1,
...                         args=(t_train, y_train))

And finally plot all the curves. We see that by selecting an appropriate loss we can get estimates close to optimal even in the presence of strong outliers. But keep in mind that generally it is recommended to try ‘soft_l1’ or ‘huber’ losses first (if at all necessary) as the other two options may cause difficulties in optimization process.

>>> t_test = np.linspace(t_min, t_max, n_points * 10)
>>> y_true = gen_data(t_test, a, b, c)
>>> y_lsq = gen_data(t_test, *res_lsq.x)
>>> y_soft_l1 = gen_data(t_test, *res_soft_l1.x)
>>> y_log = gen_data(t_test, *res_log.x)
...
>>> import matplotlib.pyplot as plt
>>> plt.plot(t_train, y_train, 'o')
>>> plt.plot(t_test, y_true, 'k', linewidth=2, label='true')
>>> plt.plot(t_test, y_lsq, label='linear loss')
>>> plt.plot(t_test, y_soft_l1, label='soft_l1 loss')
>>> plt.plot(t_test, y_log, label='cauchy loss')
>>> plt.xlabel("t")
>>> plt.ylabel("y")
>>> plt.legend()
>>> plt.show()

In the next example, we show how complex-valued residual functions of complex variables can be optimized with least_squares(). Consider the following function:

>>> def f(z):
...     return z - (0.5 + 0.5j)

We wrap it into a function of real variables that returns real residuals by simply handling the real and imaginary parts as independent variables:

>>> def f_wrap(x):
...     fx = f(x[0] + 1j*x[1])
...     return np.array([fx.real, fx.imag])

Thus, instead of the original m-dimensional complex function of n complex variables we optimize a 2m-dimensional real function of 2n real variables:

>>> from scipy.optimize import least_squares
>>> res_wrapped = least_squares(f_wrap, (0.1, 0.1), bounds=([0, 0], [1, 1]))
>>> z = res_wrapped.x[0] + res_wrapped.x[1]*1j
>>> z
(0.49999999999925893+0.49999999999925893j)

dlnpyutils.least_squares.prepare_bounds(bounds, n)

dlnpyutils.least_squares.soft_l1(z, rho, cost_only)

dlnpyutils.minpack module

dlnpyutils.minpack.curve_fit(f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False, check_finite=True, bounds=(-inf, inf), method=None, jac=None, **kwargs)

Use non-linear least squares to fit a function, f, to data.

Assumes ydata = f(xdata, *params) + eps

Parameters

fcallable

The model function, f(x, …). It must take the independent variable as the first argument and the parameters to fit as separate remaining arguments.

xdataarray_like or object

The independent variable where the data is measured. Should usually be an M-length sequence or an (k,M)-shaped array for functions with k predictors, but can actually be any object.

ydataarray_like

The dependent data, a length M array - nominally f(xdata, ...).

p0array_like, optional

Initial guess for the parameters (length N). If None, then the initial values will all be 1 (if the number of parameters for the function can be determined using introspection, otherwise a ValueError is raised).

sigmaNone or M-length sequence or MxM array, optional

Determines the uncertainty in ydata. If we define residuals as r = ydata - f(xdata, *popt), then the interpretation of sigma depends on its number of dimensions:

• A 1-d sigma should contain values of standard deviations of errors in ydata. In this case, the optimized function is chisq = sum((r / sigma) ** 2).

• A 2-d sigma should contain the covariance matrix of errors in ydata. In this case, the optimized function is chisq = r.T @ inv(sigma) @ r.

  New in version 0.19.

None (default) is equivalent of 1-d sigma filled with ones.

absolute_sigmabool, optional

If True, sigma is used in an absolute sense and the estimated parameter covariance pcov reflects these absolute values.

If False, only the relative magnitudes of the sigma values matter. The returned parameter covariance matrix pcov is based on scaling sigma by a constant factor. This constant is set by demanding that the reduced chisq for the optimal parameters popt when using the scaled sigma equals unity. In other words, sigma is scaled to match the sample variance of the residuals after the fit. Mathematically, pcov(absolute_sigma=False) = pcov(absolute_sigma=True) * chisq(popt)/(M-N)

check_finitebool, optional

If True, check that the input arrays do not contain nans of infs, and raise a ValueError if they do. Setting this parameter to False may silently produce nonsensical results if the input arrays do contain nans. Default is True.

bounds2-tuple of array_like, optional

Lower and upper bounds on parameters. Defaults to no bounds. Each element of the tuple must be either an array with the length equal to the number of parameters, or a scalar (in which case the bound is taken to be the same for all parameters.) Use np.inf with an appropriate sign to disable bounds on all or some parameters.

New in version 0.17.

method{‘lm’, ‘trf’, ‘dogbox’}, optional

Method to use for optimization. See least_squares for more details. Default is ‘lm’ for unconstrained problems and ‘trf’ if bounds are provided. The method ‘lm’ won’t work when the number of observations is less than the number of variables, use ‘trf’ or ‘dogbox’ in this case.

New in version 0.17.

jaccallable, string or None, optional

Function with signature jac(x, ...) which computes the Jacobian matrix of the model function with respect to parameters as a dense array_like structure. It will be scaled according to provided sigma. If None (default), the Jacobian will be estimated numerically. String keywords for ‘trf’ and ‘dogbox’ methods can be used to select a finite difference scheme, see least_squares.

New in version 0.18.

kwargs

Keyword arguments passed to leastsq for method='lm' or least_squares otherwise.

Returns

poptarray

Optimal values for the parameters so that the sum of the squared residuals of f(xdata, *popt) - ydata is minimized

pcov2d array

The estimated covariance of popt. The diagonals provide the variance of the parameter estimate. To compute one standard deviation errors on the parameters use perr = np.sqrt(np.diag(pcov)).

How the sigma parameter affects the estimated covariance depends on absolute_sigma argument, as described above.

If the Jacobian matrix at the solution doesn’t have a full rank, then ‘lm’ method returns a matrix filled with np.inf, on the other hand ‘trf’ and ‘dogbox’ methods use Moore-Penrose pseudoinverse to compute the covariance matrix.

Raises

ValueError

if either ydata or xdata contain NaNs, or if incompatible options are used.

RuntimeError

if the least-squares minimization fails.

OptimizeWarning

if covariance of the parameters can not be estimated.

See also

least_squares

Minimize the sum of squares of nonlinear functions.

scipy.stats.linregress

Calculate a linear least squares regression for two sets of measurements.

Notes

With method='lm', the algorithm uses the Levenberg-Marquardt algorithm through leastsq. Note that this algorithm can only deal with unconstrained problems.

Box constraints can be handled by methods ‘trf’ and ‘dogbox’. Refer to the docstring of least_squares for more information.

Examples

>>> import matplotlib.pyplot as plt
>>> from scipy.optimize import curve_fit

>>> def func(x, a, b, c):
...     return a * np.exp(-b * x) + c

Define the data to be fit with some noise:

>>> xdata = np.linspace(0, 4, 50)
>>> y = func(xdata, 2.5, 1.3, 0.5)
>>> np.random.seed(1729)
>>> y_noise = 0.2 * np.random.normal(size=xdata.size)
>>> ydata = y + y_noise
>>> plt.plot(xdata, ydata, 'b-', label='data')

Fit for the parameters a, b, c of the function func:

>>> popt, pcov = curve_fit(func, xdata, ydata)
>>> popt
array([ 2.55423706,  1.35190947,  0.47450618])
>>> plt.plot(xdata, func(xdata, *popt), 'r-',
...          label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt))

Constrain the optimization to the region of 0 <= a <= 3, 0 <= b <= 1 and 0 <= c <= 0.5:

>>> popt, pcov = curve_fit(func, xdata, ydata, bounds=(0, [3., 1., 0.5]))
>>> popt
array([ 2.43708906,  1.        ,  0.35015434])
>>> plt.plot(xdata, func(xdata, *popt), 'g--',
...          label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt))

>>> plt.xlabel('x')
>>> plt.ylabel('y')
>>> plt.legend()
>>> plt.show()

dlnpyutils.minpack.fixed_point(func, x0, args=(), xtol=1e-08, maxiter=500, method='del2')

Find a fixed point of the function.

Given a function of one or more variables and a starting point, find a fixed-point of the function: i.e. where func(x0) == x0.

Parameters

funcfunction

Function to evaluate.

x0array_like

Fixed point of function.

argstuple, optional

Extra arguments to func.

xtolfloat, optional

Convergence tolerance, defaults to 1e-08.

maxiterint, optional

Maximum number of iterations, defaults to 500.

method{“del2”, “iteration”}, optional

Method of finding the fixed-point, defaults to “del2” which uses Steffensen’s Method with Aitken’s Del^2 convergence acceleration [1]. The “iteration” method simply iterates the function until convergence is detected, without attempting to accelerate the convergence.

References

1

Burden, Faires, “Numerical Analysis”, 5th edition, pg. 80

Examples

>>> from scipy import optimize
>>> def func(x, c1, c2):
...    return np.sqrt(c1/(x+c2))
>>> c1 = np.array([10,12.])
>>> c2 = np.array([3, 5.])
>>> optimize.fixed_point(func, [1.2, 1.3], args=(c1,c2))
array([ 1.4920333 ,  1.37228132])

dlnpyutils.minpack.fsolve(func, x0, args=(), fprime=None, full_output=0, col_deriv=0, xtol=1.49012e-08, maxfev=0, band=None, epsfcn=None, factor=100, diag=None)

Find the roots of a function.

Return the roots of the (non-linear) equations defined by func(x) = 0 given a starting estimate.

Parameters

funccallable f(x, *args)

A function that takes at least one (possibly vector) argument, and returns a value of the same length.

x0ndarray

The starting estimate for the roots of func(x) = 0.

argstuple, optional

Any extra arguments to func.

fprimecallable f(x, *args), optional

A function to compute the Jacobian of func with derivatives across the rows. By default, the Jacobian will be estimated.

full_outputbool, optional

If True, return optional outputs.

col_derivbool, optional

Specify whether the Jacobian function computes derivatives down the columns (faster, because there is no transpose operation).

xtolfloat, optional

The calculation will terminate if the relative error between two consecutive iterates is at most xtol.

maxfevint, optional

The maximum number of calls to the function. If zero, then 100*(N+1) is the maximum where N is the number of elements in x0.

bandtuple, optional

If set to a two-sequence containing the number of sub- and super-diagonals within the band of the Jacobi matrix, the Jacobi matrix is considered banded (only for fprime=None).

epsfcnfloat, optional

A suitable step length for the forward-difference approximation of the Jacobian (for fprime=None). If epsfcn is less than the machine precision, it is assumed that the relative errors in the functions are of the order of the machine precision.

factorfloat, optional

A parameter determining the initial step bound (factor * || diag * x||). Should be in the interval (0.1, 100).

diagsequence, optional

N positive entries that serve as a scale factors for the variables.

Returns

xndarray

The solution (or the result of the last iteration for an unsuccessful call).

infodictdict

A dictionary of optional outputs with the keys:

nfev

number of function calls

njev

number of Jacobian calls

fvec

function evaluated at the output

fjac

the orthogonal matrix, q, produced by the QR factorization of the final approximate Jacobian matrix, stored column wise

r

upper triangular matrix produced by QR factorization of the same matrix

qtf

the vector (transpose(q) * fvec)

ierint

An integer flag. Set to 1 if a solution was found, otherwise refer to mesg for more information.

mesgstr

If no solution is found, mesg details the cause of failure.

See also

root

Interface to root finding algorithms for multivariate functions. See the method=='hybr' in particular.

Notes

fsolve is a wrapper around MINPACK’s hybrd and hybrj algorithms.

dlnpyutils.minpack.leastsq(func, x0, args=(), Dfun=None, full_output=0, col_deriv=0, ftol=1.49012e-08, xtol=1.49012e-08, gtol=0.0, maxfev=0, epsfcn=None, factor=100, diag=None)

Minimize the sum of squares of a set of equations.

x = arg min(sum(func(y)**2,axis=0))
         y

Parameters

funccallable

should take at least one (possibly length N vector) argument and returns M floating point numbers. It must not return NaNs or fitting might fail.

x0ndarray

The starting estimate for the minimization.

argstuple, optional

Any extra arguments to func are placed in this tuple.

Dfuncallable, optional

A function or method to compute the Jacobian of func with derivatives across the rows. If this is None, the Jacobian will be estimated.

full_outputbool, optional

non-zero to return all optional outputs.

col_derivbool, optional

non-zero to specify that the Jacobian function computes derivatives down the columns (faster, because there is no transpose operation).

ftolfloat, optional

Relative error desired in the sum of squares.

xtolfloat, optional

Relative error desired in the approximate solution.

gtolfloat, optional

Orthogonality desired between the function vector and the columns of the Jacobian.

maxfevint, optional

The maximum number of calls to the function. If Dfun is provided then the default maxfev is 100*(N+1) where N is the number of elements in x0, otherwise the default maxfev is 200*(N+1).

epsfcnfloat, optional

A variable used in determining a suitable step length for the forward- difference approximation of the Jacobian (for Dfun=None). Normally the actual step length will be sqrt(epsfcn)*x If epsfcn is less than the machine precision, it is assumed that the relative errors are of the order of the machine precision.

factorfloat, optional

A parameter determining the initial step bound (factor * || diag * x||). Should be in interval (0.1, 100).

diagsequence, optional

N positive entries that serve as a scale factors for the variables.

Returns

xndarray

The solution (or the result of the last iteration for an unsuccessful call).

cov_xndarray

The inverse of the Hessian. fjac and ipvt are used to construct an estimate of the Hessian. A value of None indicates a singular matrix, which means the curvature in parameters x is numerically flat. To obtain the covariance matrix of the parameters x, cov_x must be multiplied by the variance of the residuals – see curve_fit.

infodictdict

a dictionary of optional outputs with the keys:

nfev

The number of function calls

fvec

The function evaluated at the output

fjac

A permutation of the R matrix of a QR factorization of the final approximate Jacobian matrix, stored column wise. Together with ipvt, the covariance of the estimate can be approximated.

ipvt

An integer array of length N which defines a permutation matrix, p, such that fjac*p = q*r, where r is upper triangular with diagonal elements of nonincreasing magnitude. Column j of p is column ipvt(j) of the identity matrix.

qtf

The vector (transpose(q) * fvec).

mesgstr

A string message giving information about the cause of failure.

ierint

An integer flag. If it is equal to 1, 2, 3 or 4, the solution was found. Otherwise, the solution was not found. In either case, the optional output variable ‘mesg’ gives more information.

See also

least_squares

Newer interface to solve nonlinear least-squares problems with bounds on the variables. See method=='lm' in particular.

Notes

“leastsq” is a wrapper around MINPACK’s lmdif and lmder algorithms.

cov_x is a Jacobian approximation to the Hessian of the least squares objective function. This approximation assumes that the objective function is based on the difference between some observed target data (ydata) and a (non-linear) function of the parameters f(xdata, params)

func(params) = ydata - f(xdata, params)

so that the objective function is

  min   sum((ydata - f(xdata, params))**2, axis=0)
params

The solution, x, is always a 1D array, regardless of the shape of x0, or whether x0 is a scalar.

dlnpyutils.mpcommon module

Functions used by least-squares algorithms.

dlnpyutils.mpcommon.CL_scaling_vector(x, g, lb, ub)

Compute Coleman-Li scaling vector and its derivatives.

Components of a vector v are defined as follows:

       | ub[i] - x[i], if g[i] < 0 and ub[i] < np.inf
v[i] = | x[i] - lb[i], if g[i] > 0 and lb[i] > -np.inf
       | 1,           otherwise

According to this definition v[i] >= 0 for all i. It differs from the definition in paper [1] (eq. (2.2)), where the absolute value of v is used. Both definitions are equivalent down the line. Derivatives of v with respect to x take value 1, -1 or 0 depending on a case.

Returns

vndarray with shape of x

Scaling vector.

dvndarray with shape of x

Derivatives of v[i] with respect to x[i], diagonal elements of v’s Jacobian.

References

1

M.A. Branch, T.F. Coleman, and Y. Li, “A Subspace, Interior, and Conjugate Gradient Method for Large-Scale Bound-Constrained Minimization Problems,” SIAM Journal on Scientific Computing, Vol. 21, Number 1, pp 1-23, 1999.

dlnpyutils.mpcommon.build_quadratic_1d(J, g, s, diag=None, s0=None)

Parameterize a multivariate quadratic function along a line.

The resulting univariate quadratic function is given as follows:

f(t) = 0.5 * (s0 + s*t).T * (J.T*J + diag) * (s0 + s*t) +
       g.T * (s0 + s*t)

Parameters

Jndarray, sparse matrix or LinearOperator shape (m, n)

Jacobian matrix, affects the quadratic term.

gndarray, shape (n,)

Gradient, defines the linear term.

sndarray, shape (n,)

Direction vector of a line.

diagNone or ndarray with shape (n,), optional

Addition diagonal part, affects the quadratic term. If None, assumed to be 0.

s0None or ndarray with shape (n,), optional

Initial point. If None, assumed to be 0.

Returns

afloat

Coefficient for t**2.

bfloat

Coefficient for t.

cfloat

Free term. Returned only if s0 is provided.

dlnpyutils.mpcommon.check_termination(dF, F, dx_norm, x_norm, ratio, ftol, xtol)

Check termination condition for nonlinear least squares.

dlnpyutils.mpcommon.compute_grad(J, f)

Compute gradient of the least-squares cost function.

dlnpyutils.mpcommon.compute_jac_scale(J, scale_inv_old=None)

Compute variables scale based on the Jacobian matrix.

dlnpyutils.mpcommon.evaluate_quadratic(J, g, s, diag=None)

Compute values of a quadratic function arising in least squares.

The function is 0.5 * s.T * (J.T * J + diag) * s + g.T * s.

Parameters

Jndarray, sparse matrix or LinearOperator, shape (m, n)

Jacobian matrix, affects the quadratic term.

gndarray, shape (n,)

Gradient, defines the linear term.

sndarray, shape (k, n) or (n,)

Array containing steps as rows.

diagndarray, shape (n,), optional

Addition diagonal part, affects the quadratic term. If None, assumed to be 0.

Returns

valuesndarray with shape (k,) or float

Values of the function. If s was 2-dimensional then ndarray is returned, otherwise float is returned.

dlnpyutils.mpcommon.find_active_constraints(x, lb, ub, rtol=1e-10)

Determine which constraints are active in a given point.

The threshold is computed using rtol and the absolute value of the closest bound.

Returns

activendarray of int with shape of x

Each component shows whether the corresponding constraint is active:

• 0 - a constraint is not active.

• -1 - a lower bound is active.

• 1 - a upper bound is active.

dlnpyutils.mpcommon.in_bounds(x, lb, ub)

Check if a point lies within bounds.

dlnpyutils.mpcommon.intersect_trust_region(x, s, Delta)

Find the intersection of a line with the boundary of a trust region.

This function solves the quadratic equation with respect to t ||(x + s*t)||**2 = Delta**2.

Returns

t_neg, t_postuple of float

Negative and positive roots.

Raises

ValueError

If s is zero or x is not within the trust region.

dlnpyutils.mpcommon.left_multiplied_operator(J, d)

Return diag(d) J as LinearOperator.

dlnpyutils.mpcommon.left_multiply(J, d, copy=True)

Compute diag(d) J.

If copy is False, J is modified in place (unless being LinearOperator).

dlnpyutils.mpcommon.make_strictly_feasible(x, lb, ub, rstep=1e-10)

Shift a point to the interior of a feasible region.

Each element of the returned vector is at least at a relative distance rstep from the closest bound. If rstep=0 then np.nextafter is used.

dlnpyutils.mpcommon.minimize_quadratic_1d(a, b, lb, ub, c=0)

Minimize a 1-d quadratic function subject to bounds.

The free term c is 0 by default. Bounds must be finite.

Returns

tfloat

Minimum point.

yfloat

Minimum value.

dlnpyutils.mpcommon.print_header_linear()

dlnpyutils.mpcommon.print_header_nonlinear()

dlnpyutils.mpcommon.print_iteration_linear(iteration, cost, cost_reduction, step_norm, optimality)

dlnpyutils.mpcommon.print_iteration_nonlinear(iteration, nfev, cost, cost_reduction, step_norm, optimality)

dlnpyutils.mpcommon.reflective_transformation(y, lb, ub)

Compute reflective transformation and its gradient.

dlnpyutils.mpcommon.regularized_lsq_operator(J, diag)

Return a matrix arising in regularized least squares as LinearOperator.

The matrix is

[ J ] [ D ]

where D is diagonal matrix with elements from diag.

dlnpyutils.mpcommon.right_multiplied_operator(J, d)

Return J diag(d) as LinearOperator.

dlnpyutils.mpcommon.right_multiply(J, d, copy=True)

Compute J diag(d).

If copy is False, J is modified in place (unless being LinearOperator).

dlnpyutils.mpcommon.scale_for_robust_loss_function(J, f, rho)

Scale Jacobian and residuals for a robust loss function.

Arrays are modified in place.

dlnpyutils.mpcommon.solve_lsq_trust_region(n, m, uf, s, V, Delta, initial_alpha=None, rtol=0.01, max_iter=10)

Solve a trust-region problem arising in least-squares minimization.

This function implements a method described by J. J. More [1] and used in MINPACK, but it relies on a single SVD of Jacobian instead of series of Cholesky decompositions. Before running this function, compute: U, s, VT = svd(J, full_matrices=False).

Parameters

nint

Number of variables.

mint

Number of residuals.

ufndarray

Computed as U.T.dot(f).

sndarray

Singular values of J.

Vndarray

Transpose of VT.

Deltafloat

Radius of a trust region.

initial_alphafloat, optional

Initial guess for alpha, which might be available from a previous iteration. If None, determined automatically.

rtolfloat, optional

Stopping tolerance for the root-finding procedure. Namely, the solution p will satisfy abs(norm(p) - Delta) < rtol * Delta.

max_iterint, optional

Maximum allowed number of iterations for the root-finding procedure.

Returns

pndarray, shape (n,)

Found solution of a trust-region problem.

alphafloat

Positive value such that (J.T*J + alpha*I)*p = -J.T*f. Sometimes called Levenberg-Marquardt parameter.

n_iterint

Number of iterations made by root-finding procedure. Zero means that Gauss-Newton step was selected as the solution.

References

1

More, J. J., “The Levenberg-Marquardt Algorithm: Implementation and Theory,” Numerical Analysis, ed. G. A. Watson, Lecture Notes in Mathematics 630, Springer Verlag, pp. 105-116, 1977.

dlnpyutils.mpcommon.solve_trust_region_2d(B, g, Delta)

Solve a general trust-region problem in 2 dimensions.

The problem is reformulated as a 4-th order algebraic equation, the solution of which is found by numpy.roots.

Parameters

Bndarray, shape (2, 2)

Symmetric matrix, defines a quadratic term of the function.

gndarray, shape (2,)

Defines a linear term of the function.

Deltafloat

Radius of a trust region.

Returns

pndarray, shape (2,)

Found solution.

newton_stepbool

Whether the returned solution is the Newton step which lies within the trust region.

dlnpyutils.mpcommon.step_size_to_bound(x, s, lb, ub)

Compute a min_step size required to reach a bound.

The function computes a positive scalar t, such that x + s * t is on the bound.

Returns

stepfloat

Computed step. Non-negative value.

hitsndarray of int with shape of x

Each element indicates whether a corresponding variable reaches the bound:

• 0 - the bound was not hit.

• -1 - the lower bound was hit.

• 1 - the upper bound was hit.

dlnpyutils.mpcommon.update_tr_radius(Delta, actual_reduction, predicted_reduction, step_norm, bound_hit)

Update the radius of a trust region based on the cost reduction.

Returns

Deltafloat

New radius.

ratiofloat

Ratio between actual and predicted reductions.

dlnpyutils.plotting module

dlnpyutils.robust module

This is from the LWA Software Library (LSL) https://github.com/lwa-project/lsl Redistributed under the GNU license.

Small collection of robust statistical estimators based on functions from Henry Freudenriech (Hughes STX) statistics library (called ROBLIB) that have been incorporated into the AstroIDL User’s Library. Function included are:

• biweight_mean - biweighted mean estimator

• mean - robust estimator of the mean of a data set

• mode - robust estimate of the mode of a data set using the half-sample

  method

• std - robust estimator of the standard deviation of a data set

• checkfit - return the standard deviation and biweights for a fit in order

  to determine its quality

• linefit - outlier resistant fit of a line to data

• polyfit - outlier resistant fit of a polynomial to data

For the fitting routines, the coefficients are returned in the same order as numpy.polyfit, i.e., with the coefficient of the highest power listed first.

For additional information about the original IDL routines, see: http://idlastro.gsfc.nasa.gov/contents.html#C17

dlnpyutils.robust.biweight_mean(inputData, axis=None, dtype=None)

Calculate the mean of a data set using bisquare weighting.

Based on the biweight_mean routine from the AstroIDL User’s Library.

Changed in version 1.0.3: Added the ‘axis’ and ‘dtype’ keywords to make this function more compatible with numpy.mean()

dlnpyutils.robust.checkfit(inputData, inputFit, epsilon, delta, bisquare_limit=6.0)

Determine the quality of a fit and biweights. Returns a tuple with elements:

0. Status

1. Robust standard deviation analog

2. Fractional median absolute deviation of the residuals

3. Number of input points given non-zero weight in the calculation

4. Bisquare weights of the input points

5. Residual values scaled by sigma

This function is based on the rob_checkfit routine from the AstroIDL User’s Library.

dlnpyutils.robust.linefit(inputX, inputY, max_iter=25, bisector=False, bisquare_limit=6.0, close_factor=0.03)

Outlier resistance two-variable linear regression function.

Based on the robust_linefit routine in the AstroIDL User’s Library.

dlnpyutils.robust.mean(inputData, cut=3.0, axis=None, dtype=None)

Robust estimator of the mean of a data set. Based on the resistant_mean function from the AstroIDL User’s Library.

Changed in version 1.2.1: Added a ValueError if the distriubtion is too strange

Changed in version 1.0.3: Added the ‘axis’ and ‘dtype’ keywords to make this function more compatible with numpy.mean()

dlnpyutils.robust.mode(inputData, axis=None, dtype=None)

Robust estimator of the mode of a data set using the half-sample mode.

dlnpyutils.robust.polyfit(inputX, inputY, order, max_iter=25)

Outlier resistance two-variable polynomial function fitter.

Based on the robust_poly_fit routine in the AstroIDL User’s Library.

dlnpyutils.robust.std(inputData, zero=False, axis=None, dtype=None)

Robust estimator of the standard deviation of a data set.

Based on the robust_sigma function from the AstroIDL User’s Library.

Changed in version 1.2.1: Added a ValueError if the distriubtion is too strange

Changed in version 1.0.3: Added the ‘axis’ and ‘dtype’ keywords to make this function more compatible with numpy.std()

dlnpyutils.spec module

dlnpyutils.trf module

Trust Region Reflective algorithm for least-squares optimization.

The algorithm is based on ideas from paper [STIR]. The main idea is to account for presence of the bounds by appropriate scaling of the variables (or equivalently changing a trust-region shape). Let’s introduce a vector v:

ub[i] - x[i], if g[i] < 0 and ub[i] < np.inf

v[i] = | x[i] - lb[i], if g[i] > 0 and lb[i] > -np.inf

1, otherwise

where g is the gradient of a cost function and lb, ub are the bounds. Its components are distances to the bounds at which the anti-gradient points (if this distance is finite). Define a scaling matrix D = diag(v**0.5). First-order optimality conditions can be stated as

D^2 g(x) = 0.

Meaning that components of the gradient should be zero for strictly interior variables, and components must point inside the feasible region for variables on the bound.

Now consider this system of equations as a new optimization problem. If the point x is strictly interior (not on the bound) then the left-hand side is differentiable and the Newton step for it satisfies

(D^2 H + diag(g) Jv) p = -D^2 g

where H is the Hessian matrix (or its J^T J approximation in least squares), Jv is the Jacobian matrix of v with components -1, 1 or 0, such that all elements of matrix C = diag(g) Jv are non-negative. Introduce the change of the variables x = D x_h (_h would be “hat” in LaTeX). In the new variables we have a Newton step satisfying

B_h p_h = -g_h,

where B_h = D H D + C, g_h = D g. In least squares B_h = J_h^T J_h, where J_h = J D. Note that J_h and g_h are proper Jacobian and gradient with respect to “hat” variables. To guarantee global convergence we formulate a trust-region problem based on the Newton step in the new variables:

0.5 * p_h^T B_h p + g_h^T p_h -> min, ||p_h|| <= Delta

In the original space B = H + D^{-1} C D^{-1}, and the equivalent trust-region problem is

0.5 * p^T B p + g^T p -> min, ||D^{-1} p|| <= Delta

Here the meaning of the matrix D becomes more clear: it alters the shape of a trust-region, such that large steps towards the bounds are not allowed. In the implementation the trust-region problem is solved in “hat” space, but handling of the bounds is done in the original space (see below and read the code).

The introduction of the matrix D doesn’t allow to ignore bounds, the algorithm must keep iterates strictly feasible (to satisfy aforementioned differentiability), the parameter theta controls step back from the boundary (see the code for details).

The algorithm does another important trick. If the trust-region solution doesn’t fit into the bounds, then a reflected (from a firstly encountered bound) search direction is considered. For motivation and analysis refer to [STIR] paper (and other papers of the authors). In practice it doesn’t need a lot of justifications, the algorithm simply chooses the best step among three: a constrained trust-region step, a reflected step and a constrained Cauchy step (a minimizer along -g_h in “hat” space, or -D^2 g in the original space).

Another feature is that a trust-region radius control strategy is modified to account for appearance of the diagonal C matrix (called diag_h in the code).

Note, that all described peculiarities are completely gone as we consider problems without bounds (the algorithm becomes a standard trust-region type algorithm very similar to ones implemented in MINPACK).

The implementation supports two methods of solving the trust-region problem. The first, called ‘exact’, applies SVD on Jacobian and then solves the problem very accurately using the algorithm described in [JJMore]. It is not applicable to large problem. The second, called ‘lsmr’, uses the 2-D subspace approach (sometimes called “indefinite dogleg”), where the problem is solved in a subspace spanned by the gradient and the approximate Gauss-Newton step found by scipy.sparse.linalg.lsmr. A 2-D trust-region problem is reformulated as a 4-th order algebraic equation and solved very accurately by numpy.roots. The subspace approach allows to solve very large problems (up to couple of millions of residuals on a regular PC), provided the Jacobian matrix is sufficiently sparse.

References

STIR(1,2)

Branch, M.A., T.F. Coleman, and Y. Li, “A Subspace, Interior, and Conjugate Gradient Method for Large-Scale Bound-Constrained Minimization Problems,” SIAM Journal on Scientific Computing, Vol. 21, Number 1, pp 1-23, 1999.

JJMore

More, J. J., “The Levenberg-Marquardt Algorithm: Implementation and Theory,” Numerical Analysis, ed. G. A. Watson, Lecture

dlnpyutils.trf.select_step(x, J_h, diag_h, g_h, p, p_h, d, Delta, lb, ub, theta)

Select the best step according to Trust Region Reflective algorithm.

dlnpyutils.trf.trf(fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale, loss_function, tr_solver, tr_options, verbose, dx_lim=None)

dlnpyutils.trf.trf_bounds(fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale, loss_function, tr_solver, tr_options, verbose, dx_lim=None)

dlnpyutils.trf.trf_no_bounds(fun, jac, x0, f0, J0, ftol, xtol, gtol, max_nfev, x_scale, loss_function, tr_solver, tr_options, verbose)

dlnpyutils.utils module

Module contents