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


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.

Return type

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


  • 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.

  • 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.


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.



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


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


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

See also


Solve linear systems with prior information on the solution.


>>> 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"])