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
, andrelevance-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
; whenTrue
, 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 (requirescvxpy
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
- 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
; whenTrue
, 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 (requirescvxpy
to be installed).
- Returns
The Solver subclass specified in
method
.- Return type
- 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
; whenTrue
, 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
; whenTrue
, 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}\). Thesolve()
is provided viasolver_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
; whenTrue
, 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 (ifcvxpy
is installed).- Parameters
solver_args (dict) –
Optional arguments for the solver:
verbose (bool): Default value of this input is set to
False
; whenTrue
, the solver prints information to screen.alpha (float): The L1/L2 pentalty blending parameter \(\alpha\). Default:
alpha=1.0
. Ifalpha=0
is desired, useleast_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
(ifpath=True
). Default:lambda_eps=1e-3
.lambda_max (float): Maximum lambda value (if
path=True
). If notNone
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 (ifpath=True
). Default:crit='CV'
.lambda (float): The penalty scaling parameter (if
path=False
). Default: ``lambda=0.01`.
References
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
Zou, H., Hastie, T., Tibshirani, R., (2007) On the “degrees of freedom” of the lasso. The Annals of Statistics, 35(5), 2173–2192. Paper
- 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
; whenTrue
, 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
; whenTrue
, 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 (requirescvxpy
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
; whenTrue
, 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
; whenTrue
, 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
; whenTrue
, 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
; whenTrue
, the solver prints information to screen.max_iter (int): Maximum number of iterations. Default is 1000.