spectre.least_squares.common ============================ .. py:module:: spectre.least_squares.common .. autoapi-nested-parse:: Functions used by least-squares algorithms. Module Contents --------------- .. py:data:: EPS .. py:function:: 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_pos : tuple of float Negative and positive roots. Raises ------ ValueError If `s` is zero or `x` is not within the trust region. .. py:function:: 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\ :footcite:p:`More1977` 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 ---------- n : int Number of variables. m : int Number of residuals. uf : ndarray Computed as U.T.dot(f). s : ndarray Singular values of J. V : ndarray Transpose of VT. Delta : float Radius of a trust region. initial_alpha : float, optional Initial guess for alpha, which might be available from a previous iteration. If None, determined automatically. rtol : float, optional Stopping tolerance for the root-finding procedure. Namely, the solution ``p`` will satisfy ``abs(norm(p) - Delta) < rtol * Delta``. max_iter : int, optional Maximum allowed number of iterations for the root-finding procedure. Returns ------- p : ndarray, shape (n,) Found solution of a trust-region problem. alpha : float Positive value such that (J.T*J + alpha*I)*p = -J.T*f. Sometimes called Levenberg-Marquardt parameter. n_iter : int Number of iterations made by root-finding procedure. Zero means that Gauss-Newton step was selected as the solution. References ---------- .. footbibliography:: .. py:function:: 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 ---------- B : ndarray, shape (2, 2) Symmetric matrix, defines a quadratic term of the function. g : ndarray, shape (2,) Defines a linear term of the function. Delta : float Radius of a trust region. Returns ------- p : ndarray, shape (2,) Found solution. newton_step : bool Whether the returned solution is the Newton step which lies within the trust region. .. py:function:: 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 ------- Delta : float New radius. ratio : float Ratio between actual and predicted reductions. .. py:function:: 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 ---------- J : ndarray, sparse matrix or LinearOperator shape (m, n) Jacobian matrix, affects the quadratic term. g : ndarray, shape (n,) Gradient, defines the linear term. s : ndarray, shape (n,) Direction vector of a line. diag : None or ndarray with shape (n,), optional Addition diagonal part, affects the quadratic term. If None, assumed to be 0. s0 : None or ndarray with shape (n,), optional Initial point. If None, assumed to be 0. Returns ------- a : float Coefficient for t**2. b : float Coefficient for t. c : float Free term. Returned only if `s0` is provided. .. py:function:: 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 ------- t : float Minimum point. y : float Minimum value. .. py:function:: 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 ---------- J : ndarray, sparse matrix or LinearOperator, shape (m, n) Jacobian matrix, affects the quadratic term. g : ndarray, shape (n,) Gradient, defines the linear term. s : ndarray, shape (k, n) or (n,) Array containing steps as rows. diag : ndarray, shape (n,), optional Addition diagonal part, affects the quadratic term. If None, assumed to be 0. Returns ------- values : ndarray with shape (k,) or float Values of the function. If `s` was 2-D, then ndarray is returned, otherwise, float is returned. .. py:function:: in_bounds(x, lb, ub) Check if a point lies within bounds. .. py:function:: 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 ------- step : float Computed step. Non-negative value. hits : ndarray 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. .. py:function:: 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 ------- active : ndarray 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. .. py:function:: 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. .. py:function:: 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 :footcite:t:`Branch1999` (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 ------- v : ndarray with shape of x Scaling vector. dv : ndarray with shape of x Derivatives of v[i] with respect to x[i], diagonal elements of v's Jacobian. References ---------- .. footbibliography:: .. py:function:: reflective_transformation(y, lb, ub) Compute reflective transformation and its gradient. .. py:function:: compute_grad(J, f) Compute gradient of the least-squares cost function. .. py:function:: compute_jac_scale(J, scale_inv_old=None) Compute variables scale based on the Jacobian matrix. .. py:function:: left_multiplied_operator(J, d) Return diag(d) J as LinearOperator. .. py:function:: right_multiplied_operator(J, d) Return J diag(d) as LinearOperator. .. py:function:: 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`. .. py:function:: right_multiply(J, d, copy=True) Compute J diag(d). If `copy` is False, `J` is modified in place (unless being LinearOperator). .. py:function:: left_multiply(J, d, copy=True) Compute diag(d) J. If `copy` is False, `J` is modified in place (unless being LinearOperator). .. py:function:: check_termination(dF, F, dx_norm, x_norm, ratio, ftol, xtol) Check termination condition for nonlinear least squares. .. py:function:: scale_for_robust_loss_function(J, f, rho) Scale Jacobian and residuals for a robust loss function. Arrays are modified in place. .. py:class:: MemoizeJac(fun) Decorator that caches the return values of a function returning `(fun, grad)` each time it is called. .. py:attribute:: fun .. py:attribute:: jac :value: None .. py:attribute:: x :value: None .. py:method:: derivative(x, *args)