# 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[ndarray, spmatrix, LinearOperator, RandomVariable[Union[ndarray, spmatrix, ForwardRef]]]) – shape=(n, n) – A square linear operator (or matrix). Only matrix-vector products $$v \mapsto Av$$ are used internally.

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

• A0 (Union[ndarray, spmatrix, LinearOperator, RandomVariable[Union[ndarray, spmatrix, ForwardRef]], None]) – shape=(n, n) – A square matrix, linear operator or random variable representing the prior belief about the linear operator $$A$$.

• Ainv0 (Union[ndarray, spmatrix, LinearOperator, RandomVariable[Union[ndarray, spmatrix, ForwardRef]], None]) – 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 (Union[ndarray, RandomVariable[ndarray], None]) – 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) – Optional keyword arguments passed onto the solver iteration.

Return type
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.

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

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