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

Tuple[RandomVariable[ndarray], RandomVariable[LinearOperator], RandomVariable[LinearOperator], Dict]

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

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