Showing 22 of 52 files from the diff.
Other files ignored by Codecov
setup.cfg has changed.

@@ -22,6 +22,8 @@
Loading
22 22
23 23
_ValueType = TypeVar("ValueType")
24 24
25 +
# pylint: disable="too-many-lines"
26 +
25 27
26 28
class RandomVariable(Generic[_ValueType]):
27 29
    """Random variables represent uncertainty about a value.

@@ -2,21 +2,18 @@
Loading
2 2
3 3
This package implements probabilistic numerical methods for the solution
4 4
of problems arising in linear algebra, such as the solution of linear
5 -
systems.
5 +
systems :math:`Ax=b`.
6 6
"""
7 -
from probnum.linalg.linearsolvers import *
7 +
from probnum.linalg._problinsolve import bayescg, problinsolve
8 +
9 +
from . import solvers
8 10
9 11
# Public classes and functions. Order is reflected in documentation.
10 12
__all__ = [
11 13
    "problinsolve",
12 14
    "bayescg",
13 -
    "ProbabilisticLinearSolver",
14 -
    "MatrixBasedSolver",
15 -
    "AsymmetricMatrixBasedSolver",
16 -
    "SymmetricMatrixBasedSolver",
17 -
    "SolutionBasedSolver",
18 15
]
19 16
20 17
# Set correct module paths. Corrects links and module paths in documentation.
21 -
ProbabilisticLinearSolver.__module__ = "probnum.linalg"
22 -
MatrixBasedSolver.__module__ = "probnum.linalg"
18 +
problinsolve.__module__ = "probnum.linalg"
19 +
bayescg.__module__ = "probnum.linalg"
23 20
imilarity index 86%
24 21
ename from src/probnum/linalg/linearsolvers/linearsolvers.py
25 22
ename to src/probnum/linalg/_problinsolve.py

@@ -0,0 +1,9 @@
Loading
1 +
"""Belief classes for the quantities of interest of a linear system."""
2 +
3 +
from ._linear_system_belief import LinearSystemBelief
4 +
5 +
# Public classes and functions. Order is reflected in documentation.
6 +
__all__ = ["LinearSystemBelief"]
7 +
8 +
# Set correct module paths. Corrects links and module paths in documentation.
9 +
LinearSystemBelief.__module__ = "probnum.linalg.solvers.beliefs"

@@ -0,0 +1,43 @@
Loading
1 +
"""Probabilistic linear solvers.
2 +
3 +
Iterative probabilistic numerical methods solving linear systems :math:`Ax = b`.
4 +
"""
5 +
6 +
7 +
class ProbabilisticLinearSolver:
8 +
    r"""Compose a custom probabilistic linear solver.
9 +
10 +
    Class implementing probabilistic linear solvers. Such (iterative) solvers infer
11 +
    solutions to problems of the form
12 +
13 +
    .. math:: Ax=b,
14 +
15 +
    where :math:`A \in \mathbb{R}^{n \times n}` and :math:`b \in \mathbb{R}^{n}`.
16 +
    They return a probability measure which quantifies uncertainty in the output arising
17 +
    from finite computational resources or stochastic input. This class unifies and
18 +
    generalizes probabilistic linear solvers as described in the literature. [1]_ [2]_
19 +
    [3]_ [4]_
20 +
21 +
    Parameters
22 +
    ----------
23 +
24 +
25 +
    References
26 +
    ----------
27 +
    .. [1] Hennig, P., Probabilistic Interpretation of Linear Solvers, *SIAM Journal on
28 +
       Optimization*, 2015, 25, 234-260
29 +
    .. [2] Cockayne, J. et al., A Bayesian Conjugate Gradient Method, *Bayesian
30 +
       Analysis*, 2019, 14, 937-1012
31 +
    .. [3] Bartels, S. et al., Probabilistic Linear Solvers: A Unifying View,
32 +
       *Statistics and Computing*, 2019
33 +
    .. [4] Wenger, J. and Hennig, P., Probabilistic Linear Solvers for Machine Learning,
34 +
       *Advances in Neural Information Processing Systems (NeurIPS)*, 2020
35 +
36 +
    See Also
37 +
    --------
38 +
    problinsolve : Solve linear systems in a Bayesian framework.
39 +
    bayescg : Solve linear systems with prior information on the solution.
40 +
41 +
    Examples
42 +
    --------
43 +
    """

@@ -13,6 +13,8 @@
Loading
13 13
    "LinearOperator", ScalarArgType, np.ndarray, scipy.sparse.spmatrix
14 14
]
15 15
16 +
# pylint: disable="too-many-lines"
17 +
16 18
17 19
class LinearOperator:
18 20
    r"""Composite base class for finite-dimensional linear operators.

@@ -182,8 +182,8 @@
Loading
182 182
    >>> import numpy as np
183 183
    >>> A = np.eye(3)
184 184
    >>> b = np.arange(3)
185 -
    >>> lin_sys = LinearSystem(A, b)
186 -
    >>> lin_sys
185 +
    >>> linsys = LinearSystem(A, b)
186 +
    >>> linsys
187 187
    LinearSystem(A=array([[1., 0., 0.],
188 188
           [0., 1., 0.],
189 189
           [0., 0., 1.]]), b=array([0, 1, 2]), solution=None)

@@ -0,0 +1,154 @@
Loading
1 +
"""Linear system belief.
2 +
3 +
Class defining a belief about the quantities of interest of a linear
4 +
system such as its solution or the matrix inverse and any associated
5 +
hyperparameters.
6 +
"""
7 +
8 +
from typing import Mapping, Optional
9 +
10 +
from probnum import randvars
11 +
12 +
try:
13 +
    # functools.cached_property is only available in Python >=3.8
14 +
    from functools import cached_property
15 +
except ImportError:
16 +
    from cached_property import cached_property
17 +
18 +
# pylint: disable="invalid-name"
19 +
20 +
21 +
class LinearSystemBelief:
22 +
    r"""Belief about quantities of interest of a linear system.
23 +
24 +
    Random variables :math:`(\mathsf{x}, \mathsf{A}, \mathsf{H}, \mathsf{b})`
25 +
    modelling the solution :math:`x`, the system matrix :math:`A`, its (pseudo-)inverse
26 +
    :math:`H=A^{-1}` and the right hand side :math:`b` of a linear system :math:`Ax=b`, as well as any associated hyperparameters.
27 +
28 +
    For instantiation either a belief about the solution or the inverse and right hand side must be provided. Note that if both are specified, their consistency is not checked and depending on the algorithm either may be used.
29 +
30 +
    Parameters
31 +
    ----------
32 +
    x :
33 +
        Belief about the solution.
34 +
    Ainv :
35 +
        Belief about the (pseudo-)inverse of the system matrix.
36 +
    A :
37 +
        Belief about the system matrix.
38 +
    b :
39 +
        Belief about the right hand side.
40 +
    hyperparams :
41 +
        Hyperparameters of the belief.
42 +
    """
43 +
44 +
    def __init__(
45 +
        self,
46 +
        x: Optional[randvars.RandomVariable] = None,
47 +
        Ainv: Optional[randvars.RandomVariable] = None,
48 +
        A: Optional[randvars.RandomVariable] = None,
49 +
        b: Optional[randvars.RandomVariable] = None,
50 +
        hyperparams: Optional[Mapping[str, randvars.RandomVariable]] = None,
51 +
    ):
52 +
53 +
        if x is None and Ainv is None:
54 +
            raise TypeError(
55 +
                "Belief over the solution x and the inverse Ainv cannot both be None."
56 +
            )
57 +
58 +
        # Check shapes and their compatibility
59 +
        def dim_mismatch_error(**kwargs):
60 +
            argnames = list(kwargs.keys())
61 +
            return ValueError(
62 +
                f"Dimension mismatch. The shapes of {argnames[0]} : {kwargs[argnames[0]].shape} "
63 +
                f"and {argnames[1]} : {kwargs[argnames[1]].shape} must match."
64 +
            )
65 +
66 +
        if x is not None:
67 +
            if x.ndim > 2 or x.ndim < 1:
68 +
                raise ValueError(
69 +
                    f"Belief over solution must have either one or two dimensions, but has {x.ndim}."
70 +
                )
71 +
            if A is not None:
72 +
                if A.shape[1] != x.shape[0]:
73 +
                    raise dim_mismatch_error(A=A, x=x)
74 +
75 +
            if x.ndim > 1:
76 +
                if x.shape[1] != b.shape[1]:
77 +
                    raise dim_mismatch_error(x=x, b=b)
78 +
            elif b is not None:
79 +
                if b.ndim > 1:
80 +
                    raise dim_mismatch_error(x=x, b=b)
81 +
82 +
        if Ainv is not None:
83 +
            if Ainv.ndim != 2:
84 +
                raise ValueError(
85 +
                    f"Belief over the inverse system matrix may have at most two dimensions, but has {A.ndim}."
86 +
                )
87 +
            if A is not None:
88 +
                if A.shape != Ainv.shape:
89 +
                    raise dim_mismatch_error(A=A, Ainv=Ainv)
90 +
91 +
        if A is not None:
92 +
            if A.ndim != 2:
93 +
                raise ValueError(
94 +
                    f"Belief over the system matrix may have at most two dimensions, but has {A.ndim}."
95 +
                )
96 +
            if b is not None:
97 +
                if A.shape[0] != b.shape[0]:
98 +
                    raise dim_mismatch_error(A=A, b=b)
99 +
100 +
        if b is not None:
101 +
            if b.ndim > 2 or b.ndim < 1:
102 +
                raise ValueError(
103 +
                    f"Belief over right-hand-side may have either one or two dimensions but has {b.ndim}."
104 +
                )
105 +
106 +
        self._x = x
107 +
        self._A = A
108 +
        self._Ainv = Ainv
109 +
        self._b = b
110 +
        if hyperparams is None:
111 +
            hyperparams = {}
112 +
        self._hyperparams = hyperparams
113 +
114 +
    def hyperparameter(self, key: str) -> randvars.RandomVariable:
115 +
        """Hyperparameter of the linear system belief.
116 +
117 +
        Parameters
118 +
        ----------
119 +
        key :
120 +
            Hyperparameter key.
121 +
        """
122 +
        return self._hyperparams[key]
123 +
124 +
    @cached_property
125 +
    def x(self) -> randvars.RandomVariable:
126 +
        """Belief about the solution."""
127 +
        if self._x is None:
128 +
            return self._induced_x()
129 +
        else:
130 +
            return self._x
131 +
132 +
    @property
133 +
    def A(self) -> randvars.RandomVariable:
134 +
        """Belief about the system matrix."""
135 +
        return self._A
136 +
137 +
    @property
138 +
    def Ainv(self) -> Optional[randvars.RandomVariable]:
139 +
        """Belief about the (pseudo-)inverse of the system matrix."""
140 +
        return self._Ainv
141 +
142 +
    @property
143 +
    def b(self) -> randvars.RandomVariable:
144 +
        """Belief about the right hand side."""
145 +
        return self._b
146 +
147 +
    def _induced_x(self) -> randvars.RandomVariable:
148 +
        r"""Induced belief about the solution from a belief about the inverse.
149 +
150 +
        Computes the induced belief about the solution given by (an approximation
151 +
        to) the random variable :math:`x=Hb`. This assumes independence between
152 +
        :math:`H` and :math:`b`.
153 +
        """
154 +
        return self.Ainv @ self.b
0 155
imilarity index 70%
1 156
ename from src/probnum/linalg/linearsolvers/matrixbased.py
2 157
ename to src/probnum/linalg/solvers/matrixbased.py

@@ -339,17 +339,17 @@
Loading
339 339
        return np.array([y1, y2])
340 340
341 341
    def df(t, x):
342 -
        x1, x2 = x
342 +
        x1, _ = x
343 343
        y1 = [1, step]
344 344
        y2 = [-g * np.cos(x1) * step, 1]
345 345
        return np.array([y1, y2])
346 346
347 347
    def h(t, x):
348 -
        x1, x2 = x
348 +
        x1, _ = x
349 349
        return np.array([np.sin(x1)])
350 350
351 351
    def dh(t, x):
352 -
        x1, x2 = x
352 +
        x1, _ = x
353 353
        return np.array([[np.cos(x1), 0.0]])
354 354
355 355
    process_noise_cov = (

@@ -0,0 +1,23 @@
Loading
1 +
"""Policy returning randomly drawn standard unit vectors."""
2 +
import numpy as np
3 +
4 +
import probnum  # pylint: disable="unused-import"
5 +
6 +
from . import _linear_solver_policy
7 +
8 +
9 +
class RandomUnitVectorPolicy(_linear_solver_policy.LinearSolverPolicy):
10 +
    r"""Policy returning randomly drawn standard unit vectors.
11 +
12 +
    Draw a standard unit vector :math:`e_i` at random and return it. This policy corresponds to selecting columns of the matrix as observations.
13 +
    """
14 +
15 +
    def __call__(
16 +
        self, solver_state: "probnum.linalg.solvers.ProbabilisticLinearSolverState"
17 +
    ) -> np.ndarray:
18 +
19 +
        n = solver_state.problem.A.shape[1]
20 +
        idx = solver_state.rng.choice(n, 1)
21 +
        action = np.zeros(n)
22 +
        action[idx] = 1.0
23 +
        return action

@@ -8,91 +8,12 @@
Loading
8 8
9 9
import numpy as np
10 10
11 -
import probnum
12 11
from probnum import linops, randvars
13 12
13 +
# pylint: disable="too-many-branches,too-many-lines,too-complex,too-many-statements,redefined-builtin,arguments-differ,abstract-method,unused-argument"
14 14
15 -
class ProbabilisticLinearSolver(abc.ABC):
16 -
    """An abstract base class for probabilistic linear solvers.
17 15
18 -
    This class is designed to be subclassed with new (probabilistic) linear solvers,
19 -
    which implement a ``.solve()`` method. Objects of this type are instantiated in
20 -
    wrapper functions such as :meth:``problinsolve``.
21 -
22 -
    Parameters
23 -
    ----------
24 -
    A : array-like or LinearOperator or RandomVariable, shape=(n,n)
25 -
        A square matrix or linear operator. A prior distribution can be provided as a
26 -
        :class:`~randvars.RandomVariable`. If an array or linear operator is given,
27 -
        a prior distribution is chosen automatically.
28 -
    b : RandomVariable, shape=(n,) or (n, nrhs)
29 -
        Right-hand side vector, matrix or RandomVariable of :math:`A x = b`.
30 -
    """
31 -
32 -
    def __init__(self, A, b):
33 -
        self.A = A
34 -
        self.b = b
35 -
        self.n = A.shape[1]
36 -
37 -
    def has_converged(self, iter, maxiter, **kwargs):
38 -
        """Check convergence of a linear solver.
39 -
40 -
        Evaluates a set of convergence criteria based on its input arguments to decide
41 -
        whether the iteration has converged.
42 -
43 -
        Parameters
44 -
        ----------
45 -
        iter : int
46 -
            Current iteration of solver.
47 -
        maxiter : int
48 -
            Maximum number of iterations
49 -
50 -
        Returns
51 -
        -------
52 -
        has_converged : bool
53 -
            True if the method has converged.
54 -
        convergence_criterion : str
55 -
            Convergence criterion which caused termination.
56 -
        """
57 -
        # maximum iterations
58 -
        if iter >= maxiter:
59 -
            warnings.warn(
60 -
                "Iteration terminated. Solver reached the maximum number of iterations."
61 -
            )
62 -
            return True, "maxiter"
63 -
        else:
64 -
            return False, ""
65 -
66 -
    def solve(self, callback=None, **kwargs):
67 -
        """Solve the linear system :math:`Ax=b`.
68 -
69 -
        Parameters
70 -
        ----------
71 -
        callback : function, optional
72 -
            User-supplied function called after each iteration of the linear solver. It
73 -
            is called as ``callback(xk, Ak, Ainvk, sk, yk, alphak, resid, **kwargs)``
74 -
            and can be used to return quantities from the iteration. Note that depending
75 -
            on the function supplied, this can slow down the solver.
76 -
        kwargs
77 -
            Key-word arguments adjusting the behaviour of the ``solve`` iteration. These
78 -
            are usually convergence criteria.
79 -
80 -
        Returns
81 -
        -------
82 -
        x : RandomVariable, shape=(n,) or (n, nrhs)
83 -
            Approximate solution :math:`x` to the linear system. Shape of the return
84 -
            matches the shape of ``b``.
85 -
        A : RandomVariable, shape=(n,n)
86 -
            Posterior belief over the linear operator.
87 -
        Ainv : RandomVariable, shape=(n,n)
88 -
            Posterior belief over the linear operator inverse :math:`H=A^{-1}`.
89 -
        info : dict
90 -
            Information on convergence of the solver.
91 -
        """
92 -
        raise NotImplementedError
93 -
94 -
95 -
class MatrixBasedSolver(ProbabilisticLinearSolver, abc.ABC):
16 +
class MatrixBasedSolver(abc.ABC):
96 17
    """Abstract class for matrix-based probabilistic linear solvers.
97 18
98 19
    Parameters
@@ -110,7 +31,9 @@
Loading
110 31
111 32
    def __init__(self, A, b, x0=None):
112 33
        self.x0 = x0
113 -
        super().__init__(A=A, b=b)
34 +
        self.A = A
35 +
        self.b = b
36 +
        self.n = A.shape[1]
114 37
115 38
    def _get_prior_params(self, A0, Ainv0, x0, b):
116 39
        """Find the parameters of the prior distribution.
@@ -201,31 +124,60 @@
Loading
201 124
        return A0_mean, Ainv0_mean
202 125
203 126
    def has_converged(self, iter, maxiter, **kwargs):
204 -
        raise NotImplementedError
205 -
206 -
    def solve(self, callback=None, maxiter=None, atol=None):
207 -
        raise NotImplementedError
127 +
        """Check convergence of a linear solver.
208 128
129 +
        Evaluates a set of convergence criteria based on its input arguments to decide
130 +
        whether the iteration has converged.
209 131
210 -
class AsymmetricMatrixBasedSolver(ProbabilisticLinearSolver):
211 -
    """Asymmetric matrix-based probabilistic linear solver.
132 +
        Parameters
133 +
        ----------
134 +
        iter : int
135 +
            Current iteration of solver.
136 +
        maxiter : int
137 +
            Maximum number of iterations
212 138
213 -
    Parameters
214 -
    ----------
215 -
    A : array-like or LinearOperator or RandomVariable, shape=(n,n)
216 -
        The square matrix or linear operator of the linear system.
217 -
    b : array_like, shape=(n,) or (n, nrhs)
218 -
        Right-hand side vector or matrix in :math:`A x = b`.
219 -
    """
139 +
        Returns
140 +
        -------
141 +
        has_converged : bool
142 +
            True if the method has converged.
143 +
        convergence_criterion : str
144 +
            Convergence criterion which caused termination.
145 +
        """
146 +
        # maximum iterations
147 +
        if iter >= maxiter:
148 +
            warnings.warn(
149 +
                "Iteration terminated. Solver reached the maximum number of iterations."
150 +
            )
151 +
            return True, "maxiter"
152 +
        else:
153 +
            return False, ""
220 154
221 -
    def __init__(self, A, b, x0):
222 -
        self.x0 = x0
223 -
        super().__init__(A=A, b=b)
155 +
    def solve(self, callback=None, **kwargs):
156 +
        """Solve the linear system :math:`Ax=b`.
224 157
225 -
    def has_converged(self, iter, maxiter, **kwargs):
226 -
        raise NotImplementedError
158 +
        Parameters
159 +
        ----------
160 +
        callback : function, optional
161 +
            User-supplied function called after each iteration of the linear solver. It
162 +
            is called as ``callback(xk, Ak, Ainvk, sk, yk, alphak, resid, **kwargs)``
163 +
            and can be used to return quantities from the iteration. Note that depending
164 +
            on the function supplied, this can slow down the solver.
165 +
        kwargs
166 +
            Key-word arguments adjusting the behaviour of the ``solve`` iteration. These
167 +
            are usually convergence criteria.
227 168
228 -
    def solve(self, callback=None, maxiter=None, atol=None):
169 +
        Returns
170 +
        -------
171 +
        x : RandomVariable, shape=(n,) or (n, nrhs)
172 +
            Approximate solution :math:`x` to the linear system. Shape of the return
173 +
            matches the shape of ``b``.
174 +
        A : RandomVariable, shape=(n,n)
175 +
            Posterior belief over the linear operator.
176 +
        Ainv : RandomVariable, shape=(n,n)
177 +
            Posterior belief over the linear operator inverse :math:`H=A^{-1}`.
178 +
        info : dict
179 +
            Information on convergence of the solver.
180 +
        """
229 181
        raise NotImplementedError
230 182
231 183
@@ -1017,345 +969,3 @@
Loading
1017 969
        }
1018 970
1019 971
        return x, A, Ainv, info
1020 -
1021 -
1022 -
class NoisySymmetricMatrixBasedSolver(MatrixBasedSolver):
1023 -
    """Solver iteration of the noisy symmetric probabilistic linear solver.
1024 -
1025 -
    Implements the solve iteration of the symmetric matrix-based probabilistic linear
1026 -
    solver taking into account noisy matrix-vector products :math:`y_k = (A + E_k)s_k`
1027 -
    as described in [1]_ and [2]_.
1028 -
1029 -
    Parameters
1030 -
    ----------
1031 -
    A : LinearOperator or RandomVariable, shape=(n,n)
1032 -
        The square matrix or linear operator of the linear system.
1033 -
    b : array_like, shape=(n,) or (n, nrhs)
1034 -
        Right-hand side vector or matrix in :math:`A x = b`.
1035 -
    A0 : array-like or LinearOperator or RandomVariable, shape=(n, n), optional
1036 -
        A square matrix, linear operator or random variable representing the prior
1037 -
        belief over the linear operator :math:`A`. If an array or linear operator is
1038 -
        given, a prior distribution is chosen automatically.
1039 -
    Ainv0 : array-like or LinearOperator or RandomVariable, shape=(n,n), optional
1040 -
        A square matrix, linear operator or random variable representing the prior
1041 -
        belief over the inverse :math:`H=A^{-1}`. This can be viewed as taking the form
1042 -
        of a pre-conditioner. If an array or linear operator is given, a prior
1043 -
        distribution is chosen automatically.
1044 -
    x0 : array-like, or RandomVariable, shape=(n,) or (n, nrhs)
1045 -
        Optional. Prior belief for the solution of the linear system. Will be ignored if
1046 -
        ``Ainv0`` is given.
1047 -
1048 -
    Returns
1049 -
    -------
1050 -
    A : RandomVariable
1051 -
        Posterior belief over the linear operator.
1052 -
    Ainv : RandomVariable
1053 -
        Posterior belief over the inverse linear operator.
1054 -
    x : RandomVariable
1055 -
        Posterior belief over the solution of the linear system.
1056 -
    info : dict
1057 -
        Information about convergence and the solution.
1058 -
1059 -
    References
1060 -
    ----------
1061 -
    .. [1] Wenger, J., de Roos, F. and Hennig, P., Probabilistic Solution of Noisy
1062 -
       Linear Systems, 2020
1063 -
    .. [2] Hennig, P., Probabilistic Interpretation of Linear Solvers, *SIAM Journal on
1064 -
       Optimization*, 2015, 25, 234-260
1065 -
1066 -
    See Also
1067 -
    --------
1068 -
    SymmetricMatrixBasedSolver :
1069 -
        Class implementing the symmetric probabilistic linear solver.
1070 -
    """
1071 -
1072 -
    def __init__(self, A, b, A0=None, Ainv0=None, x0=None):
1073 -
1074 -
        # Transform right hand side to random variable
1075 -
        if not isinstance(b, randvars.RandomVariable):
1076 -
            _b = probnum.asrandvar(b)
1077 -
        else:
1078 -
            _b = b
1079 -
1080 -
        super().__init__(A=A, b=_b, x0=x0)
1081 -
1082 -
        # Get or initialize prior parameters
1083 -
        (
1084 -
            A0_mean,
1085 -
            A0_covfactor,
1086 -
            Ainv0_mean,
1087 -
            Ainv0_covfactor,
1088 -
            b_mean,
1089 -
        ) = self._get_prior_params(A0=A0, Ainv0=Ainv0, x0=x0, b=_b)
1090 -
1091 -
        # Matrix prior parameters
1092 -
        self.A0_mean = linops.aslinop(A0_mean)
1093 -
        self.A_mean = linops.aslinop(A0_mean)
1094 -
        self.A0_covfactor = A0_covfactor
1095 -
        self.Ainv0_mean = linops.aslinop(Ainv0_mean)
1096 -
        self.Ainv_mean = linops.aslinop(Ainv0_mean)
1097 -
        self.Ainv0_covfactor = Ainv0_covfactor
1098 -
        self.b_mean = b_mean
1099 -
1100 -
        # Induced distribution on x via Ainv
1101 -
        # Exp = x = A^-1 b, Cov = 1/2 (W b'Wb + Wbb'W)
1102 -
        Wb = Ainv0_covfactor @ self.b_mean
1103 -
        bWb = np.squeeze(Wb.T @ self.b_mean)
1104 -
1105 -
        def _matmul(x):
1106 -
            return 0.5 * (bWb * Ainv0_covfactor @ x + Wb @ (Wb.T @ x))
1107 -
1108 -
        self.x_cov = linops.LinearOperator(
1109 -
            shape=(self.n, self.n),
1110 -
            dtype=np.result_type(bWb.dtype, Ainv0_covfactor.dtype, Wb.dtype),
1111 -
            matmul=_matmul,
1112 -
        )
1113 -
        if isinstance(x0, np.ndarray):
1114 -
            self.x_mean = x0
1115 -
        elif x0 is None:
1116 -
            self.x_mean = Ainv0_mean @ self.b_mean
1117 -
        else:
1118 -
            raise NotImplementedError
1119 -
        self.x0 = self.x_mean
1120 -
1121 -
    def _get_prior_params(self, A0, Ainv0, x0, b):
1122 -
        """Get the parameters of the matrix priors on A and H.
1123 -
1124 -
        Retrieves and / or initializes prior parameters of ``A0`` and ``Ainv0``.
1125 -
1126 -
        Parameters
1127 -
        ----------
1128 -
        A0 : array-like or LinearOperator or RandomVariable, shape=(n,n), optional
1129 -
            A square matrix, linear operator or random variable representing the prior
1130 -
            belief over the linear operator :math:`A`. If an array or linear operator is
1131 -
            given, a prior distribution is chosen automatically.
1132 -
        Ainv0 : array-like or LinearOperator or RandomVariable, shape=(n,n), optional
1133 -
            A square matrix, linear operator or random variable representing the prior
1134 -
            belief over the inverse :math:`H=A^{-1}`. This can be viewed as taking the
1135 -
            form of a pre-conditioner. If an array or linear operator is given, a prior
1136 -
            distribution is chosen automatically.
1137 -
        x0 : array-like, or RandomVariable, shape=(n,)
1138 -
            Optional. Prior belief for the solution of the linear system. Will be
1139 -
            ignored if ``A0`` or ``Ainv0`` is given.
1140 -
        b : RandomVariable, shape=(n,) or (n, nrhs)
1141 -
            Right-hand side random variable `b` in :math:`A x = b`.
1142 -
1143 -
        Returns
1144 -
        -------
1145 -
        A0_mean : array-like or LinearOperator, shape=(n,n)
1146 -
            Prior mean of the linear operator :math:`A`.
1147 -
        A0_covfactor : array-like or LinearOperator, shape=(n,n)
1148 -
            Factor :math:`W^A` of the symmetric Kronecker product prior covariance
1149 -
            :math:`W^A \\otimes_s W^A` of :math:`A`.
1150 -
        Ainv0_mean : array-like or LinearOperator, shape=(n,n)
1151 -
            Prior mean of the linear operator :math:`H`.
1152 -
        Ainv0_covfactor : array-like or LinearOperator, shape=(n,n)
1153 -
            Factor :math:`W^H` of the symmetric Kronecker product prior covariance
1154 -
            :math:`W^H \\otimes_s W^H` of :math:`H`.
1155 -
        b_mean : array-like, shape=(n,nrhs)
1156 -
            Prior mean of the right hand side :math:`b`.
1157 -
        """
1158 -
1159 -
        # Right hand side mean
1160 -
        b_mean = b.sample(1)  # TODO: build prior model for rhs and change to b.mean
1161 -
1162 -
        # No matrix priors specified
1163 -
        if A0 is None and Ainv0 is None:
1164 -
            # No prior information given
1165 -
            if x0 is None:
1166 -
                Ainv0_mean = linops.Identity(shape=self.n)
1167 -
                Ainv0_covfactor = linops.Identity(shape=self.n)
1168 -
                # Standard normal covariance
1169 -
                A0_mean = linops.Identity(shape=self.n)
1170 -
                A0_covfactor = linops.Identity(shape=self.n)
1171 -
                # TODO: should this be a sample from A to achieve symm. posterior
1172 -
                # correspondence?
1173 -
                return A0_mean, A0_covfactor, Ainv0_mean, Ainv0_covfactor, b_mean
1174 -
            # Construct matrix priors from initial guess x0
1175 -
            elif isinstance(x0, np.ndarray):
1176 -
                # Sample from linear operator for prior construction
1177 -
                if isinstance(self.A, randvars.RandomVariable):
1178 -
                    _A = self.A.sample([1])[0]
1179 -
                else:
1180 -
                    _A = self.A
1181 -
                A0_mean, Ainv0_mean = self._construct_symmetric_matrix_prior_means(
1182 -
                    A=_A, x0=x0, b=b_mean
1183 -
                )
1184 -
                Ainv0_covfactor = Ainv0_mean
1185 -
                # Standard normal covariance
1186 -
                A0_covfactor = linops.Identity(shape=self.n)
1187 -
                # TODO: should this be a sample from A to achieve symm. posterior
1188 -
                # correspondence?
1189 -
                return A0_mean, A0_covfactor, Ainv0_mean, Ainv0_covfactor, b_mean
1190 -
            elif isinstance(x0, randvars.RandomVariable):
1191 -
                raise NotImplementedError
1192 -
1193 -
        # Prior on Ainv specified
1194 -
        if not isinstance(A0, randvars.RandomVariable) and Ainv0 is not None:
1195 -
            if isinstance(Ainv0, randvars.RandomVariable):
1196 -
                Ainv0_mean = Ainv0.mean
1197 -
                Ainv0_covfactor = Ainv0.cov.A
1198 -
            else:
1199 -
                Ainv0_mean = Ainv0
1200 -
                Ainv0_covfactor = Ainv0  # Symmetric posterior correspondence
1201 -
            try:
1202 -
                if A0 is not None:
1203 -
                    A0_mean = A0
1204 -
                elif isinstance(Ainv0, randvars.RandomVariable):
1205 -
                    A0_mean = Ainv0.mean.inv()
1206 -
                else:
1207 -
                    A0_mean = Ainv0.inv()
1208 -
            except AttributeError:
1209 -
                warnings.warn(
1210 -
                    "Prior specified only for Ainv. Inverting prior mean naively. "
1211 -
                    "This operation is computationally costly! Specify an inverse "
1212 -
                    "prior (mean) instead."
1213 -
                )
1214 -
                A0_mean = np.linalg.inv(Ainv0.mean)
1215 -
            except NotImplementedError:
1216 -
                A0_mean = linops.Identity(self.n)
1217 -
                warnings.warn(
1218 -
                    "Prior specified only for Ainv. Automatic prior mean inversion "
1219 -
                    "not implemented, falling back to standard normal prior."
1220 -
                )
1221 -
            # Standard normal covariance
1222 -
            A0_covfactor = linops.Identity(shape=self.n)
1223 -
            # TODO: should this be a sample from A to achieve symm. posterior
1224 -
            # correspondence?
1225 -
            return A0_mean, A0_covfactor, Ainv0_mean, Ainv0_covfactor, b_mean
1226 -
1227 -
        # Prior on A specified
1228 -
        elif A0 is not None and not isinstance(Ainv0, randvars.RandomVariable):
1229 -
            if isinstance(A0, randvars.RandomVariable):
1230 -
                A0_mean = A0.mean
1231 -
                A0_covfactor = A0.cov.A
1232 -
            else:
1233 -
                A0_mean = A0
1234 -
                A0_covfactor = A0  # Symmetric posterior correspondence
1235 -
            try:
1236 -
                if Ainv0 is not None:
1237 -
                    Ainv0_mean = Ainv0
1238 -
                elif isinstance(A0, randvars.RandomVariable):
1239 -
                    Ainv0_mean = A0.mean.inv()
1240 -
                else:
1241 -
                    Ainv0_mean = A0.inv()
1242 -
            except AttributeError:
1243 -
                warnings.warn(
1244 -
                    "Prior specified only for A. Inverting prior mean naively. "
1245 -
                    "This operation is computationally costly! Specify an inverse "
1246 -
                    "prior (mean) instead."
1247 -
                )
1248 -
                Ainv0_mean = np.linalg.inv(A0.mean)
1249 -
            except NotImplementedError:
1250 -
                Ainv0_mean = linops.Identity(self.n)
1251 -
                warnings.warn(
1252 -
                    "Prior specified only for A. Automatic prior mean inversion "
1253 -
                    "failed, falling back to standard normal prior."
1254 -
                )
1255 -
            # Symmetric posterior correspondence
1256 -
            Ainv0_covfactor = Ainv0_mean
1257 -
            return A0_mean, A0_covfactor, Ainv0_mean, Ainv0_covfactor, b_mean
1258 -
        # Both matrix priors on A and H specified via random variables
1259 -
        elif isinstance(A0, randvars.RandomVariable) and isinstance(
1260 -
            Ainv0, randvars.RandomVariable
1261 -
        ):
1262 -
            A0_mean = A0.mean
1263 -
            A0_covfactor = A0.cov.A
1264 -
            Ainv0_mean = Ainv0.mean
1265 -
            Ainv0_covfactor = Ainv0.cov.A
1266 -
            return A0_mean, A0_covfactor, Ainv0_mean, Ainv0_covfactor, b_mean
1267 -
        else:
1268 -
            raise NotImplementedError
1269 -
1270 -
    def has_converged(self, iter, maxiter, atol=None, rtol=None):
1271 -
        """Check convergence of a linear solver.
1272 -
1273 -
        Evaluates a set of convergence criteria based on its input arguments to decide
1274 -
        whether the iteration has converged.
1275 -
1276 -
        Parameters
1277 -
        ----------
1278 -
        iter : int
1279 -
            Current iteration of solver.
1280 -
        maxiter : int
1281 -
            Maximum number of iterations
1282 -
        atol : float
1283 -
            Absolute tolerance for the uncertainty about the solution estimate. Stops if
1284 -
            :math:`\\sqrt{\\text{tr}(\\Sigma)}  \\leq \\text{atol}`, where
1285 -
            :math:`\\Sigma` is the covariance of the solution :math:`x`.
1286 -
        rtol : float
1287 -
            Relative tolerance for the uncertainty about the solution estimate. Stops if
1288 -
            :math:`\\sqrt{\\text{tr}(\\Sigma)} \\leq \\text{rtol} \\lVert x_i \\rVert`,
1289 -
            where :math:`\\Sigma` is the covariance of the solution :math`x` and
1290 -
            :math:`x_i` its mean.
1291 -
1292 -
        Returns
1293 -
        -------
1294 -
        has_converged : bool
1295 -
            True if the method has converged.
1296 -
        convergence_criterion : str
1297 -
            Convergence criterion which caused termination.
1298 -
        """
1299 -
        # maximum iterations
1300 -
        if iter >= maxiter:
1301 -
            warnings.warn(
1302 -
                "Iteration terminated. Solver reached the maximum number of iterations."
1303 -
            )
1304 -
            return True, "maxiter"
1305 -
        # uncertainty-based
1306 -
        if isinstance(self.x_cov, linops.LinearOperator):
1307 -
            sqrttracecov = np.sqrt(self.x_cov.trace())
1308 -
        else:
1309 -
            sqrttracecov = np.sqrt(np.trace(self.x_cov))
1310 -
        if sqrttracecov <= atol:
1311 -
            return True, "covar_atol"
1312 -
        elif sqrttracecov <= rtol * np.linalg.norm(self.x_mean):
1313 -
            return True, "covar_rtol"
1314 -
        else:
1315 -
            return False, ""
1316 -
1317 -
    def solve(
1318 -
        self,
1319 -
        callback=None,
1320 -
        maxiter=None,
1321 -
        atol=10 ** -6,
1322 -
        rtol=10 ** -6,
1323 -
        noise_scale=None,
1324 -
        **kwargs
1325 -
    ):
1326 -
        """Solve the linear system :math:`Ax=b`.
1327 -
1328 -
        Parameters
1329 -
        ----------
1330 -
        callback : function, optional
1331 -
            User-supplied function called after each iteration of the linear solver. It
1332 -
            is called as ``callback(xk, Ak, Ainvk, sk, yk, alphak, resid, noise_scale)``
1333 -
            and can be used to return quantities from the iteration. Note that depending
1334 -
            on the function supplied, this can slow down the solver.
1335 -
        maxiter : int
1336 -
            Maximum number of iterations
1337 -
        atol : float
1338 -
            Absolute tolerance for the uncertainty about the solution estimate. Stops if
1339 -
            :math:`\\sqrt{\\text{tr}(\\Sigma)}  \\leq \\text{atol}`, where
1340 -
            :math:`\\Sigma` is the covariance of the solution :math:`x`.
1341 -
        rtol : float
1342 -
            Relative tolerance for the uncertainty about the solution estimate. Stops if
1343 -
            :math:`\\sqrt{\\text{tr}(\\Sigma)} \\leq \\text{rtol} \\lVert x_i \\rVert`,
1344 -
            where :math:`\\Sigma` is the covariance of the solution :math`x` and
1345 -
            :math:`x_i` its mean.
1346 -
        noise_scale : float
1347 -
            Assumed (initial) noise scale :math:`\\varepsilon^2`.
1348 -
1349 -
        Returns
1350 -
        -------
1351 -
        x : RandomVariable, shape=(n,) or (n, nrhs)
1352 -
            Approximate solution :math:`x` to the linear system. Shape of the return
1353 -
            matches the shape of ``b``.
1354 -
        A : RandomVariable, shape=(n,n)
1355 -
            Posterior belief over the linear operator.
1356 -
        Ainv : RandomVariable, shape=(n,n)
1357 -
            Posterior belief over the linear operator inverse :math:`H=A^{-1}`.
1358 -
        info : dict
1359 -
            Information on convergence of the solver.
1360 -
        """
1361 -
        raise NotImplementedError

@@ -17,6 +17,7 @@
Loading
17 17
from typing import Iterable, Tuple, Union
18 18
19 19
import numpy as np
20 +
import scipy.sparse
20 21
21 22
########################################################################################
22 23
# API Types
@@ -46,6 +47,13 @@
Loading
46 47
should always be converted into :class:`np.generic` using the function
47 48
:func:`probnum.utils.as_scalar` before further internal processing."""
48 49
50 +
LinearOperatorArgType = Union[
51 +
    np.ndarray,
52 +
    scipy.sparse.spmatrix,
53 +
    "probnum.linops.LinearOperator",
54 +
]
55 +
"""Type of a public API argument for supplying a matrix or finite-dimensional linear operator."""
56 +
49 57
ArrayLikeGetitemArgType = Union[
50 58
    int,
51 59
    slice,
52 60
imilarity index 100%
53 61
ename from tests/test_linalg/test_linearsolvers/__init__.py
54 62
ename to tests/test_linalg/cases/__init__.py

@@ -0,0 +1,26 @@
Loading
1 +
"""Probabilistic linear solvers.
2 +
3 +
Compositional implementation of probabilistic linear solvers. The
4 +
classes and methods in this subpackage allow the creation of custom
5 +
iterative methods for the solution of linear systems.
6 +
"""
7 +
8 +
from probnum.linalg.solvers.matrixbased import (
9 +
    MatrixBasedSolver,
10 +
    SymmetricMatrixBasedSolver,
11 +
)
12 +
13 +
from . import beliefs, policies
14 +
from ._probabilistic_linear_solver import ProbabilisticLinearSolver
15 +
from ._state import ProbabilisticLinearSolverState
16 +
17 +
# Public classes and functions. Order is reflected in documentation.
18 +
__all__ = [
19 +
    "ProbabilisticLinearSolver",
20 +
    "MatrixBasedSolver",
21 +
    "SymmetricMatrixBasedSolver",
22 +
    "ProbabilisticLinearSolverState",
23 +
]
24 +
25 +
# Set correct module paths. Corrects links and module paths in documentation.
26 +
ProbabilisticLinearSolver.__module__ = "probnum.linalg.solvers"

@@ -13,13 +13,13 @@
Loading
13 13
    dim: IntArgType,
14 14
    spectrum: Sequence = None,
15 15
) -> np.ndarray:
16 -
    """Random symmetric positive definite matrix.
16 +
    r"""Random symmetric positive definite matrix.
17 17
18 18
    Constructs a random symmetric positive definite matrix from a given spectrum. An
19 -
    orthogonal matrix :math:`Q` with :math:`\\operatorname{det}(Q)` (a rotation) is
19 +
    orthogonal matrix :math:`Q` with :math:`\operatorname{det}(Q)` (a rotation) is
20 20
    sampled with respect to the Haar measure and the diagonal matrix
21 21
    containing the eigenvalues is rotated accordingly resulting in :math:`A=Q
22 -
    \\operatorname{diag}(\\lambda_1, \\dots, \\lambda_n)Q^\\top`. If no spectrum is
22 +
    \operatorname{diag}(\lambda_1, \dots, \lambda_n)Q^\top`. If no spectrum is
23 23
    provided, one is randomly drawn from a Gamma distribution.
24 24
25 25
    Parameters
@@ -37,8 +37,8 @@
Loading
37 37
38 38
    Examples
39 39
    --------
40 -
    >>> from probnum.problems.zoo.linalg import random_spd_matrix
41 40
    >>> import numpy as np
41 +
    >>> from probnum.problems.zoo.linalg import random_spd_matrix
42 42
    >>> rng = np.random.default_rng(1)
43 43
    >>> mat = random_spd_matrix(rng, dim=5)
44 44
    >>> mat
@@ -97,7 +97,7 @@
Loading
97 97
    chol_entry_min: float = 0.1,
98 98
    chol_entry_max: float = 1.0,
99 99
    format="csr",  # pylint: disable="redefined-builtin"
100 -
) -> np.ndarray:
100 +
) -> scipy.sparse.spmatrix:
101 101
    r"""Random sparse symmetric positive definite matrix.
102 102
103 103
    Constructs a random sparse symmetric positive definite matrix for a given degree
@@ -128,8 +128,8 @@
Loading
128 128
129 129
    Examples
130 130
    --------
131 -
    >>> from probnum.problems.zoo.linalg import random_sparse_spd_matrix
132 131
    >>> import numpy as np
132 +
    >>> from probnum.problems.zoo.linalg import random_sparse_spd_matrix
133 133
    >>> rng = np.random.default_rng(42)
134 134
    >>> sparsemat = random_sparse_spd_matrix(rng, dim=5, density=0.1)
135 135
    >>> sparsemat
@@ -151,7 +151,6 @@
Loading
151 151
    num_nonzero_entries = int(num_off_diag_cholesky * density)
152 152
153 153
    if num_nonzero_entries > 0:
154 -
155 154
        sparse_matrix = scipy.sparse.rand(
156 155
            m=dim,
157 156
            n=dim,

@@ -0,0 +1,31 @@
Loading
1 +
"""Policy returning :math:`A`-conjugate actions."""
2 +
3 +
import numpy as np
4 +
5 +
import probnum  # pylint: disable="unused-import"
6 +
7 +
from . import _linear_solver_policy
8 +
9 +
10 +
class ConjugateGradientPolicy(_linear_solver_policy.LinearSolverPolicy):
11 +
    r"""Policy returning :math:`A`-conjugate actions.
12 +
13 +
    Selects the negative gradient / residual as an initial action :math:`s_0 = b - A x_0` and then successively generates :math:`A`-conjugate actions, i.e. the actions satisfy :math:`s_i^\top A s_j = 0` iff :math:`i \neq j`.
14 +
    """
15 +
16 +
    def __call__(
17 +
        self, solver_state: "probnum.linalg.solvers.ProbabilisticLinearSolverState"
18 +
    ) -> np.ndarray:
19 +
20 +
        action = -solver_state.residual.copy()
21 +
22 +
        if solver_state.step > 0:
23 +
            # A-conjugacy correction (in exact arithmetic)
24 +
            beta = (
25 +
                np.linalg.norm(solver_state.residual)
26 +
                / np.linalg.norm(solver_state.residuals[solver_state.step - 1])
27 +
            ) ** 2
28 +
29 +
            action += beta * solver_state.actions[solver_state.step - 1]
30 +
31 +
        return action

@@ -27,6 +27,9 @@
Loading
27 27
_ValueType = Union[np.floating, np.ndarray, linops.LinearOperator]
28 28
29 29
30 +
# pylint: disable="too-complex"
31 +
32 +
30 33
class Normal(_random_variable.ContinuousRandomVariable[_ValueType]):
31 34
    """Random variable with a normal distribution.
32 35

@@ -12,27 +12,33 @@
Loading
12 12
import numpy as np
13 13
import scipy.sparse
14 14
15 +
import probnum  # pylint: disable=unused-import
15 16
from probnum import linops, randvars, utils
16 -
from probnum.linalg.linearsolvers.matrixbased import (
17 -
    AsymmetricMatrixBasedSolver,
18 -
    NoisySymmetricMatrixBasedSolver,
19 -
    SymmetricMatrixBasedSolver,
20 -
)
21 -
from probnum.linalg.linearsolvers.solutionbased import SolutionBasedSolver
17 +
from probnum.linalg.solvers.matrixbased import SymmetricMatrixBasedSolver
18 +
from probnum.typing import LinearOperatorArgType
22 19
23 -
# Type aliases
24 -
SquareLinOp = Union[
25 -
    np.ndarray, scipy.sparse.spmatrix, linops.LinearOperator, "randvars.RandomVariable"
26 -
]
27 -
RandomVecMat = Union[np.ndarray, "randvars.RandomVariable"]
20 +
# pylint: disable=too-many-branches
28 21
29 22
30 23
def problinsolve(
31 -
    A: SquareLinOp,
32 -
    b: RandomVecMat,
33 -
    A0: Optional[SquareLinOp] = None,
34 -
    Ainv0: Optional[SquareLinOp] = None,
35 -
    x0: Optional[RandomVecMat] = None,
24 +
    A: Union[
25 +
        LinearOperatorArgType,
26 +
        "randvars.RandomVariable[LinearOperatorArgType]",
27 +
    ],
28 +
    b: Union[np.ndarray, "randvars.RandomVariable[np.ndarray]"],
29 +
    A0: Optional[
30 +
        Union[
31 +
            LinearOperatorArgType,
32 +
            "randvars.RandomVariable[LinearOperatorArgType]",
33 +
        ]
34 +
    ] = None,
35 +
    Ainv0: Optional[
36 +
        Union[
37 +
            LinearOperatorArgType,
38 +
            "randvars.RandomVariable[LinearOperatorArgType]",
39 +
        ]
40 +
    ] = None,
41 +
    x0: Optional[Union[np.ndarray, "randvars.RandomVariable[np.ndarray]"]] = None,
36 42
    assume_A: str = "sympos",
37 43
    maxiter: Optional[int] = None,
38 44
    atol: float = 10 ** -6,
@@ -40,45 +46,40 @@
Loading
40 46
    callback: Optional[Callable] = None,
41 47
    **kwargs
42 48
) -> Tuple[
43 -
    "randvars.RandomVariable",
44 -
    "randvars.RandomVariable",
45 -
    "randvars.RandomVariable",
49 +
    "randvars.RandomVariable[np.ndarray]",
50 +
    "randvars.RandomVariable[linops.LinearOperator]",
51 +
    "randvars.RandomVariable[linops.LinearOperator]",
46 52
    Dict,
47 53
]:
48 -
    """Infer a solution to the linear system :math:`A x = b` in a Bayesian framework.
54 +
    r"""Solve the linear system :math:`A x = b` in a Bayesian framework.
49 55
50 56
    Probabilistic linear solvers infer solutions to problems of the form
51 57
52 58
    .. math:: Ax=b,
53 59
54 -
    where :math:`A \\in \\mathbb{R}^{n \\times n}` and :math:`b \\in \\mathbb{R}^{n}`.
60 +
    where :math:`A \in \mathbb{R}^{n \times n}` and :math:`b \in \mathbb{R}^{n}`.
55 61
    They return a probability measure which quantifies uncertainty in the output arising
56 -
    from finite computational resources. This solver can take prior information either
57 -
    on the linear operator :math:`A` or its inverse :math:`H=A^{-1}` in the form of a
58 -
    random variable ``A0`` or ``Ainv0`` and outputs a posterior belief over :math:`A` or
59 -
    :math:`H`. This code implements the method described in Wenger et al. [1]_ based on
60 -
    the work in Hennig et al. [2]_.
62 +
    from finite computational resources or stochastic input. This solver can take prior
63 +
    information either on the linear operator :math:`A` or its inverse :math:`H=A^{
64 +
    -1}` in the form of a random variable ``A0`` or ``Ainv0`` and outputs a posterior
65 +
    belief about :math:`A` or :math:`H`. This code implements the method described in
66 +
    Wenger et al. [1]_ based on the work in Hennig et al. [2]_.
61 67
62 68
    Parameters
63 69
    ----------
64 70
    A :
65 71
        *shape=(n, n)* -- A square linear operator (or matrix). Only matrix-vector
66 -
        products :math:`v \\mapsto Av` are used internally.
72 +
        products :math:`v \mapsto Av` are used internally.
67 73
    b :
68 74
        *shape=(n, ) or (n, nrhs)* -- Right-hand side vector, matrix or random
69 -
        variable in :math:`A x = b`. For multiple right hand sides, ``nrhs`` problems
70 -
        are solved sequentially with the posteriors over the matrices acting as priors
71 -
        for subsequent solves. If the right-hand-side is assumed to be noisy, every
72 -
        iteration of the solver samples a realization from ``b``.
75 +
        variable in :math:`A x = b`.
73 76
    A0 :
74 77
        *shape=(n, n)* -- A square matrix, linear operator or random variable
75 -
        representing the prior belief over the linear operator :math:`A`. If an array or
76 -
        linear operator is given, a prior distribution is chosen automatically.
78 +
        representing the prior belief about the linear operator :math:`A`.
77 79
    Ainv0 :
78 80
        *shape=(n, n)* -- A square matrix, linear operator or random variable
79 -
        representing the prior belief over the inverse :math:`H=A^{-1}`. This can be
80 -
        viewed as taking the form of a pre-conditioner. If an array or linear operator
81 -
        is given, a prior distribution is chosen automatically.
81 +
        representing the prior belief about the inverse :math:`H=A^{-1}`. This can be
82 +
        viewed as a preconditioner.
82 83
    x0 :
83 84
        *shape=(n, ) or (n, nrhs)* -- Prior belief for the solution of the linear
84 85
        system. Will be ignored if ``Ainv0`` is given.
@@ -274,26 +275,7 @@
Loading
274 275
    --------
275 276
    problinsolve : Solve linear systems in a Bayesian framework.
276 277
    """
277 -
    # Check linear system for type and dimension mismatch
278 -
    _check_linear_system(A=A, b=b, x0=x0)
279 -
280 -
    # Preprocess linear system
281 -
    A, b, x0 = _preprocess_linear_system(A=A, b=b, x0=x0)
282 -
283 -
    # Set default convergence parameters
284 -
    n = A.shape[0]
285 -
    if maxiter is None:
286 -
        maxiter = n * 10
287 -
288 -
    # Solve linear system
289 -
    x, info = SolutionBasedSolver(A=A, b=b, x0=x0).solve(
290 -
        callback=callback, maxiter=maxiter, atol=atol, rtol=rtol
291 -
    )
292 -
293 -
    # Check result and issue warnings (e.g. singular or ill-conditioned matrix)
294 -
    _postprocess(info=info, A=A)
295 -
296 -
    return x, info
278 +
    raise NotImplementedError
297 279
298 280
299 281
def _check_linear_system(A, b, A0=None, Ainv0=None, x0=None):
@@ -483,18 +465,16 @@
Loading
483 465
484 466
    # Solution-based view
485 467
    if isinstance(x0, randvars.RandomVariable):
486 -
        return SolutionBasedSolver(A=A, b=b, x0=x0)
468 +
        raise NotImplementedError
487 469
    # Matrix-based view
488 470
    else:
489 471
        if "sym" in assume_A and "pos" in assume_A:
490 472
            if "noise" in assume_A:
491 -
                return NoisySymmetricMatrixBasedSolver(
492 -
                    A=A, b=b, x0=x0, A0=A0, Ainv0=Ainv0
493 -
                )
473 +
                raise NotImplementedError
494 474
            else:
495 475
                return SymmetricMatrixBasedSolver(A=A, b=b, x0=x0, A0=A0, Ainv0=Ainv0)
496 476
        elif "sym" not in assume_A and "pos" in assume_A:
497 -
            return AsymmetricMatrixBasedSolver(A=A, b=b, x0=x0)
477 +
            raise NotImplementedError
498 478
        else:
499 479
            raise NotImplementedError
500 480

@@ -0,0 +1,32 @@
Loading
1 +
"""Base class for policies of probabilistic linear solvers returning actions."""
2 +
import abc
3 +
4 +
import numpy as np
5 +
6 +
import probnum  # pylint: disable="unused-import"
7 +
8 +
9 +
class LinearSolverPolicy(abc.ABC):
10 +
    r"""Policy of a (probabilistic) linear solver.
11 +
12 +
    The policy :math:`\pi(s \mid \mathsf{A}, \mathsf{H}, \mathsf{x}, A, b)` of a
13 +
    linear solver returns a vector to probe the linear system with, typically via
14 +
    multiplication, resulting in an observation. Policies can either be deterministic or
15 +
    stochastic depending on the application.
16 +
17 +
    See Also
18 +
    --------
19 +
    ConjugateDirectionsPolicy : Policy returning :math:`A`-conjugate actions.
20 +
    """
21 +
22 +
    def __call__(
23 +
        self, solver_state: "probnum.linalg.solvers.ProbabilisticLinearSolverState"
24 +
    ) -> np.ndarray:
25 +
        """Return an action for a given solver state.
26 +
27 +
        Parameters
28 +
        ----------
29 +
        solver_state :
30 +
            Current state of the linear solver.
31 +
        """
32 +
        raise NotImplementedError

@@ -0,0 +1,95 @@
Loading
1 +
"""Generate random linear systems as test problems."""
2 +
3 +
from typing import Any, Callable, Optional, Union
4 +
5 +
import numpy as np
6 +
import scipy.sparse
7 +
8 +
from probnum import linops, problems, randvars
9 +
from probnum.typing import LinearOperatorArgType
10 +
11 +
12 +
def random_linear_system(
13 +
    rng: np.random.Generator,
14 +
    matrix: Union[
15 +
        LinearOperatorArgType,
16 +
        Callable[
17 +
            [np.random.Generator, Optional[Any]],
18 +
            Union[np.ndarray, scipy.sparse.spmatrix, linops.LinearOperator],
19 +
        ],
20 +
    ],
21 +
    solution_rv: Optional[randvars.RandomVariable] = None,
22 +
    **kwargs,
23 +
) -> problems.LinearSystem:
24 +
    """Random linear system.
25 +
26 +
    Generate a random linear system from a (random) matrix. If ``matrix`` is a callable instead of a matrix or
27 +
    linear operator, the system matrix is sampled by passing the random generator instance ``rng``. The solution
28 +
    of the linear system is set to a realization from ``solution_rv``. If ``None`` the solution is drawn from a
29 +
    standard normal distribution with iid components.
30 +
31 +
    Parameters
32 +
    ----------
33 +
    rng
34 +
        Random number generator.
35 +
    matrix
36 +
        Matrix, linear operator or callable returning either for a given random number generator instance.
37 +
    solution_rv
38 +
        Random variable from which the solution of the linear system is sampled.
39 +
    kwargs
40 +
        Miscellaneous arguments passed onto the matrix-generating callable ``matrix``.
41 +
42 +
    See Also
43 +
    --------
44 +
    random_spd_matrix : Generate a random symmetric positive-definite matrix.
45 +
    random_sparse_spd_matrix : Generate a random sparse symmetric positive-definite matrix.
46 +
47 +
    Examples
48 +
    --------
49 +
    >>> import numpy as np
50 +
    >>> from probnum.problems.zoo.linalg import random_linear_system
51 +
    >>> rng = np.random.default_rng(42)
52 +
53 +
    Linear system with given system matrix.
54 +
55 +
    >>> import scipy.stats
56 +
    >>> unitary_matrix = scipy.stats.unitary_group.rvs(dim=5, random_state=rng)
57 +
    >>> linsys_unitary = random_linear_system(rng, unitary_matrix)
58 +
    >>> np.abs(np.linalg.det(linsys_unitary.A))
59 +
    1.0
60 +
61 +
    Linear system with random symmetric positive-definite matrix.
62 +
63 +
    >>> from probnum.problems.zoo.linalg import random_spd_matrix
64 +
    >>> linsys_spd = random_linear_system(rng, random_spd_matrix, dim=2)
65 +
    >>> linsys_spd
66 +
    LinearSystem(A=array([[ 9.62543582,  3.14955953],
67 +
           [ 3.14955953, 13.28720426]]), b=array([-2.7108139 ,  1.10779288]), solution=array([-0.33488503,  0.16275307]))
68 +
69 +
70 +
    Linear system with random sparse matrix.
71 +
72 +
    >>> import scipy.sparse
73 +
    >>> random_sparse_matrix = lambda rng,m,n: scipy.sparse.random(m=m, n=n, random_state=rng)
74 +
    >>> linsys_sparse = random_linear_system(rng, random_sparse_matrix, m=4, n=2)
75 +
    >>> isinstance(linsys_sparse.A, scipy.sparse.spmatrix)
76 +
    True
77 +
    """
78 +
    # Generate system matrix
79 +
    if isinstance(matrix, (np.ndarray, scipy.sparse.spmatrix, linops.LinearOperator)):
80 +
        A = matrix
81 +
    else:
82 +
        A = matrix(rng=rng, **kwargs)
83 +
84 +
    # Sample solution
85 +
    if solution_rv is None:
86 +
        n = A.shape[1]
87 +
        x = randvars.Normal(mean=0.0, cov=1.0).sample(size=(n,), rng=rng)
88 +
    else:
89 +
        if A.shape[1] != solution_rv.shape[0]:
90 +
            raise ValueError(
91 +
                f"Shape of the system matrix: {A.shape} must match shape of the solution: {solution_rv.shape}."
92 +
            )
93 +
        x = solution_rv.sample(size=(), rng=rng)
94 +
95 +
    return problems.LinearSystem(A=A, b=A @ x, solution=x)

@@ -0,0 +1,13 @@
Loading
1 +
"""Policies of probabilistic linear solvers returning actions."""
2 +
3 +
from ._conjugate_gradient import ConjugateGradientPolicy
4 +
from ._linear_solver_policy import LinearSolverPolicy
5 +
from ._random_unit_vector import RandomUnitVectorPolicy
6 +
7 +
# Public classes and functions. Order is reflected in documentation.
8 +
__all__ = ["LinearSolverPolicy", "ConjugateGradientPolicy", "RandomUnitVectorPolicy"]
9 +
10 +
# Set correct module paths. Corrects links and module paths in documentation.
11 +
LinearSolverPolicy.__module__ = "probnum.linalg.solvers.policies"
12 +
RandomUnitVectorPolicy.__module__ = "probnum.linalg.solvers.policies"
13 +
ConjugateGradientPolicy.__module__ = "probnum.linalg.solvers.policies"

@@ -1,5 +1,6 @@
Loading
1 1
"""Test problems from linear algebra."""
2 2
3 +
from ._random_linear_system import random_linear_system
3 4
from ._random_spd_matrix import random_sparse_spd_matrix, random_spd_matrix
4 5
from ._suitesparse_matrix import SuiteSparseMatrix, suitesparse_matrix
5 6
@@ -8,6 +9,7 @@
Loading
8 9
    "random_spd_matrix",
9 10
    "random_sparse_spd_matrix",
10 11
    "suitesparse_matrix",
12 +
    "random_linear_system",
11 13
    "SuiteSparseMatrix",
12 14
]
13 15

@@ -8,6 +8,8 @@
Loading
8 8
9 9
from . import _linear_operator
10 10
11 +
# pylint: disable="too-many-statements"
12 +
11 13
12 14
class Scaling(_linear_operator.LinearOperator):
13 15
    r"""Scaling linear operator.

@@ -0,0 +1,112 @@
Loading
1 +
"""State of a probabilistic linear solver."""
2 +
3 +
import dataclasses
4 +
from typing import Any, List, Optional, Tuple
5 +
6 +
import numpy as np
7 +
8 +
import probnum  # pylint:disable="unused-import"
9 +
from probnum import problems
10 +
11 +
12 +
@dataclasses.dataclass
13 +
class ProbabilisticLinearSolverState:
14 +
    """State of a probabilistic linear solver.
15 +
16 +
    The solver state separates the state of a probabilistic linear solver from the algorithm itself, making the solver stateless. The state contains the problem to be solved, the current belief over the quantities of interest and any miscellaneous quantities computed during an iteration of a probabilistic linear solver. The solver state is passed between the different components of the solver and may be used internally to cache quantities which are used more than once.
17 +
18 +
    Parameters
19 +
    ----------
20 +
    problem
21 +
        Linear system to be solved.
22 +
    prior
23 +
        Prior belief over the quantities of interest of the linear system.
24 +
    rng
25 +
        Random number generator.
26 +
    """
27 +
28 +
    def __init__(
29 +
        self,
30 +
        problem: problems.LinearSystem,
31 +
        prior: "probnum.linalg.solvers.beliefs.LinearSystemBelief",
32 +
        rng: Optional[np.random.Generator] = None,
33 +
    ):
34 +
        self.problem = problem
35 +
        self.belief = prior
36 +
37 +
        self.step = 0
38 +
39 +
        self._actions: List[np.ndarray] = [None]
40 +
        self._observations: List[Any] = [None]
41 +
        self._residuals: List[np.ndarray] = [
42 +
            self.problem.A @ self.belief.x.mean - self.problem.b,
43 +
            None,
44 +
        ]
45 +
        self.rng = rng
46 +
47 +
    def __repr__(self) -> str:
48 +
        return f"{self.__class__.__name__}(step={self.step})"
49 +
50 +
    @property
51 +
    def action(self) -> Optional[np.ndarray]:
52 +
        """Action of the solver for the current step.
53 +
54 +
        Is ``None`` at the beginning of a step and will be set by the
55 +
        policy.
56 +
        """
57 +
        return self._actions[self.step]
58 +
59 +
    @action.setter
60 +
    def action(self, value: np.ndarray) -> None:
61 +
        assert self._actions[self.step] is None
62 +
        self._actions[self.step] = value
63 +
64 +
    @property
65 +
    def observation(self) -> Optional[Any]:
66 +
        """Observation of the solver for the current step.
67 +
68 +
        Is ``None`` at the beginning of a step, will be set by the
69 +
        observation model for a given action.
70 +
        """
71 +
        return self._observations[self.step]
72 +
73 +
    @observation.setter
74 +
    def observation(self, value: Any) -> None:
75 +
        assert self._observations[self.step] is None
76 +
        self._observations[self.step] = value
77 +
78 +
    @property
79 +
    def actions(self) -> Tuple[np.ndarray, ...]:
80 +
        """Actions taken by the solver."""
81 +
        return tuple(self._actions)
82 +
83 +
    @property
84 +
    def observations(self) -> Tuple[Any, ...]:
85 +
        """Observations of the problem by the solver."""
86 +
        return tuple(self._observations)
87 +
88 +
    @property
89 +
    def residual(self) -> np.ndarray:
90 +
        r"""Cached residual :math:`Ax_i-b` for the current solution estimate :math:`x_i`."""
91 +
        if self._residuals[self.step] is None:
92 +
            self._residuals[self.step] = (
93 +
                self.problem.A @ self.belief.x.mean - self.problem.b
94 +
            )
95 +
        return self._residuals[self.step]
96 +
97 +
    @property
98 +
    def residuals(self) -> Tuple[np.ndarray, ...]:
99 +
        r"""Residuals :math:`\{Ax_i - b\}_i`."""
100 +
        return tuple(self._residuals)
101 +
102 +
    def next_step(self) -> None:
103 +
        """Advance the solver state to the next solver step.
104 +
105 +
        Called after a completed step / iteration of the probabilistic
106 +
        linear solver.
107 +
        """
108 +
        self._actions.append(None)
109 +
        self._observations.append(None)
110 +
        self._residuals.append(None)
111 +
112 +
        self.step += 1
Files Coverage
src/probnum 86.40%
Project Totals (124 files) 86.40%
1
coverage:
2
  precision: 2
3
  status:
4
    project:
5
      default:
6
        target: auto
7
        threshold: 1%
8
    patch:
9
      default:
10
        target: 90%
11
        threshold: 1%
12

13
comment:
14
  layout: "reach, diff, files"
15
  behavior: default
16
  require_changes: true
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading