problinsolve

probnum.linalg.problinsolve(A, b, A0=None, Ainv0=None, x0=None, assume_A='sympos', maxiter=None, atol=1e-06, rtol=1e-06, callback=None, **kwargs)

Solve the linear system \(A x = b\) in a Bayesian framework.

Probabilistic linear solvers infer solutions to problems of the form

\[Ax=b,\]

where \(A \in \mathbb{R}^{n \times n}\) and \(b \in \mathbb{R}^{n}\). They return a probability measure which quantifies uncertainty in the output arising from finite computational resources or stochastic input. This solver can take prior information either on the linear operator \(A\) or its inverse \(H=A^{ -1}\) in the form of a random variable A0 or Ainv0 and outputs a posterior belief about \(A\) or \(H\). This code implements the method described in Wenger et al. 1 based on the work in Hennig et al. 2.

Parameters

  • A (Union[LinearOperatorLike, 'randvars.RandomVariable[LinearOperatorLike]']) – shape=(n, n) – A square linear operator (or matrix). Only matrix-vector products \(v \mapsto Av\) are used internally.

  • b (Union[np.ndarray, 'randvars.RandomVariable[np.ndarray]']) – shape=(n, ) or (n, nrhs) – Right-hand side vector, matrix or random variable in \(A x = b\).

  • A0 (Optional[Union[LinearOperatorLike, 'randvars.RandomVariable[LinearOperatorLike]']]) – shape=(n, n) – A square matrix, linear operator or random variable representing the prior belief about the linear operator \(A\).

  • Ainv0 (Optional[Union[LinearOperatorLike, 'randvars.RandomVariable[LinearOperatorLike]']]) – shape=(n, n) – A square matrix, linear operator or random variable representing the prior belief about the inverse \(H=A^{-1}\). This can be viewed as a preconditioner.

  • x0 (Optional[Union[np.ndarray, 'randvars.RandomVariable[np.ndarray]']]) – shape=(n, ) or (n, nrhs) – Prior belief for the solution of the linear system. Will be ignored if Ainv0 is given.

  • assume_A (str) –

    Assumptions on the linear operator which can influence solver choice and behavior. The available options are (combinations of)

    generic matrix

    gen

    symmetric

    sym

    positive definite

    pos

    (additive) noise

    noise

  • maxiter (Optional[int]) – Maximum number of iterations. Defaults to \(10n\), where \(n\) is the dimension of \(A\).

  • atol (float) – Absolute convergence tolerance.

  • rtol (float) – Relative convergence tolerance.

  • callback (Optional[Callable]) – User-supplied function called after each iteration of the linear solver. It is called as callback(xk, Ak, Ainvk, sk, yk, alphak, resid, **kwargs) and can be used to return quantities from the iteration. Note that depending on the function supplied, this can slow down the solver considerably.

  • kwargs – Optional keyword arguments passed onto the solver iteration.

Returns

  • x – Approximate solution \(x\) to the linear system. Shape of the return matches the shape of b.

  • A – Posterior belief over the linear operator.

  • Ainv – Posterior belief over the linear operator inverse \(H=A^{-1}\).

  • info – Information on convergence of the solver.

Raises

  • ValueError – If size mismatches detected or input matrices are not square.

  • LinAlgError – If the matrix A is singular.

  • LinAlgWarning – If an ill-conditioned input A is detected.

Return type

Tuple[‘randvars.RandomVariable[np.ndarray]’, ‘randvars.RandomVariable[linops.LinearOperator]’, ‘randvars.RandomVariable[linops.LinearOperator]’, Dict]

Notes

For a specific class of priors the posterior mean of \(x_k=Hb\) coincides with the iterates of the conjugate gradient method. The matrix-based view taken here recovers the solution-based inference of bayescg() 3.

References

1

Wenger, J. and Hennig, P., Probabilistic Linear Solvers for Machine Learning, Advances in Neural Information Processing Systems (NeurIPS), 2020

2

Hennig, P., Probabilistic Interpretation of Linear Solvers, SIAM Journal on Optimization, 2015, 25, 234-260

3

Bartels, S. et al., Probabilistic Linear Solvers: A Unifying View, Statistics and Computing, 2019

See also

bayescg

Solve linear systems with prior information on the solution.

Examples

>>> import numpy as np
>>> np.random.seed(1)
>>> n = 20
>>> A = np.random.rand(n, n)
>>> A = 0.5 * (A + A.T) + 5 * np.eye(n)
>>> b = np.random.rand(n)
>>> x, A, Ainv, info = problinsolve(A=A, b=b)
>>> print(info["iter"])
9