# 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
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