capfit 2.6.5

CapFit: linearly-constrained non-linear least-squares optimization

The CapFit Package

CapFit: Linearly-Constrained Non-Linear Least-Squares Optimization

The CapFit package provides a Python implementation of a linearly-constrained non-linear least-squares optimization method. This method is described in Section 3.2 of the paper by Cappellari (2023).

Contents

• The CapFit Package

  • Attribution

  • Installation

  • Usage Examples

• CapFit

  • CapFit Purpose

  • Calling Sequence

  • Parameters

  • Other Parameters

  • Returns

• lsq_box

  • Parameters

  • Returns

• lsq_lin

  • Parameters

  • Returns

• License

  • Changelog

Attribution

If you use this software for your research, please cite the Cappellari (2023) paper where the algorithm was introduced. The BibTeX entry for the paper is:

@ARTICLE{Cappellari2023,
    author = {{Cappellari}, M.},
    title = "{Full spectrum fitting with photometry in PPXF: stellar population
        versus dynamical masses, non-parametric star formation history and
        metallicity for 3200 LEGA-C galaxies at redshift $z\approx0.8$}",
    journal = {MNRAS},
    eprint = {2208.14974},
    year = 2023,
    volume = 526,
    pages = {3273-3300},
    doi = {10.1093/mnras/stad2597}
}

Installation

To install CapFit, use the following command:

pip install capfit

If you don’t have write access to the global site-packages directory, install it with:

pip install --user capfit

To upgrade to the latest version, run:

pip install --upgrade capfit

Usage Examples

Explore how to use the CapFit package by referring to the capfit_examples.py file located within the capfit/examples directory. You can find this file in the main CapFit package installation folder inside the site-packages directory. Detailed documentation is available in the docstring of the capfit.py file or on PyPi.

---

CapFit

CapFit Purpose

CapFit solves linearly-constrained non-linear least-squares optimization problems.

Linear inequality/equality constraints and bound constraints are supported. Moreover, one can easily tie (i.e. enforce non-linear constraints) or fix some parameters without having to modify the fitting function.

CapFit implements Algorithm 2 of Cappellari (2023).

CapFit combines two successful ideas:

1. The Sequential Quadratic Programming (SQP) approach, specialized for the case of linear constraints;

2. The Levenberg-Marquardt (LM) method.

It was designed for the common situations where the user function is not a simple analytic function but is the result of some complex calculations and is more expensive to compute than the small quadratic subproblem.

CapFit performance is generally better than the best unconstrained or bound-constrained least-squares algorithms currently available, in terms of robustness and number of function evaluations. However, in addition CapFit allows for more general constraints.

Given a function of n model parameters x_k returning the m model residuals f_j(x), CapFit finds a local minimum of the cost function:

G(x) = sum[f_j(x)^2]

subject to:

A_ineq @ x <= b_ineq            # Linear Inequality Constraints
A_eq @ x == b_eq                # Linear Equality Constraints
bounds[0] <= x <= bounds[1]     # Bounds
x_k = f(x)                      # Tied Parameters
x_k = a_k                       # Fixed Parameters

Calling Sequence

from capfit.capfit import capfit

res = capfit(func, x0, A_eq=None, A_ineq=None, abs_step=None,
        b_eq=None, b_ineq=None, bounds=(-np.inf, np.inf),
        diff_step=1e-4, fixed=None, ftol=1e-4, linear_method='lsq_lin',
        max_nfev=None, monitor=None, rcond=None, tied=None, verbose=0,
        x_scale='jac', xtol=1e-4, args=(), kwargs={})

print(f"solution: {res.x}")

Parameters

fun: callable

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 1-d darray of shape (n,). The function must return a 1-d array of shape (m,).

x0: array_like with shape (n,) or float

Initial guess on independent variables. For guaranteed convergence, the initial guess must be feasible (satsfies the constraints) and for this reason an error is returned if it is not the case.

Other Parameters

abs_step: None, scalar or array_like, optional

Determines the absolute step size for the finite difference approximation of the Jacobian. If abs_step is given, then the value of diff_step is ignored and the finite differences use this fixed absolute step.

A_eq: array_like with shape (q, n), optional

Defines the linear equality constraints on the fitted parameters:

A_eq @ x == b_eq

The same result can be achieved using the tied keyword, which also allows for non-linear equality constraints.

A_ineq: array_like with shape (p, n), optional

Defines the linear inequality constraints on the fitted parameters:

A_ineq @ x <= b_ineq

b_ineq: array_like with shape (p,), optional

See description of A_ineq.

b_eq: array_like with shape (q,), optional

See description of A_eq.

bounds: 2-tuple of array_like, optional

Lower and upper bounds (lb, ub) 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.

diff_step: None, scalar or array_like, optional

Determines the relative step size for the finite difference approximation of the Jacobian. The actual step is computed as diff_step*maximum(1, abs(x)) (default diff_step=1e-4)

fixed: boolean array_like with shape (n,), optional

Boolean vector set to True where a given parameter has to be held fixed with the value given in x0.

ftol: float or None, optional

Tolerance for termination by the change of the cost function (default is 1e-4). The optimization process is stopped when both prered < ftol and abs(actred) < ftol and additionally actred <= 2*prered, where actred and prered are the actual and predicted relative changes respectively (as described in More’ et al. 1980).

linear_method: {‘lsq_lin’, ‘cvxopt’} str, optional

Method used for the solution of the linear least-squares sub-problem. When using linear constraints, this can be either ‘lsq_lin’ (default) or ‘cvxopt’. In the latter case, the cvxopt package must be installed. With only bounds, this keyword is ignored as capfit uses ‘lsq_box’.

max_nfev: None or int, optional

Maximum number of function evaluations before the termination (default is 100*n).

monitor: dict, optional

Dictionary {"fun": monitor_fun, "num": monitor_num} containing the name of a function with signature monitor_fun(x, niter, chi2) to be called every monitor_num succesful iterations niter of the CapFit algorithm, to monitor the progression of the optimization. The function is also called on the first/last iterations. In most cases, the verbose keyword is an easier alternative, but there are situations where this keyword is needed.

rcond: float, optional

Cutoff for small singular values used to determine the effective rank of the Jacobian. Singular values smaller than rcond*largest_singular_value are considered zero.

tied: array_like with shape (n,), optional

A list of string expressions. Each expression “ties” the parameter to other free or fixed parameters. Any expression involving constants and the parameter array p[j] are permitted. Since they are totally constrained, tied parameters are considered to be fixed; no errors are computed for them.

This is a vector with the same dimensions as x0. In practice, for every element of x0 one needs to specify either an empty string '' implying that the parameter is free, or a string expression involving some of the variables p[j], where j represents the index of the vector of parameters. See usage in capfit/examples.

verbose: {0, 1, 2}, optional

Level of algorithm’s verbosity:

• 0: Works silently (default).

• 1: display a termination report.

• 2: display progress during iterations.

x_scale: array_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 initial size/2 of the box trust region along j-th dimension is given by x_scale[j] and the box ratains its shape during the optimization. Improved convergence is 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 More’ 1978).

xtol: float or None, optional

Tolerance for termination by the change dx of the independent variables (default is 1e-4). The condition is norm(dx) < xtol*(xtol + norm(xs)) where xs is the value of x scaled according to x_scale parameter (see below). If None, the termination by this condition is disabled.

args, kwargs: tuple and dict, optional

Additional arguments passed to fun, empty by default. The calling signature is fun(x, *args, **kwargs).

Returns

The following attributes of the capfit class are defined in output:

.x: ndarray, shape (n,)

Solution found.

.x_err: ndarray, shape (n,)

Formal uncertainties of the solution estimated from diagonal elements of the covariance matrix (i.e. ignoring covariance). This is only meaningful if the fun values represent residuals divided by their 1sigma errors. Uncertainties of the non-free parameters are returned as zero.

.chi2: float

Value of the chi^2 = fun @ fun at the solution.

.cost: float

Value of the cost = 0.5*(fun @ fun) function at the solution.

.cov: ndarray, shape (nfree, nfree)

Covariance matrix for the free (not tied or fixed) parameters.

.fun: ndarray, shape (m,)

Vector of residuals at the solution.

.jac: ndarray, shape (m, nfree)

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. NB: This is the Jacobian of the free (not tied or fixed) parameters.

.grad: ndarray, shape (m,)

Gradient of the cost function at the solution.

.nfev: int

Number of function evaluations done.

.njev: int

Number of Jacobian evaluations done.

.status: int

The reason for algorithm termination:

• -1: improper input parameters status

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

• 2: ftol termination condition is satisfied.

• 3: xtol termination condition is satisfied.

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

.message: str

Verbal description of the termination reason.

.success: bool

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

---

lsq_box

Fast linear least-squares with box (bounds) constraints.

Given a matrix A and vector b, and an optional starting guess for x, find a vector x which solves:

Minimize      || A @ x - b ||
Subject to    bounds[0] <= x <= bounds[1]

This implementation is described in Section 3.2.3 of the paper by Cappellari (2023).

It is an improved and faster version of the active-set BVLS algorithm in Appendix C of Lawson & Hanson (1995).

It allows for an initial value for x, which can dramatically speed up convergence when a good guess is available. It also includes an initialization phase which significantly reduces the number of iterations in the main loop. This is Algorithm 3 of Cappellari (2023)

Parameters

A: array_like with size (m, n)

Design matrix of the linear least-squares problem.

b: array_like with size (m,)

Target vector of the linear least-squares problem.

bounds: array_like with size (2, n)

Lower and upper bounds on parameters. Defaults to no bounds.

x: array_like with size (n,), optional

Starting guess for the solution vector. A good guess can dramatically speed up convergence.

rcond: float, optional

Cutoff for small singular values used to determine the effective rank of A. Singular values smaller than rcond*largest_singular_value are considered zero.

ftol: float, optional

Safety criterion to stop iterating. The minimum is generally computed to machine accuracy regardless of this value.

verbose: bool, optional

Set verbose=True to print the intermediate results.

Returns

Set as attributes of the lsq_box class.

.x: array_like with size (n,), optional

Solution vector.

.active_mask: array_like of int with shape (n,)

Each component shows whether a corresponding constraint is active.

.cost: float

Value of the cost function chi^2/2 at the solution.

---

lsq_lin

Linear least-squares with linear equality and/or inequality constraints.

Given a matrix A and vector b, find x which solves:

Minimize      || A @ x - b ||
Subject to    A_ineq @ x <= b_ineq
and           A_eq @ x == b_eq
and           bounds[0] <= x <= bounds[1]

This implementation is described in Section 3.2.3 of the paper by Cappellari (2023).

The main loop of this algorithm is a generalization of the active-set NNLS Algorithm (23.10) in Lawson & Hanson (1995). The generalization to linear constraints follows the ideas presented in Sec. 16.5 of Nocedal & Wright (2006).

Parameters

A: array_like with size (m, n)

Design matrix of the linear least-squares problem.

b: array_like with size (m,)

Target vector of the linear least-squares problem.

A_ineq: array_like with size (q, n)

Matrix defining the q linear inequality constraints.

b_ineq: array_like with size (q,)

Vector defining the q linear inequality constraints.

A_eq: array_like with size (p, n)

Matrix defining the q linear equality constraints.

b_eq: array_like with size (p,)

Vector defining the p linear equality constraints.

bounds: array_like with size (2, n)

Lower and upper bounds on parameters. Defaults to no bounds.

x: array_like with size (n,), optional

Starting guess for the solution vector. A good guess can dramatically speed up convergence.

rcond: float, optional

Cutoff for small singular values used to determine the effective rank of A. Singular values smaller than rcond*largest_singular_value are considered zero.

ftol: float, optional

Safety criterion to stop iterating. The minimum is generally computed to machine accuracy regardless of this value.

verbose: bool, optional

Set verbose=True to print the intermediate results.

Returns

Set as attributes of the lsq_lin class.

.x: array_like with size (n,), optional

Solution vector.

.active_mask: array_like of int with shape (n,)

Each component shows whether a corresponding constraint is active.

.cost: float

Value of the cost function chi^2/2 at the solution.

---

License

Other/Proprietary License

Copyright (c) 2017-2024 Michele Cappellari

This software is provided as is with no warranty. You may use it for non-commercial purposes and modify it for personal or internal use, as long as you include this copyright and disclaimer in all copies. You may not redistribute the code.

---

Changelog

V2.6.5: MC, Oxford, 05 August 2024

• Removed capfit from the ppxf package and made it a separated package.

• capfit: Introduced a new keyword, monitor to define a user function for optimization monitoring.

V2.5.1: MC, Oxford, 22 May 2023

• capfit: Relaxed tolerance when checking initial guess feasibility.

V2.5.0: MC, Oxford, 16 August 2022

• Uses scipy.optimize.linprog to find feasible starting point in lsq_lin.

• Set default linear_method='lsq_lin' in capfit. This eliminates the need to install cvxopt when using general linear constraints.

V2.4.0: MC, Oxford, 04 March 2022

• Remove the non-free variables before the optimization. This reduces the degeneracy of the Jacobian.

V2.3.0: MC, Oxford, 20 December 2020

• New linear_method keyword to select cvxopt or lsq_lin, for cases where the latter stops, when using linear constraints. Thanks to Kyle Westfall (UCO Lick) for a detailed bug report.

V2.2.1: MC, Oxford, 11 September 2020

• Fixed possible infinite loop in lsq_box and lsq_lin. Thanks to Shravan Shetty (pku.edu.cn) for the detailed report.

• Use NumPy rather than SciPy version of linalg.lstsq to avoid a Scipy bug in the default criterion for rank deficiency.

• Pass rcond keyword to cov_err for consistency.

V2.2.0: MC, Oxford, 20 August 2020

• New function lsq_lin implementing a robust linear least-squares linearly-constrained algorithm which works with a rank-deficient matrix and allows for a starting guess. lsq_lin supersedes the former lsqlin.

• Renamed lsqbox to lsq_box and revised its interface.

V2.1.0: MC, Oxford, 09 July 2020

• New function lsqbox implementing a fast linear least-squares box-constrained (bounds) algorithm which allows for a starting guess.

V2.0.2: MC, Oxford, 20 June 2020

• capfit: new keyword cond (passed to lsqlin).

• capfit: Use bvls to solve quadratic subproblem with only bounds.

V2.0.1: MC, Oxford, 24 January 2020

• New keyword cond for lsqlin.

• Relaxed assertion for inconsistent inequalities in lsqlin to avoid false positives. Thanks to Kyle Westfall (UCO Lick) for a detailed bug report.

V2.0.0: MC, Oxford, 10 January 2020

• Use the new general linear least-squares optimization function lsqlin to solve the quadratic sub-problem.

• Allow for linear inequality/equality constraints A_ineq, b_ineq and A_eq, b_eq

V1.0.7: MC, Oxford, 10 October 2019

• Included complete documentation.

• Improved print formatting.

• Return .message attribute.

• Improved xtol convergence test.

• Only accept the final move if chi2 decreases.

• Strictly satisfy bounds during Jacobian computation.

V1.0.6: MC, Oxford, 11 June 2019

• Use only free parameters for small-step convergence tests.

• Explain in words convergence status when verbose != 0

• Fixed program stops when abs_step is undefined.

• Fixed capfit ignoring optional max_nfev.

V1.0.5: MC, Oxford, 28 March 2019

• Raise an error if the user function returns non-finite values.

V1.0.4: MC, Oxford, 30 November 2018

• Allow for a scalar abs_step.

V1.0.3: MC, Oxford, 20 September 2018

• Raise an error if one tries to tie parameters to themselves. Thanks to Kyle Westfall (Univ. Santa Cruz) for the feedback.

• Use Python 3.6 f-strings.

V1.0.2: MC, Oxford, 10 May 2018

• Dropped legacy Python 2.7 support.

V1.0.1: MC, Oxford, 13 February 2018

• Make output errors of non-free variables exactly zero.

V1.0.0: MC, Oxford, 15 June 2017

• Written by Michele Cappellari