spectre.least_squares.common

Functions used by least-squares algorithms.

Module Contents

spectre.least_squares.common.EPS

spectre.least_squares.common.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.

spectre.least_squares.common.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]

J. J. Moré. The levenberg-marquardt algorithm: implementation and theory. In G. A. Watson, editor, Numerical Analysis, volume 630 of Lecture Notes in Mathematics, pages 105–116. Springer, 1977.

spectre.least_squares.common.solve_trust_region_2d(B, g, Delta)

Solve a general trust-region problem in 2 dimensions.

The problem is reformulated as a 4th 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.

spectre.least_squares.common.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.

spectre.least_squares.common.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.

spectre.least_squares.common.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.

spectre.least_squares.common.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-D, then ndarray is returned, otherwise, float is returned.

spectre.least_squares.common.in_bounds(x, lb, ub)

Check if a point lies within bounds.

spectre.least_squares.common.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.

spectre.least_squares.common.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.

spectre.least_squares.common.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.

spectre.least_squares.common.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 Branch et al.[2] (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

[2]

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, 21(1):1–23, 1999.

spectre.least_squares.common.reflective_transformation(y, lb, ub)

Compute reflective transformation and its gradient.

spectre.least_squares.common.compute_grad(J, f)

Compute gradient of the least-squares cost function.

spectre.least_squares.common.compute_jac_scale(J, scale_inv_old=None)

Compute variables scale based on the Jacobian matrix.

spectre.least_squares.common.left_multiplied_operator(J, d)

Return diag(d) J as LinearOperator.

spectre.least_squares.common.right_multiplied_operator(J, d)

Return J diag(d) as LinearOperator.

spectre.least_squares.common.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.

spectre.least_squares.common.right_multiply(J, d, copy=True)

Compute J diag(d).

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

spectre.least_squares.common.left_multiply(J, d, copy=True)

Compute diag(d) J.

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

spectre.least_squares.common.check_termination(dF, F, dx_norm, x_norm, ratio, ftol, xtol)

Check termination condition for nonlinear least squares.

spectre.least_squares.common.scale_for_robust_loss_function(J, f, rho)

Scale Jacobian and residuals for a robust loss function.

Arrays are modified in place.

class spectre.least_squares.common.MemoizeJac(fun)

Decorator that caches the return values of a function returning (fun, grad) each time it is called.

fun

jac = None

x = None

derivative(x, *args)