Linear system solvers#

Solvers for computation of linear systems of the form \(\mathbf{A}\mathbf{x}=\mathbf{b}\).

These solvers are used under-the-hood by Poly when set_model() is called, with the solver specified by the method argument. The solvers can also be used independently (see examples) for debugging and experimentation purposes.

class equadratures.solver.Solver(method, solver_args={})[source]#

The parent solver class. All other solver subclasses inherit attributes and methods from this parent class.

Note

This parent class should not be used directly, but can be used to select the desired solver subclass (see examples).

Parameters

  • method (string) – The method used for solving the linear system. Options include: compressed-sensing, least-squares, minimum-norm, numerical-integration, least-squares-with-gradients, least-absolute-residual, huber, elastic-net, and relevance-vector-machine.

  • solver_args (dict, optional) –

    Dictionary containing optional arguments centered around the specific solver. Arguments are specific to the requested solver, except for the following generic arguments:

    • verbose (bool): Default value of this input is set to False; when True, the solver prints information to screen.

    • opt (str): Specifies the underlying optimisation method to use. scipy for scipy, and osqp to use the OSQP optimiser (requires cvxpy to be installed).

Examples

Accessing the underlying solver used by a Poly instance
>>> # Fit a 2D polynomial using the elastic-net solver (with some specific solver_args passed)
>>> x = np.linspace(-1,1,80)
>>> xx = np.vstack([x,x]).T
>>> y = xx[:,0]**3 + xx[0,1]**2 + np.random.normal(0,0.1,80)
>>> myparam = eq.Parameter(distribution='uniform',order=3)
>>> mybasis = eq.Basis('total-order')
>>> mypoly = eq.Poly([myparam,myparam],mybasis,method='elastic-net',
>>>                 sampling_args={'mesh': 'user-defined','sample-points':xx,'sample-outputs':y},
>>>                 solver_args={'verbose':True,'crit':'AIC'})
>>> mypoly.set_model()
>>> # Plot the regularisation path using the solver's method
>>> mypoly.solver.plot_regpath()
Using solver subclasses independently
>>> # Example data
>>> A = np.array([[1,2,3],[1,1,1],[2,1,1]])
>>> b = np.array([1,0,1])
>>>
>>> # # select_solver will return the requested solver subclass
>>> mysolver = eq.Solver.select_solver('least-squares',solver_args={'verbose':False})
>>>
>>> # Or you can select the solversubclass directly
>>> mysolver = solver.least_squares(solver_args={'verbose':True})
>>>
>>> # You can manually run solve and then get_coefficients
>>> mysolver.solve(A,b)
>>> x = mysolver.get_coefficients()
>>>
>>> # Or providing get_coefficients() with A and b will automatically run solve()
>>> x = mysolver.get_coefficients(A,b)
The condition number of the matrix is 29.82452477932781.
The condition number of the matrix is 29.82452477932781.
get_coefficients(*args)[source]#

Get the coefficients the solver has calculated. If the A and b matrices are provied Solver.solve() will be run first. If the solver uses gradients, C and d can also be given.

Parameters

  • A (numpy.ndarray, optional) – Array with shape (number_of_observations, cardinality).

  • b (numpy.ndarray, optional) – Array with shape (number_of_observations, 1) or (number_of_observations).

  • C (numpy.ndarray, optional) – Array with shape (number_of_observations, cardinality).

  • d (numpy.ndarray, optional) – Array with shape (number_of_observations, 1) or (number_of_observations).

Returns

Array containing the coefficients x computed by the solver.

Return type

numpy.ndarray

static select_solver(method, solver_args={})[source]#

Method to return a solver subclass.

Parameters

  • method (str) – The name of the solver to return.

  • solver_args (dict, optional) –

    Dictionary containing optional arguments centered around the specific solver. Arguments are specific to the requested solver, except for the following generic arguments:

    • verbose (bool): Default value of this input is set to False; when True, the solver prints information to screen.

    • opt (str): Specifies the underlying optimisation method to use. scipy for scipy, and osqp to use the OSQP optimiser (requires cvxpy to be installed).

Returns

The Solver subclass specified in method.

Return type

Solver

class equadratures.solver.compressed_sensing(solver_args={})[source]#

Compressed sensing solver.

Parameters

solver_args (dict) –

Optional arguments for the solver:

  • verbose (bool): Default value of this input is set to False; when True, the solver prints information to screen.

  • noise-level (float)

class equadratures.solver.constrained_least_squares(solver_args={})[source]#

Constrained least squares regression.

This solves an orthogonal linear system i.e. simply c=A^T.b.

Parameters

solver_args (dict) –

Optional arguments for the solver:

  • verbose (bool): Default value of this input is set to False; when True, the solver prints information to screen.

class equadratures.solver.custom_solver(solver_args={})[source]#

Custom solver class.

This class allows you to enter a custom solve() function to solve \(\mathbf{A}\mathbf{x}=\mathbf{b}\). The solve() is provided via solver_args, and should accept the numpy.ndarray’s \(\mathbf{A}\) and \(\mathbf{b}\), and return a numpy.ndarray containing the coefficients \(\mathbf{x}\).

Parameters

solver_args (dict) –

Optional arguments for the solver. Optional arguments for the solver. Additional arguments are passed through to the custom solve() function as kwargs.

  • verbose (bool): Default value of this input is set to False; when True, the solver prints information to screen.

  • solve (Callable): The custom solve() solve function. This must accept A and b, and return x.

Example

>>> # Create a new solve method - the Kaczmarz iterative solver
>>> def kaczmarz(A, b, verbose=False):
>>>     m, n = A.shape
>>>     x = np.random.rand(n, 1) # initial guess
>>>     MAXITER = 50
>>>     for iters in range(0, MAXITER):
>>>         if verbose: print('Iteration %d/%d' %(iters,MAXITER))
>>>         for i in range(0, m):
>>>             a_row = A[i, :].reshape(1, n)
>>>             term_3 = float(1.0/(np.dot(a_row , a_row.T)))
>>>             term_1 = float(b[i] - float(np.dot(a_row , x)))
>>>             x = x + (term_1 * term_3 * a_row.T)
>>>     return x
>>>
>>> # Wrap this function in the custom_solver subclass and use it to fit a polynomial
>>> mypoly = eq.Poly(myparam, mybasis, method='custom-solver',
>>>             sampling_args={'mesh':'user-defined', 'sample-points':X, 'sample-outputs':y},
>>>               solver_args={'solve':lsq,'verbose':True})
class equadratures.solver.elastic_net(solver_args={})[source]#

Elastic net solver.

The elastic net solver minimises

\(\frac{1}{2}||\mathbf{A}\mathbf{x}-\mathbf{b}||_2 + \lambda (\alpha||\mathbf{x}||_1 + \frac{1}{2}(1-\alpha)||\mathbf{x}||_2^2)\).

where \(\lambda\) controls the strength of the regularisation (with OLS returned when \(\lambda=0\)), and \(\alpha\) controlling the blending between the L1 and L2 penalty terms.

If path=True, the elastic net is solved via coordinate descent [1]. The full regularisation path (\(\lambda\) path) is computed (for a given \(\alpha\)), and the set of coefficients with the lowest model selection criteria is then selected according to [2]. Otherwise, the elastic net is solved for a given \(\alpha\) and \(\lambda\) using the OSQP quadratic program optimiser (if cvxpy is installed).

Parameters

solver_args (dict) –

Optional arguments for the solver:

  • verbose (bool): Default value of this input is set to False; when True, the solver prints information to screen.

  • alpha (float): The L1/L2 pentalty blending parameter \(\alpha\). Default: alpha=1.0. If alpha=0 is desired, use least_squares.

  • path (bool): Computes the full regularisation path if True, otherwise solves for a single \(\lambda\) value. Default: path=True.

  • max_iter (int): Max number of iterations for the coordinate descent solver (if path=True). Default: max_iter=100.

  • n_lambdas (int): Number of lambda values along the regularisation path (if path=True). Default: n_lambdas=100.

  • lambda_eps (float): Minimum lambda value, as a fraction of lambda_max (if path=True). Default: lambda_eps=1e-3.

  • lambda_max (float): Maximum lambda value (if path=True). If not None this overrides the automatically calculated value. Default: lambda_max=None.

  • tol (float): Convergence tolerance criterion for coordinate descent solver (if path=True). Default: tol=1e-6.

  • crit (str): Information criterion to select optimal lambda from regularisation path. 'AIC' is Akaike information criterion, 'BIC' is Bayesian information criterion, 'CV' is 5-fold cross-validated RSS error (if path=True). Default: crit='CV'.

  • lambda (float): The penalty scaling parameter (if path=False). Default: ``lambda=0.01`.

References

  1. Friedman J., Hastie T., Tibshirani R., (2010) Regularization Paths for Generalized Linear Models via Coordinate Descent. Journal of Statistical Software, 33(1), 1-22. Paper

  2. Zou, H., Hastie, T., Tibshirani, R., (2007) On the “degrees of freedom” of the lasso. The Annals of Statistics, 35(5), 2173–2192. Paper

plot_regpath(nplot=None, show=True)[source]#

Plots the regularisation path. See plot.plot_regpath for full description.

class equadratures.solver.huber(solver_args={})[source]#

Huber regression solver. Minimises the Huber loss function.

This function is identical to the least squares (L2) penalty for small residuals (i.e. \(||\mathbf{A}\mathbf{x}-\mathbf{b}||^2\le M\)). But on large residuals (\(||\mathbf{A}\mathbf{x}-\mathbf{b}||^2 > M\)), its penalty is lower (L1) and increases linearly rather than quadratically. It is thus more forgiving of outliers.

Parameters

solver_args (dict, optional) –

Optional arguments for the solver:

  • verbose (bool): Default value of this input is set to False; when True, the solver prints information to screen.

  • M (float): The huber threshold. Default M=1.

class equadratures.solver.least_absolute_residual(solver_args={})[source]#

Least absolute residual solver. Minimises the L1 norm (absolute residuals).

Parameters

solver_args (dict) –

Optional arguments for the solver:

  • verbose (bool): Default value of this input is set to False; when True, the solver prints information to screen.

  • opt (str): Specifies the underlying optimisation method to use. scipy for scipy, and osqp to use the OSQP optimiser (requires cvxpy to be installed).

class equadratures.solver.least_squares(solver_args={})[source]#

Ordinary least squares solver.

Parameters

solver_args (dict) –

Optional arguments for the solver:

  • verbose (bool): Default value of this input is set to False; when True, the solver prints information to screen.

class equadratures.solver.minimum_norm(solver_args={})[source]#

Minimum norm solver, used for uniquely or overdetermined systems (i.e. m>=n).

Parameters

solver_args (dict) –

Optional arguments for the solver:

  • verbose (bool): Default value of this input is set to False; when True, the solver prints information to screen.

class equadratures.solver.numerical_integration(solver_args={})[source]#

Numerical integration solver. This solves an orthogonal linear system i.e. simply \(\mathbf{x}=\mathbf{A}^T\mathbf{b}\).

Parameters

solver_args (dict) –

Optional arguments for the solver:

  • verbose (bool): Default value of this input is set to False; when True, the solver prints information to screen.

class equadratures.solver.rvm(solver_args={})[source]#

Relevence vector machine solver.

Parameters

solver_args (dict) –

Optional arguments for the solver:

  • verbose (bool): Default value of this input is set to False; when True, the solver prints information to screen.

  • max_iter (int): Maximum number of iterations. Default is 1000.