From c24b0ef857030979d3d8a69a7d82f4a2c037d959 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Mon, 26 Jul 2021 18:06:00 +0200 Subject: [PATCH 01/74] Scaling arithmetic --- src/probnum/linops/_arithmetic.py | 5 +++++ src/probnum/linops/_scaling.py | 36 +++++++++++++++++++++++++++++++ 2 files changed, 41 insertions(+) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index 52953dbc5..cdf8291ab 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -57,6 +57,11 @@ def matmul(op1: LinearOperator, op2: LinearOperator) -> LinearOperator: _mul_fns: _BinaryOperatorRegistryType = {} _matmul_fns: _BinaryOperatorRegistryType = {} +_add_fns[(Scaling, Scaling)] = Scaling._add_scaling +_sub_fns[(Scaling, Scaling)] = Scaling._sub_scaling +_mul_fns[(Scaling, Scaling)] = Scaling._mul_scaling +_matmul_fns[(Scaling, Scaling)] = Scaling._matmul_scaling + def _apply( op_registry: _BinaryOperatorRegistryType, diff --git a/src/probnum/linops/_scaling.py b/src/probnum/linops/_scaling.py index 21274b6ab..0fdeaf2be 100644 --- a/src/probnum/linops/_scaling.py +++ b/src/probnum/linops/_scaling.py @@ -201,6 +201,42 @@ def is_isotropic(self) -> bool: """Whether scaling is uniform / isotropic.""" return self._scalar is not None + def _add_scaling(self, other: "Scaling") -> "Scaling": + if other.shape != self.shape: + raise ValueError( + "Addition of two Scaling LinearOperators is only " + "possible if both operands have the same shape." + ) + + return Scaling( + factors=self.factors + other.factors, shape=self.shape, dtype=self.dtype + ) + + def _sub_scaling(self, other: "Scaling") -> "Scaling": + if other.shape != self.shape: + raise ValueError( + "Subtraction of two Scaling LinearOperators is only " + "possible if both operands have the same shape." + ) + + return Scaling( + factors=self.factors - other.factors, shape=self.shape, dtype=self.dtype + ) + + def _mul_scaling(self, other: "Scaling") -> "Scaling": + if other.shape != self.shape: + raise ValueError( + "Multiplication of two Scaling LinearOperators is only " + "possible if both operands have the same shape." + ) + + return Scaling( + factors=self.factors * other.factors, shape=self.shape, dtype=self.dtype + ) + + def _matmul_scaling(self, other: "Scaling") -> "Scaling": + return self._mul_scaling(other) + def _astype(self, dtype, order, casting, copy) -> "Scaling": if self.dtype == dtype and not copy: return self From c18e51a7a3bf9e7e0fc5908ddfea34e86d487a34 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Mon, 26 Jul 2021 18:06:38 +0200 Subject: [PATCH 02/74] Kronecker arithmetic --- src/probnum/linops/_arithmetic.py | 4 ++++ src/probnum/linops/_kronecker.py | 18 ++++++++++++++++++ 2 files changed, 22 insertions(+) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index cdf8291ab..07f50d55a 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -9,6 +9,7 @@ import probnum.utils from probnum.typing import ScalarArgType, ShapeArgType +from ._kronecker import Kronecker from ._linear_operator import ( # pylint: disable=cyclic-import BinaryOperandType, LinearOperator, @@ -62,6 +63,9 @@ def matmul(op1: LinearOperator, op2: LinearOperator) -> LinearOperator: _mul_fns[(Scaling, Scaling)] = Scaling._mul_scaling _matmul_fns[(Scaling, Scaling)] = Scaling._matmul_scaling +_matmul_fns[(Kronecker, Kronecker)] = Kronecker._matmul_kronecker +_mul_fns[(Kronecker, Kronecker)] = Kronecker._mul_kronecker + def _apply( op_registry: _BinaryOperatorRegistryType, diff --git a/src/probnum/linops/_kronecker.py b/src/probnum/linops/_kronecker.py index f50fd4c05..615bd073c 100644 --- a/src/probnum/linops/_kronecker.py +++ b/src/probnum/linops/_kronecker.py @@ -161,6 +161,24 @@ def _cond(self, p) -> np.inexact: return np.linalg.cond(self.todense(cache=False), p=p) + def _matmul_kronecker(self, other: "Kronecker") -> "Kronecker": + _A, _B = self.A, self.B + _C, _D = other.A, other.B + if not (_A.shape[1] == _C.shape[0] and _B.shape[1] == _D.shape[0]): + return NotImplemented + + # Using (A (x) B) @ (C (x) D) = (A @ C) (x) (B @ D) + return Kronecker(A=_A @ _C, B=_B @ _D) + + def _mul_kronecker(self, other: "Kronecker") -> "Kronecker": + _A, _B = self.A, self.B + _C, _D = other.A, other.B + if not (_A.shape == _C.shape and _B.shape == _D.shape): + return NotImplemented + + # Using (A (x) B) o (C (x) D) = (A o C) (x) (B o D) + return Kronecker(A=_A * _C, B=_B * _D) + def _kronecker_matmul( A: _linear_operator.LinearOperator, From 4e5d5b6323cdca0adc20afa6a86712ec6494973c Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Mon, 26 Jul 2021 19:12:21 +0200 Subject: [PATCH 03/74] Identity arithmetic, outsource fallbacks --- src/probnum/linops/_arithmetic.py | 270 +++++--------------- src/probnum/linops/_arithmetic_fallbacks.py | 207 +++++++++++++++ 2 files changed, 273 insertions(+), 204 deletions(-) create mode 100644 src/probnum/linops/_arithmetic_fallbacks.py diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index 07f50d55a..87b293bea 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -1,22 +1,52 @@ """Linear operator arithmetic.""" -import functools -import operator +import numbers from typing import Any, Callable, Dict, Optional, Tuple, Union import numpy as np import scipy.sparse -import probnum.utils -from probnum.typing import ScalarArgType, ShapeArgType +from probnum.typing import ShapeArgType -from ._kronecker import Kronecker -from ._linear_operator import ( # pylint: disable=cyclic-import +from ._arithmetic_fallbacks import ( + NegatedLinearOperator, + ProductLinearOperator, + ScaledLinearOperator, + SumLinearOperator, + _matmul_fallback, + _mul_fallback, +) +from ._kronecker import Kronecker, SymmetricKronecker, Symmetrize +from ._linear_operator import ( + AdjointLinearOperator, BinaryOperandType, + Identity, LinearOperator, Matrix, + TransposedLinearOperator, + _ConjugateLinearOperator, + _InverseLinearOperator, + _TypeCastLinearOperator, ) from ._scaling import Scaling +_AnyLinOp = [ + NegatedLinearOperator, + ProductLinearOperator, + ScaledLinearOperator, + SumLinearOperator, + AdjointLinearOperator, + Identity, + LinearOperator, + Matrix, + TransposedLinearOperator, + SymmetricKronecker, + Symmetrize, + _ConjugateLinearOperator, + _InverseLinearOperator, + _TypeCastLinearOperator, + Scaling, +] + def add(op1: BinaryOperandType, op2: BinaryOperandType) -> LinearOperator: op1, op2 = _operands_to_compatible_linops(op1, op2) @@ -58,14 +88,44 @@ def matmul(op1: LinearOperator, op2: LinearOperator) -> LinearOperator: _mul_fns: _BinaryOperatorRegistryType = {} _matmul_fns: _BinaryOperatorRegistryType = {} +# Scaling _add_fns[(Scaling, Scaling)] = Scaling._add_scaling _sub_fns[(Scaling, Scaling)] = Scaling._sub_scaling _mul_fns[(Scaling, Scaling)] = Scaling._mul_scaling _matmul_fns[(Scaling, Scaling)] = Scaling._matmul_scaling +# Kronecker _matmul_fns[(Kronecker, Kronecker)] = Kronecker._matmul_kronecker _mul_fns[(Kronecker, Kronecker)] = Kronecker._mul_kronecker +# Identity +for op_type in _AnyLinOp: + _matmul_fns[(Identity, op_type)] = lambda idty, other: other + _matmul_fns[(op_type, Identity)] = lambda other, idty: other + + +def _mul_id(arg1, arg2): + if isinstance(arg1, Identity): + if not isinstance(arg2, (int, float, complex, np.number, numbers.Number)): + return NotImplemented + if np.ndim(arg2) != 0: + return NotImplemented + return Scaling(factors=arg2, shape=arg1.shape, dtype=arg1.dtype) + + if isinstance(arg2, Identity): + if not isinstance(arg1, (int, float, complex, np.number, numbers.Number)): + return NotImplemented + if np.ndim(arg1) != 0: + return NotImplemented + return Scaling(factors=arg1, shape=arg2.shape, dtype=arg2.dtype) + + return NotImplemented + + +for sc_type in [int, float, complex, np.number, numbers.Number]: + _mul_fns[(Identity, sc_type)] = _mul_id + _mul_fns[(sc_type, Identity)] = _mul_id + def _apply( op_registry: _BinaryOperatorRegistryType, @@ -123,201 +183,3 @@ def _operands_to_compatible_linops( pass return op1, op2 - - -######################################################################################## -# Generic Linear Operator Arithmetic (Fallbacks) -######################################################################################## - - -class ScaledLinearOperator(LinearOperator): - """Linear operator scaled with a scalar.""" - - def __init__(self, linop: LinearOperator, scalar: ScalarArgType): - if not isinstance(linop, LinearOperator): - raise TypeError("`linop` must be a `LinearOperator`") - - if np.ndim(scalar) != 0: - raise TypeError("`scalar` must be a scalar.") - - dtype = np.result_type(linop.dtype, scalar) - - self._linop = linop - self._scalar = probnum.utils.as_numpy_scalar(scalar, dtype) - - super().__init__( - self._linop.shape, - dtype=dtype, - matmul=lambda x: self._scalar * (self._linop @ x), - rmatmul=lambda x: self._scalar * (x @ self._linop), - todense=lambda: self._scalar * self._linop.todense(cache=False), - transpose=lambda: self._scalar * self._linop.T, - adjoint=lambda: np.conj(self._scalar) * self._linop.H, - inverse=self._inv, - trace=lambda: self._scalar * self._linop.trace(), - ) - - def _inv(self) -> "ScaledLinearOperator": - if self._scalar == 0: - raise np.linalg.LinAlgError("The operator is not invertible") - - return ScaledLinearOperator(self._linop.inv(), 1.0 / self._scalar) - - -class NegatedLinearOperator(ScaledLinearOperator): - def __init__(self, linop: LinearOperator): - super().__init__(linop, scalar=probnum.utils.as_numpy_scalar(-1, linop.dtype)) - - def __neg__(self) -> "LinearOperator": - return self._linop - - -class SumLinearOperator(LinearOperator): - """Sum of two linear operators.""" - - def __init__(self, *summands: LinearOperator): - if not all(isinstance(summand, LinearOperator) for summand in summands): - raise TypeError("All summands must be `LinearOperator`s") - - if len(summands) < 2: - raise ValueError("There must be at least two summands") - - if not all(summand.shape == summands[0].shape for summand in summands): - raise ValueError("All summands must have the same shape") - - self._summands = SumLinearOperator._expand_sum_ops(*summands) - - super().__init__( - shape=summands[0].shape, - dtype=np.find_common_type( - [summand.dtype for summand in self._summands], [] - ), - matmul=lambda x: functools.reduce( - operator.add, (summand @ x for summand in self._summands) - ), - rmatmul=lambda x: functools.reduce( - operator.add, (x @ summand for summand in self._summands) - ), - todense=lambda: functools.reduce( - operator.add, - (summand.todense(cache=False) for summand in self._summands), - ), - transpose=lambda: SumLinearOperator( - *(summand.T for summand in self._summands) - ), - adjoint=lambda: SumLinearOperator( - *(summand.H for summand in self._summands) - ), - trace=lambda: functools.reduce( - operator.add, (summand.trace() for summand in self._summands) - ), - ) - - def __neg__(self): - return SumLinearOperator(*(-summand for summand in self._summands)) - - @staticmethod - def _expand_sum_ops(*summands: LinearOperator) -> Tuple[LinearOperator, ...]: - expanded_summands = [] - - for summand in summands: - if isinstance(summand, SumLinearOperator): - expanded_summands.extend(summand._summands) - else: - expanded_summands.append(summand) - - return tuple(expanded_summands) - - -def _mul_fallback( - op1: BinaryOperandType, op2: BinaryOperandType -) -> Union[LinearOperator, type(NotImplemented)]: - res = NotImplemented - - if isinstance(op1, LinearOperator) and isinstance(op2, LinearOperator): - pass # TODO: Implement generic Hadamard product - elif isinstance(op1, LinearOperator): - if np.ndim(op2) == 0: - res = ScaledLinearOperator(op1, op2) - elif isinstance(op2, LinearOperator): - if np.ndim(op1) == 0: - res = ScaledLinearOperator(op2, op1) - else: - raise TypeError("At least one of the two operands must be a `LinearOperator`.") - - return res - - -class ProductLinearOperator(LinearOperator): - """(Operator) Product of two linear operators.""" - - def __init__(self, *factors: LinearOperator): - if not all(isinstance(factor, LinearOperator) for factor in factors): - raise TypeError("All factors must be `LinearOperator`s") - - if len(factors) < 2: - raise ValueError("There must be at least two factors") - - if not all( - lfactor.shape[1] == rfactor.shape[0] - for lfactor, rfactor in zip(factors[:-1], factors[1:]) - ): - raise ValueError( - f"Shape mismatch: Cannot multiply linear operators with shapes " - f"{' x '.join(factor.shape for factor in factors)}." - ) - - self._factors = ProductLinearOperator._expand_prod_ops(*factors) - - super().__init__( - shape=(self._factors[0].shape[0], self._factors[-1].shape[1]), - dtype=np.find_common_type([factor.dtype for factor in self._factors], []), - matmul=lambda x: functools.reduce( - lambda vec, op: op @ vec, reversed(self._factors), x - ), - rmatmul=lambda x: functools.reduce( - lambda vec, op: vec @ op, self._factors, x - ), - todense=lambda: functools.reduce( - operator.matmul, - (factor.todense(cache=False) for factor in self._factors), - ), - transpose=lambda: ProductLinearOperator( - *(factor.T for factor in reversed(self._factors)) - ), - adjoint=lambda: ProductLinearOperator( - *(factor.H for factor in reversed(self._factors)) - ), - inverse=lambda: ProductLinearOperator( - *(factor.inv() for factor in reversed(self._factors)) - ), - det=lambda: functools.reduce( - operator.mul, (factor.det() for factor in self._factors) - ), - logabsdet=lambda: functools.reduce( - operator.add, (factor.logabsdet() for factor in self._factors) - ), - ) - - @staticmethod - def _expand_prod_ops(*factors: LinearOperator) -> Tuple[LinearOperator, ...]: - expanded_factors = [] - - for factor in factors: - if isinstance(factor, ProductLinearOperator): - expanded_factors.extend(factor._factors) - else: - expanded_factors.append(factor) - - return tuple(expanded_factors) - - -def _matmul_fallback( - op1: BinaryOperandType, op2: BinaryOperandType -) -> Union[LinearOperator, type(NotImplemented)]: - res = NotImplemented - - if isinstance(op1, LinearOperator) and isinstance(op2, LinearOperator): - res = ProductLinearOperator(op1, op2) - - return res diff --git a/src/probnum/linops/_arithmetic_fallbacks.py b/src/probnum/linops/_arithmetic_fallbacks.py new file mode 100644 index 000000000..7833b486c --- /dev/null +++ b/src/probnum/linops/_arithmetic_fallbacks.py @@ -0,0 +1,207 @@ +import functools +import operator +from typing import Tuple, Union + +import numpy as np + +import probnum.utils +from probnum.typing import ScalarArgType + +from ._linear_operator import BinaryOperandType, LinearOperator + +######################################################################################## +# Generic Linear Operator Arithmetic (Fallbacks) +######################################################################################## + + +class ScaledLinearOperator(LinearOperator): + """Linear operator scaled with a scalar.""" + + def __init__(self, linop: LinearOperator, scalar: ScalarArgType): + if not isinstance(linop, LinearOperator): + raise TypeError("`linop` must be a `LinearOperator`") + + if np.ndim(scalar) != 0: + raise TypeError("`scalar` must be a scalar.") + + dtype = np.result_type(linop.dtype, scalar) + + self._linop = linop + self._scalar = probnum.utils.as_numpy_scalar(scalar, dtype) + + super().__init__( + self._linop.shape, + dtype=dtype, + matmul=lambda x: self._scalar * (self._linop @ x), + rmatmul=lambda x: self._scalar * (x @ self._linop), + todense=lambda: self._scalar * self._linop.todense(cache=False), + transpose=lambda: self._scalar * self._linop.T, + adjoint=lambda: np.conj(self._scalar) * self._linop.H, + inverse=self._inv, + trace=lambda: self._scalar * self._linop.trace(), + ) + + def _inv(self) -> "ScaledLinearOperator": + if self._scalar == 0: + raise np.linalg.LinAlgError("The operator is not invertible") + + return ScaledLinearOperator(self._linop.inv(), 1.0 / self._scalar) + + +class NegatedLinearOperator(ScaledLinearOperator): + def __init__(self, linop: LinearOperator): + super().__init__(linop, scalar=probnum.utils.as_numpy_scalar(-1, linop.dtype)) + + def __neg__(self) -> "LinearOperator": + return self._linop + + +class SumLinearOperator(LinearOperator): + """Sum of two linear operators.""" + + def __init__(self, *summands: LinearOperator): + if not all(isinstance(summand, LinearOperator) for summand in summands): + raise TypeError("All summands must be `LinearOperator`s") + + if len(summands) < 2: + raise ValueError("There must be at least two summands") + + if not all(summand.shape == summands[0].shape for summand in summands): + raise ValueError("All summands must have the same shape") + + self._summands = SumLinearOperator._expand_sum_ops(*summands) + + super().__init__( + shape=summands[0].shape, + dtype=np.find_common_type( + [summand.dtype for summand in self._summands], [] + ), + matmul=lambda x: functools.reduce( + operator.add, (summand @ x for summand in self._summands) + ), + rmatmul=lambda x: functools.reduce( + operator.add, (x @ summand for summand in self._summands) + ), + todense=lambda: functools.reduce( + operator.add, + (summand.todense(cache=False) for summand in self._summands), + ), + transpose=lambda: SumLinearOperator( + *(summand.T for summand in self._summands) + ), + adjoint=lambda: SumLinearOperator( + *(summand.H for summand in self._summands) + ), + trace=lambda: functools.reduce( + operator.add, (summand.trace() for summand in self._summands) + ), + ) + + def __neg__(self): + return SumLinearOperator(*(-summand for summand in self._summands)) + + @staticmethod + def _expand_sum_ops(*summands: LinearOperator) -> Tuple[LinearOperator, ...]: + expanded_summands = [] + + for summand in summands: + if isinstance(summand, SumLinearOperator): + expanded_summands.extend(summand._summands) + else: + expanded_summands.append(summand) + + return tuple(expanded_summands) + + +def _mul_fallback( + op1: BinaryOperandType, op2: BinaryOperandType +) -> Union[LinearOperator, type(NotImplemented)]: + res = NotImplemented + + if isinstance(op1, LinearOperator) and isinstance(op2, LinearOperator): + pass # TODO: Implement generic Hadamard product + elif isinstance(op1, LinearOperator): + if np.ndim(op2) == 0: + res = ScaledLinearOperator(op1, op2) + elif isinstance(op2, LinearOperator): + if np.ndim(op1) == 0: + res = ScaledLinearOperator(op2, op1) + else: + raise TypeError("At least one of the two operands must be a `LinearOperator`.") + + return res + + +class ProductLinearOperator(LinearOperator): + """(Operator) Product of two linear operators.""" + + def __init__(self, *factors: LinearOperator): + if not all(isinstance(factor, LinearOperator) for factor in factors): + raise TypeError("All factors must be `LinearOperator`s") + + if len(factors) < 2: + raise ValueError("There must be at least two factors") + + if not all( + lfactor.shape[1] == rfactor.shape[0] + for lfactor, rfactor in zip(factors[:-1], factors[1:]) + ): + raise ValueError( + f"Shape mismatch: Cannot multiply linear operators with shapes " + f"{' x '.join(factor.shape for factor in factors)}." + ) + + self._factors = ProductLinearOperator._expand_prod_ops(*factors) + + super().__init__( + shape=(self._factors[0].shape[0], self._factors[-1].shape[1]), + dtype=np.find_common_type([factor.dtype for factor in self._factors], []), + matmul=lambda x: functools.reduce( + lambda vec, op: op @ vec, reversed(self._factors), x + ), + rmatmul=lambda x: functools.reduce( + lambda vec, op: vec @ op, self._factors, x + ), + todense=lambda: functools.reduce( + operator.matmul, + (factor.todense(cache=False) for factor in self._factors), + ), + transpose=lambda: ProductLinearOperator( + *(factor.T for factor in reversed(self._factors)) + ), + adjoint=lambda: ProductLinearOperator( + *(factor.H for factor in reversed(self._factors)) + ), + inverse=lambda: ProductLinearOperator( + *(factor.inv() for factor in reversed(self._factors)) + ), + det=lambda: functools.reduce( + operator.mul, (factor.det() for factor in self._factors) + ), + logabsdet=lambda: functools.reduce( + operator.add, (factor.logabsdet() for factor in self._factors) + ), + ) + + @staticmethod + def _expand_prod_ops(*factors: LinearOperator) -> Tuple[LinearOperator, ...]: + expanded_factors = [] + + for factor in factors: + if isinstance(factor, ProductLinearOperator): + expanded_factors.extend(factor._factors) + else: + expanded_factors.append(factor) + + return tuple(expanded_factors) + + +def _matmul_fallback( + op1: BinaryOperandType, op2: BinaryOperandType +) -> Union[LinearOperator, type(NotImplemented)]: + res = NotImplemented + + if isinstance(op1, LinearOperator) and isinstance(op2, LinearOperator): + res = ProductLinearOperator(op1, op2) + + return res From 9b6884d01c95fec0eed0c17e899a54a6ef8b110f Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Mon, 26 Jul 2021 19:25:52 +0200 Subject: [PATCH 04/74] Fix pylint --- src/probnum/linops/_arithmetic.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index 87b293bea..a8c330cf5 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -106,18 +106,18 @@ def matmul(op1: LinearOperator, op2: LinearOperator) -> LinearOperator: def _mul_id(arg1, arg2): if isinstance(arg1, Identity): - if not isinstance(arg2, (int, float, complex, np.number, numbers.Number)): - return NotImplemented - if np.ndim(arg2) != 0: - return NotImplemented - return Scaling(factors=arg2, shape=arg1.shape, dtype=arg1.dtype) + if ( + isinstance(arg2, (int, float, complex, np.number, numbers.Number)) + and np.ndim(arg2) == 0 + ): + return Scaling(factors=arg2, shape=arg1.shape, dtype=arg1.dtype) if isinstance(arg2, Identity): - if not isinstance(arg1, (int, float, complex, np.number, numbers.Number)): - return NotImplemented - if np.ndim(arg1) != 0: - return NotImplemented - return Scaling(factors=arg1, shape=arg2.shape, dtype=arg2.dtype) + if ( + isinstance(arg1, (int, float, complex, np.number, numbers.Number)) + and np.ndim(arg1) == 0 + ): + return Scaling(factors=arg1, shape=arg2.shape, dtype=arg2.dtype) return NotImplemented From 92270d510894ded63ef840824aa11909ef94a77b Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Wed, 28 Jul 2021 17:21:49 +0200 Subject: [PATCH 05/74] Kronecker distributive law --- src/probnum/linops/_kronecker.py | 33 ++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/src/probnum/linops/_kronecker.py b/src/probnum/linops/_kronecker.py index 615bd073c..4a28e1e9d 100644 --- a/src/probnum/linops/_kronecker.py +++ b/src/probnum/linops/_kronecker.py @@ -6,6 +6,7 @@ from probnum.typing import DTypeArgType from . import _linear_operator, _utils +from ._scaling import Scaling class Symmetrize(_linear_operator.LinearOperator): @@ -179,6 +180,38 @@ def _mul_kronecker(self, other: "Kronecker") -> "Kronecker": # Using (A (x) B) o (C (x) D) = (A o C) (x) (B o D) return Kronecker(A=_A * _C, B=_B * _D) + def _add_kronecker( + self, other: "Kronecker" + ) -> Union[type(NotImplemented), "Kronecker"]: + + if self.A == other.A and self.B == other.B: + return Kronecker(A=2 * self.A, B=other.B) + + if self.A == other.A: + return Kronecker(A=self.A, B=self.B + other.B) + + if self.B == other.B: + return Kronecker(A=self.A + other.A, B=other.B) + + return NotImplemented + + def _sub_kronecker( + self, other: "Kronecker" + ) -> Union[type(NotImplemented), "Kronecker"]: + + if self.A == other.A and self.B == other.B: + return Kronecker( + A=Scaling(0.0, shape=self.A.shape), B=Scaling(0.0, shape=self.B.shape) + ) + + if self.A == other.A: + return Kronecker(A=self.A, B=self.B - other.B) + + if self.B == other.B: + return Kronecker(A=self.A - other.A, B=other.B) + + return NotImplemented + def _kronecker_matmul( A: _linear_operator.LinearOperator, From d65d7bd39fa97ab9ae4a2f608e513c47f87e7db1 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Wed, 28 Jul 2021 17:21:57 +0200 Subject: [PATCH 06/74] Kronecker shape mismatch raises --- src/probnum/linops/_kronecker.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/probnum/linops/_kronecker.py b/src/probnum/linops/_kronecker.py index 4a28e1e9d..69c5a5413 100644 --- a/src/probnum/linops/_kronecker.py +++ b/src/probnum/linops/_kronecker.py @@ -166,7 +166,10 @@ def _matmul_kronecker(self, other: "Kronecker") -> "Kronecker": _A, _B = self.A, self.B _C, _D = other.A, other.B if not (_A.shape[1] == _C.shape[0] and _B.shape[1] == _D.shape[0]): - return NotImplemented + raise ValueError( + f"Matmul shape mismatch {_A.shape} x {_C.shape} " + f"or {_B.shape} x {_D.shape}" + ) # Using (A (x) B) @ (C (x) D) = (A @ C) (x) (B @ D) return Kronecker(A=_A @ _C, B=_B @ _D) @@ -175,7 +178,10 @@ def _mul_kronecker(self, other: "Kronecker") -> "Kronecker": _A, _B = self.A, self.B _C, _D = other.A, other.B if not (_A.shape == _C.shape and _B.shape == _D.shape): - return NotImplemented + raise ValueError( + f"Matmul shape mismatch {_A.shape} x {_C.shape} " + f"or {_B.shape} x {_D.shape}" + ) # Using (A (x) B) o (C (x) D) = (A o C) (x) (B o D) return Kronecker(A=_A * _C, B=_B * _D) From c9ad2902f7f7a956c614bb3520103902039dfe43 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Wed, 28 Jul 2021 17:22:27 +0200 Subject: [PATCH 07/74] Dense matrix matrix product --- src/probnum/linops/_linear_operator.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/probnum/linops/_linear_operator.py b/src/probnum/linops/_linear_operator.py index cdf47eb5c..5a7c74d4b 100644 --- a/src/probnum/linops/_linear_operator.py +++ b/src/probnum/linops/_linear_operator.py @@ -1012,6 +1012,13 @@ def _sparse_inv(self) -> "Matrix": except RuntimeError as err: raise np.linalg.LinAlgError(str(err)) from err + def _matmul_matrix(self, other: "Matrix") -> "Matrix": + # TODO: Switch via config option + if not self.shape[1] == other.shape[0]: + raise ValueError(f"Matmul shape mismatch {self.shape} x {other.shape}") + + return Matrix(A=self.A @ other.A) + class Identity(LinearOperator): """The identity operator. From 53e7b517a74ffb6efd4a9f718902b6c2eace03b2 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Wed, 28 Jul 2021 17:23:14 +0200 Subject: [PATCH 08/74] Add new ops to arithmetic --- src/probnum/linops/_arithmetic.py | 123 +++++++++++++++++++++++++----- 1 file changed, 102 insertions(+), 21 deletions(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index a8c330cf5..75da77bde 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -94,37 +94,108 @@ def matmul(op1: LinearOperator, op2: LinearOperator) -> LinearOperator: _mul_fns[(Scaling, Scaling)] = Scaling._mul_scaling _matmul_fns[(Scaling, Scaling)] = Scaling._matmul_scaling +for op_type in _AnyLinOp: + _add_fns[(op_type, Scaling)] = ( + lambda op, scaling: op if scaling._scalar == 0.0 else NotImplemented + ) + _add_fns[(Scaling, op_type)] = ( + lambda scaling, op: op if scaling._scalar == 0.0 else NotImplemented + ) + _sub_fns[(Scaling, op_type)] = ( + lambda scaling, op: op if scaling._scalar == 0.0 else NotImplemented + ) + _sub_fns[(op_type, Scaling)] = ( + lambda op, scaling: op if scaling._scalar == 0.0 else NotImplemented + ) + +# ScaledLinearOperator +def _matmul_scaled_op(scaled, anylinop): + return scaled._scalar * (scaled._linop @ anylinop) + + +def _matmul_op_scaled(anylinop, scaled): + return scaled._scalar * (anylinop @ scaled._linop) + + +for op_type in _AnyLinOp: + _matmul_fns[(ScaledLinearOperator, op_type)] = _matmul_scaled_op + _matmul_fns[(op_type, ScaledLinearOperator)] = _matmul_op_scaled + + # Kronecker + +_add_fns[(Kronecker, Kronecker)] = Kronecker._add_kronecker +_sub_fns[(Kronecker, Kronecker)] = Kronecker._sub_kronecker + + +def _matmul_scaling_kronecker(scaling: Scaling, kronecker: Kronecker) -> Kronecker: + if scaling.is_isotropic: + return Kronecker(A=scaling.scalar * kronecker.A, B=kronecker.B) + return NotImplemented + + +def _matmul_kronecker_scaling(kronecker: Kronecker, scaling: Scaling) -> Kronecker: + if scaling.is_isotropic: + return Kronecker(A=scaling.scalar * kronecker.A, B=kronecker.B) + return NotImplemented + + _matmul_fns[(Kronecker, Kronecker)] = Kronecker._matmul_kronecker _mul_fns[(Kronecker, Kronecker)] = Kronecker._mul_kronecker -# Identity -for op_type in _AnyLinOp: - _matmul_fns[(Identity, op_type)] = lambda idty, other: other - _matmul_fns[(op_type, Identity)] = lambda other, idty: other +_mul_fns[("scalar", Kronecker)] = lambda sc, kr: Kronecker(A=sc * kr.A, B=kr.B) +_mul_fns[(Kronecker, "scalar")] = lambda kr, sc: Kronecker(A=sc * kr.A, B=kr.B) +_matmul_fns[(Kronecker, Scaling)] = _matmul_kronecker_scaling +_matmul_fns[(Scaling, Kronecker)] = _matmul_scaling_kronecker + +# Matrix +def _matmul_scaling_matrix(scaling: Scaling, matrix: Matrix) -> Matrix: + if scaling.shape[1] != matrix.shape[0]: + return NotImplemented + + return Matrix(A=np.multiply(scaling.factors[:, np.newaxis], matrix.A)) + + +def _matmul_matrix_scaling(matrix: Matrix, scaling: Scaling) -> Matrix: + if matrix.shape[1] != scaling.shape[0]: + return NotImplemented + return Matrix(A=np.multiply(matrix.A, scaling.factors)) -def _mul_id(arg1, arg2): - if isinstance(arg1, Identity): - if ( - isinstance(arg2, (int, float, complex, np.number, numbers.Number)) - and np.ndim(arg2) == 0 - ): - return Scaling(factors=arg2, shape=arg1.shape, dtype=arg1.dtype) - if isinstance(arg2, Identity): - if ( - isinstance(arg1, (int, float, complex, np.number, numbers.Number)) - and np.ndim(arg1) == 0 - ): - return Scaling(factors=arg1, shape=arg2.shape, dtype=arg2.dtype) +def _mul_matrix_scalar(mat, scalar) -> Union[type(NotImplemented), Matrix]: + if np.isscalar(scalar): + return Matrix(A=scalar * mat.A) return NotImplemented -for sc_type in [int, float, complex, np.number, numbers.Number]: - _mul_fns[(Identity, sc_type)] = _mul_id - _mul_fns[(sc_type, Identity)] = _mul_id +def _mul_scalar_matrix(scalar, mat) -> Union[type(NotImplemented), Matrix]: + if np.isscalar(scalar): + return Matrix(A=scalar * mat.A) + + return NotImplemented + + +_mul_fns[(Matrix, "scalar")] = _mul_matrix_scalar +_mul_fns[("scalar", Matrix)] = _mul_scalar_matrix + +_matmul_fns[(Matrix, Matrix)] = Matrix._matmul_matrix +_matmul_fns[(Scaling, Matrix)] = _matmul_scaling_matrix +_matmul_fns[(Matrix, Scaling)] = _matmul_matrix_scaling + + +# Identity +for op_type in _AnyLinOp: + _matmul_fns[(Identity, op_type)] = lambda idty, other: other + _matmul_fns[(op_type, Identity)] = lambda other, idty: other + +_mul_fns[(Identity, "scalar")] = lambda idty, sc: Scaling( + sc, shape=idty.shape, dtype=idty.dtype +) +_mul_fns[("scalar", Identity)] = lambda sc, idty: Scaling( + sc, shape=idty.shape, dtype=idty.dtype +) def _apply( @@ -138,7 +209,17 @@ def _apply( ] ] = None, ) -> Union[LinearOperator, type(NotImplemented)]: - key = (type(op1), type(op2)) + if np.isscalar(op1): + key1 = "scalar" + else: + key1 = type(op1) + + if np.isscalar(op2): + key2 = "scalar" + else: + key2 = type(op2) + + key = (key1, key2) if key in op_registry: res = op_registry[key](op1, op2) From 0e098864717ae1d84ac36fe8a25fc08c75bc23b8 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Wed, 28 Jul 2021 18:20:23 +0200 Subject: [PATCH 09/74] Scalar Identity arithmetic --- src/probnum/linops/_arithmetic.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index 75da77bde..b3cb8fa0c 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -190,12 +190,8 @@ def _mul_scalar_matrix(scalar, mat) -> Union[type(NotImplemented), Matrix]: _matmul_fns[(Identity, op_type)] = lambda idty, other: other _matmul_fns[(op_type, Identity)] = lambda other, idty: other -_mul_fns[(Identity, "scalar")] = lambda idty, sc: Scaling( - sc, shape=idty.shape, dtype=idty.dtype -) -_mul_fns[("scalar", Identity)] = lambda sc, idty: Scaling( - sc, shape=idty.shape, dtype=idty.dtype -) +_mul_fns[(Identity, "scalar")] = lambda idty, sc: Scaling(sc, shape=idty.shape) +_mul_fns[("scalar", Identity)] = lambda sc, idty: Scaling(sc, shape=idty.shape) def _apply( From 813c7f7ebbcd0e1bf7a7a311e6e036b8ed7b6ad5 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Wed, 28 Jul 2021 18:20:35 +0200 Subject: [PATCH 10/74] Scalar scaling arithmetic --- src/probnum/linops/_arithmetic.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index b3cb8fa0c..a86147bbf 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -108,6 +108,24 @@ def matmul(op1: LinearOperator, op2: LinearOperator) -> LinearOperator: lambda op, scaling: op if scaling._scalar == 0.0 else NotImplemented ) + +def _mul_scalar_scaling(scalar: ScalarArgType, scaling: Scaling) -> Scaling: + if scaling.is_isotropic: + return Scaling(scalar * scaling.scalar, shape=scaling.shape) + + return Scaling(scalar * scaling.factors, shape=scaling.shape) + + +def _mul_scaling_scalar(scaling: Scaling, scalar: ScalarArgType) -> Scaling: + if scaling.is_isotropic: + return Scaling(scalar * scaling.scalar, shape=scaling.shape) + + return Scaling(scalar * scaling.factors, shape=scaling.shape) + + +_mul_fns[("scalar", Scaling)] = _mul_scalar_scaling +_mul_fns[(Scaling, "scalar")] = _mul_scaling_scalar + # ScaledLinearOperator def _matmul_scaled_op(scaled, anylinop): return scaled._scalar * (scaled._linop @ anylinop) From 4b2ac0fe174a4673ef258e8c8f7283b5ac5515d5 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Wed, 28 Jul 2021 18:20:53 +0200 Subject: [PATCH 11/74] rmult for scalar Kronecker --- src/probnum/linops/_arithmetic.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index a86147bbf..fe2c77788 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -5,7 +5,7 @@ import numpy as np import scipy.sparse -from probnum.typing import ShapeArgType +from probnum.typing import ScalarArgType, ShapeArgType from ._arithmetic_fallbacks import ( NegatedLinearOperator, @@ -154,7 +154,7 @@ def _matmul_scaling_kronecker(scaling: Scaling, kronecker: Kronecker) -> Kroneck def _matmul_kronecker_scaling(kronecker: Kronecker, scaling: Scaling) -> Kronecker: if scaling.is_isotropic: - return Kronecker(A=scaling.scalar * kronecker.A, B=kronecker.B) + return Kronecker(A=kronecker.A, B=kronecker.B * scaling.scalar) return NotImplemented @@ -162,7 +162,7 @@ def _matmul_kronecker_scaling(kronecker: Kronecker, scaling: Scaling) -> Kroneck _mul_fns[(Kronecker, Kronecker)] = Kronecker._mul_kronecker _mul_fns[("scalar", Kronecker)] = lambda sc, kr: Kronecker(A=sc * kr.A, B=kr.B) -_mul_fns[(Kronecker, "scalar")] = lambda kr, sc: Kronecker(A=sc * kr.A, B=kr.B) +_mul_fns[(Kronecker, "scalar")] = lambda kr, sc: Kronecker(A=kr.A, B=kr.B * sc) _matmul_fns[(Kronecker, Scaling)] = _matmul_kronecker_scaling _matmul_fns[(Scaling, Kronecker)] = _matmul_scaling_kronecker From f59d37d151af07a205c6b29dbc8c35e6f124fd11 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Wed, 28 Jul 2021 18:21:14 +0200 Subject: [PATCH 12/74] Equality for Scaling --- src/probnum/linops/_scaling.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/probnum/linops/_scaling.py b/src/probnum/linops/_scaling.py index 0fdeaf2be..6df61b03b 100644 --- a/src/probnum/linops/_scaling.py +++ b/src/probnum/linops/_scaling.py @@ -201,6 +201,18 @@ def is_isotropic(self) -> bool: """Whether scaling is uniform / isotropic.""" return self._scalar is not None + def __eq__(self, other: _linear_operator.LinearOperator) -> bool: + if self._is_type_shape_dtype_equal(other): + return False + + if self.is_isotropic and other.is_isotropic: + return self.scalar == other.scalar + + if not self.is_isotropic and not self.is_isotropic: + return np.all(self.factors == other.factors) + + return False + def _add_scaling(self, other: "Scaling") -> "Scaling": if other.shape != self.shape: raise ValueError( From 1190b760426350b70f07e351040df3ba82409452 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Wed, 28 Jul 2021 18:22:08 +0200 Subject: [PATCH 13/74] Equality --- src/probnum/linops/_linear_operator.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/probnum/linops/_linear_operator.py b/src/probnum/linops/_linear_operator.py index 5a7c74d4b..3bbaec0ec 100644 --- a/src/probnum/linops/_linear_operator.py +++ b/src/probnum/linops/_linear_operator.py @@ -563,6 +563,16 @@ def __rmul__(self, other: BinaryOperandType) -> "LinearOperator": return mul(other, self) + def _is_type_shape_dtype_equal(self, other: "LinearOperator") -> bool: + return ( + type(self) == type(other) + and self.shape == other.shape + and self.dtype == other.dtype + ) + + def __eq__(self, other: "LinearOperator") -> bool: + raise NotImplementedError + def __matmul__( self, other: BinaryOperandType ) -> Union["LinearOperator", np.ndarray]: @@ -1019,6 +1029,12 @@ def _matmul_matrix(self, other: "Matrix") -> "Matrix": return Matrix(A=self.A @ other.A) + def __eq__(self, other: LinearOperator) -> bool: + if self._is_type_shape_dtype_equal(other): + return False + + return np.all(self.A == other.A) + class Identity(LinearOperator): """The identity operator. @@ -1086,3 +1102,6 @@ def _astype( return self else: return Identity(self.shape, dtype=dtype) + + def __eq__(self, other: LinearOperator) -> bool: + return self._is_type_shape_dtype_equal(other) From 9861def9a80bd1f8e562e25527d95ee416e16df5 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 29 Jul 2021 12:01:33 +0200 Subject: [PATCH 14/74] Fix pylint issues --- src/probnum/linops/_arithmetic.py | 1 - src/probnum/linops/_linear_operator.py | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index fe2c77788..0f49cb865 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -1,5 +1,4 @@ """Linear operator arithmetic.""" -import numbers from typing import Any, Callable, Dict, Optional, Tuple, Union import numpy as np diff --git a/src/probnum/linops/_linear_operator.py b/src/probnum/linops/_linear_operator.py index 3bbaec0ec..669ca8fa2 100644 --- a/src/probnum/linops/_linear_operator.py +++ b/src/probnum/linops/_linear_operator.py @@ -565,7 +565,7 @@ def __rmul__(self, other: BinaryOperandType) -> "LinearOperator": def _is_type_shape_dtype_equal(self, other: "LinearOperator") -> bool: return ( - type(self) == type(other) + isinstance(self, type(other)) and self.shape == other.shape and self.dtype == other.dtype ) From c77c5f03ba2880b5c69d1b0bedb842c3f9547878 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 29 Jul 2021 18:16:11 +0200 Subject: [PATCH 15/74] Zero LinOp --- src/probnum/linops/__init__.py | 3 +-- src/probnum/linops/_scaling.py | 31 ++++++++++++++++++++++++++++++- 2 files changed, 31 insertions(+), 3 deletions(-) diff --git a/src/probnum/linops/__init__.py b/src/probnum/linops/__init__.py index bb79ceec0..ed707ab45 100644 --- a/src/probnum/linops/__init__.py +++ b/src/probnum/linops/__init__.py @@ -13,8 +13,7 @@ """ from ._kronecker import Kronecker, SymmetricKronecker, Symmetrize -from ._linear_operator import Identity, LinearOperator, Matrix -from ._scaling import Scaling +from ._scaling import Scaling, Zero from ._utils import LinearOperatorLike, aslinop # Public classes and functions. Order is reflected in documentation. diff --git a/src/probnum/linops/_scaling.py b/src/probnum/linops/_scaling.py index 6df61b03b..223dbf9bd 100644 --- a/src/probnum/linops/_scaling.py +++ b/src/probnum/linops/_scaling.py @@ -202,7 +202,7 @@ def is_isotropic(self) -> bool: return self._scalar is not None def __eq__(self, other: _linear_operator.LinearOperator) -> bool: - if self._is_type_shape_dtype_equal(other): + if not self._is_type_shape_dtype_equal(other): return False if self.is_isotropic and other.is_isotropic: @@ -317,3 +317,32 @@ def _cond_isotropic(self, p: Union[None, int, float, str]) -> np.inexact: ) else: return np.linalg.cond(self.todense(cache=False), p=p) + + +class Zero(_linear_operator.LinearOperator): + def __init__(self, shape, dtype=np.float64): + + zero = np.zeros(shape=(), dtype=dtype) + matmul = lambda x: zero * x + rmatmul = lambda x: zero * x + apply = lambda x, axis: zero * x + todense = lambda: np.zeros(shape=shape, dtype=dtype) + rank = lambda: np.intp(0) + eigvals = lambda: np.full(shape[0], zero, dtype=dtype) + det = lambda: (zero.astype(dtype, copy=False) ** shape[0]) + + trace = lambda: zero + + super().__init__( + shape, + dtype=dtype, + matmul=matmul, + rmatmul=rmatmul, + apply=apply, + todense=todense, + transpose=lambda: self, + rank=rank, + eigvals=eigvals, + det=det, + trace=trace, + ) From ec23359b8b9574067998721637d1e943db93dae5 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 29 Jul 2021 18:17:19 +0200 Subject: [PATCH 16/74] Selection LinOp --- src/probnum/linops/__init__.py | 1 + src/probnum/linops/_linear_operator.py | 46 +++++++++++++++++++++++++- 2 files changed, 46 insertions(+), 1 deletion(-) diff --git a/src/probnum/linops/__init__.py b/src/probnum/linops/__init__.py index ed707ab45..d14c1bf69 100644 --- a/src/probnum/linops/__init__.py +++ b/src/probnum/linops/__init__.py @@ -13,6 +13,7 @@ """ from ._kronecker import Kronecker, SymmetricKronecker, Symmetrize +from ._linear_operator import Identity, LinearOperator, Matrix, Selection from ._scaling import Scaling, Zero from ._utils import LinearOperatorLike, aslinop diff --git a/src/probnum/linops/_linear_operator.py b/src/probnum/linops/_linear_operator.py index 669ca8fa2..f48fb4db1 100644 --- a/src/probnum/linops/_linear_operator.py +++ b/src/probnum/linops/_linear_operator.py @@ -1030,7 +1030,7 @@ def _matmul_matrix(self, other: "Matrix") -> "Matrix": return Matrix(A=self.A @ other.A) def __eq__(self, other: LinearOperator) -> bool: - if self._is_type_shape_dtype_equal(other): + if not self._is_type_shape_dtype_equal(other): return False return np.all(self.A == other.A) @@ -1105,3 +1105,47 @@ def _astype( def __eq__(self, other: LinearOperator) -> bool: return self._is_type_shape_dtype_equal(other) + + +class Selection(LinearOperator): + def __init__(self, indices, shape): + self._indices = probnum.utils.as_shape(indices) # TODO + assert len(self._indices) == shape[0] or len(self._indices) == shape[1] + + super().__init__( + dtype=np.uint8, # TODO + shape=shape, + matmul=lambda x: _selection_matmul(self.indices, x), + rmatmul=lambda x: _selection_rmatmul(self.indices, x), + todense=lambda: self._todense(), + transpose=lambda: Selection( + indices, shape=(self.shape[1], self.shape[0]) + ), # TODO + ) + + @property + def indices(self): + return self._indices + + def _todense(self): + res = np.eye(self.shape[1], self.shape[1]) + return _selection_matmul(self, res) + + def _matmul_selection( + self, other: "Selection" + ) -> Union[type(NotImplemented), "Identity"]: + + if (self.shape[-1] == other.shape[-2]) and np.all( + self.indices == other.indices + ): + return Identity(shape=(self.shape[-2], other.shape[-1])) + + return NotImplemented + + +def _selection_matmul(indices, M): + return np.take(M, indices=indices, axis=-2) + + +def _selection_rmatmul(indices, M): + return np.take(M, indices=indices, axis=-1) From 52db4b34573c83dab0d60bc0e8ae54db958b1f84 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 29 Jul 2021 18:17:36 +0200 Subject: [PATCH 17/74] Zero and Selection arithmetic --- src/probnum/linops/_arithmetic.py | 56 +++++++++++++++++++++++++++++-- 1 file changed, 53 insertions(+), 3 deletions(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index 0f49cb865..ec5305b46 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -21,12 +21,13 @@ Identity, LinearOperator, Matrix, + Selection, TransposedLinearOperator, _ConjugateLinearOperator, _InverseLinearOperator, _TypeCastLinearOperator, ) -from ._scaling import Scaling +from ._scaling import Scaling, Zero _AnyLinOp = [ NegatedLinearOperator, @@ -44,6 +45,9 @@ _InverseLinearOperator, _TypeCastLinearOperator, Scaling, + Selection, + Zero, + Kronecker, ] @@ -147,13 +151,19 @@ def _matmul_op_scaled(anylinop, scaled): def _matmul_scaling_kronecker(scaling: Scaling, kronecker: Kronecker) -> Kronecker: if scaling.is_isotropic: - return Kronecker(A=scaling.scalar * kronecker.A, B=kronecker.B) + if ("scalar", type(kronecker.A)) in _mul_fns: + return Kronecker(A=scaling.scalar * kronecker.A, B=kronecker.B) + if ("scalar", type(kronecker.B)) in _mul_fns: + return Kronecker(A=kronecker.A, B=scaling.scalar * kronecker.B) return NotImplemented def _matmul_kronecker_scaling(kronecker: Kronecker, scaling: Scaling) -> Kronecker: if scaling.is_isotropic: - return Kronecker(A=kronecker.A, B=kronecker.B * scaling.scalar) + if (type(kronecker.B), "scalar") in _mul_fns: + return Kronecker(A=kronecker.A, B=kronecker.B * scaling.scalar) + if (type(kronecker.A), "scalar") in _mul_fns: + return Kronecker(A=kronecker.A, B=kronecker.B * scaling.scalar) return NotImplemented @@ -211,6 +221,46 @@ def _mul_scalar_matrix(scalar, mat) -> Union[type(NotImplemented), Matrix]: _mul_fns[("scalar", Identity)] = lambda sc, idty: Scaling(sc, shape=idty.shape) +# Selection +_matmul_fns[(Selection, Selection)] = Selection._matmul_selection + +# Zero +def _matmul_zero_anylinop(z: Zero, op: LinearOperator) -> Zero: + if z.shape[1] != op.shape[0]: + raise ValueError(f"shape mismatch") # TODO + + return Zero(shape=(z.shape[0], op.shape[1])) + + +def _matmul_anylinop_zero(op: LinearOperator, z: Zero) -> Zero: + if z.shape[0] != op.shape[1]: + raise ValueError(f"shape mismatch") # TODO + + return Zero(shape=(op.shape[0], z.shape[1])) + + +def _add_zero_anylinop(z: Zero, op: LinearOperator) -> Zero: + if z.shape != op.shape: + raise ValueError(f"shape mismatch") # TODO + + return op + + +def _add_anylinop_zero(op: LinearOperator, z: Zero) -> Zero: + if z.shape != op.shape: + raise ValueError(f"shape mismatch") # TODO + + return op + + +for op_type in _AnyLinOp: + _matmul_fns[(Zero, op_type)] = _matmul_zero_anylinop + _matmul_fns[(op_type, Zero)] = _matmul_anylinop_zero + _add_fns[(Zero, op_type)] = _add_zero_anylinop + _add_fns[(op_type, Zero)] = _add_anylinop_zero + # TODO scalar mult + + def _apply( op_registry: _BinaryOperatorRegistryType, op1: LinearOperator, From 2516f94691d3c9d06389c859e5f9ff5e56a5acc2 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 29 Jul 2021 18:18:04 +0200 Subject: [PATCH 18/74] Allow kronecker cov for vector-valued normal rvs --- src/probnum/randvars/_normal.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/probnum/randvars/_normal.py b/src/probnum/randvars/_normal.py index 29b769636..5daa700a7 100644 --- a/src/probnum/randvars/_normal.py +++ b/src/probnum/randvars/_normal.py @@ -193,20 +193,20 @@ def __init__( self._symmetric_kronecker_identical_factors_cov_cholesky ) elif isinstance(cov, linops.Kronecker): - m, n = mean.shape - - if ( - m != cov.A.shape[0] - or m != cov.A.shape[1] - or n != cov.B.shape[0] - or n != cov.B.shape[1] - ): - raise ValueError( - "Kronecker structured kernels must have factors with the same " - "shape as the mean." - ) - compute_cov_cholesky = self._kronecker_cov_cholesky + if mean.ndim == 2: + m, n = mean.shape + + if ( + m != cov.A.shape[0] + or m != cov.A.shape[1] + or n != cov.B.shape[0] + or n != cov.B.shape[1] + ): + raise ValueError( + "Kronecker structured kernels must have factors with the same " + "shape as the mean." + ) else: # This case handles all linear operators, for which no Cholesky From 3cab34d7454c6094a3dc2846acd5ae0683277692 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 29 Jul 2021 18:18:14 +0200 Subject: [PATCH 19/74] Use zero linop in constant --- src/probnum/randvars/_constant.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/probnum/randvars/_constant.py b/src/probnum/randvars/_constant.py index 5d9b775f2..1f24a0008 100644 --- a/src/probnum/randvars/_constant.py +++ b/src/probnum/randvars/_constant.py @@ -75,11 +75,7 @@ def __init__( if config.matrix_free: cov = lambda: ( - linops.Scaling( - 0.0, - shape=(self._support.size, self._support.size), - dtype=support_floating.dtype, - ) + linops.Zero(shape=((self._support.size, self._support.size))) if self._support.ndim > 0 else _utils.as_numpy_scalar(0.0, support_floating.dtype) ) From f58b9512c91907de948e315988a2d060950db3cc Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 29 Jul 2021 19:01:49 +0200 Subject: [PATCH 20/74] Remove stupid case from kronecker --- src/probnum/linops/_kronecker.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/probnum/linops/_kronecker.py b/src/probnum/linops/_kronecker.py index 69c5a5413..6a1dfa149 100644 --- a/src/probnum/linops/_kronecker.py +++ b/src/probnum/linops/_kronecker.py @@ -190,9 +190,6 @@ def _add_kronecker( self, other: "Kronecker" ) -> Union[type(NotImplemented), "Kronecker"]: - if self.A == other.A and self.B == other.B: - return Kronecker(A=2 * self.A, B=other.B) - if self.A == other.A: return Kronecker(A=self.A, B=self.B + other.B) From 9fc34f702ad5b1ee0c475d148e2f83de345760b3 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 29 Jul 2021 19:02:00 +0200 Subject: [PATCH 21/74] Fix arithmetic --- src/probnum/linops/_arithmetic.py | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index ec5305b46..aeeec8181 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -97,20 +97,6 @@ def matmul(op1: LinearOperator, op2: LinearOperator) -> LinearOperator: _mul_fns[(Scaling, Scaling)] = Scaling._mul_scaling _matmul_fns[(Scaling, Scaling)] = Scaling._matmul_scaling -for op_type in _AnyLinOp: - _add_fns[(op_type, Scaling)] = ( - lambda op, scaling: op if scaling._scalar == 0.0 else NotImplemented - ) - _add_fns[(Scaling, op_type)] = ( - lambda scaling, op: op if scaling._scalar == 0.0 else NotImplemented - ) - _sub_fns[(Scaling, op_type)] = ( - lambda scaling, op: op if scaling._scalar == 0.0 else NotImplemented - ) - _sub_fns[(op_type, Scaling)] = ( - lambda op, scaling: op if scaling._scalar == 0.0 else NotImplemented - ) - def _mul_scalar_scaling(scalar: ScalarArgType, scaling: Scaling) -> Scaling: if scaling.is_isotropic: From 350c58d11f3f9c72b30cf55268c0d9212828f4b4 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Fri, 30 Jul 2021 12:12:11 +0200 Subject: [PATCH 22/74] Selection Embedding --- src/probnum/linops/_linear_operator.py | 76 +++++++++++++++++++------- 1 file changed, 56 insertions(+), 20 deletions(-) diff --git a/src/probnum/linops/_linear_operator.py b/src/probnum/linops/_linear_operator.py index f48fb4db1..b0d470b7e 100644 --- a/src/probnum/linops/_linear_operator.py +++ b/src/probnum/linops/_linear_operator.py @@ -7,7 +7,7 @@ import scipy.sparse.linalg import probnum.utils -from probnum.typing import DTypeArgType, ScalarArgType, ShapeArgType +from probnum.typing import DTypeArgType, NotImplementedType, ScalarArgType, ShapeArgType BinaryOperandType = Union[ "LinearOperator", ScalarArgType, np.ndarray, scipy.sparse.spmatrix @@ -785,6 +785,9 @@ def _astype( ) -> "LinearOperator": return self._linop.astype(dtype, order=order, casting=casting, copy=copy).T + def __repr__(self) -> str: + return f"Transpose of {self._linop}" + class AdjointLinearOperator(LinearOperator): def __init__( @@ -870,6 +873,9 @@ def __init__(self, linop: LinearOperator): logabsdet=lambda: -self._linop.logabsdet(), ) + def __repr__(self) -> str: + return f"Inverse of {self._linop}" + @property def factorization(self): if self.__factorization is None: @@ -1109,18 +1115,22 @@ def __eq__(self, other: LinearOperator) -> bool: class Selection(LinearOperator): def __init__(self, indices, shape): + if np.ndim(indices) > 1: + raise ValueError(f"bla") # TODO self._indices = probnum.utils.as_shape(indices) # TODO - assert len(self._indices) == shape[0] or len(self._indices) == shape[1] + assert shape[0] <= shape[1] + assert len(self._indices) == shape[0] super().__init__( dtype=np.uint8, # TODO shape=shape, matmul=lambda x: _selection_matmul(self.indices, x), - rmatmul=lambda x: _selection_rmatmul(self.indices, x), todense=lambda: self._todense(), - transpose=lambda: Selection( - indices, shape=(self.shape[1], self.shape[0]) - ), # TODO + transpose=lambda: Embedding( + take_indices=np.arange(len(self._indices)), + put_indices=self._indices, + shape=(self.shape[1], self.shape[0]), + ), ) @property @@ -1129,23 +1139,49 @@ def indices(self): def _todense(self): res = np.eye(self.shape[1], self.shape[1]) - return _selection_matmul(self, res) - - def _matmul_selection( - self, other: "Selection" - ) -> Union[type(NotImplemented), "Identity"]: - - if (self.shape[-1] == other.shape[-2]) and np.all( - self.indices == other.indices - ): - return Identity(shape=(self.shape[-2], other.shape[-1])) - - return NotImplemented + return _selection_matmul(self.indices, res) def _selection_matmul(indices, M): return np.take(M, indices=indices, axis=-2) -def _selection_rmatmul(indices, M): - return np.take(M, indices=indices, axis=-1) +class Embedding(LinearOperator): + def __init__(self, take_indices, put_indices, shape, fill_value=0.0): + self._take_indices = probnum.utils.as_shape(take_indices) # TODO + self._put_indices = probnum.utils.as_shape(put_indices) # TODO + self._fill_value = fill_value # TODO + + super().__init__( + dtype=np.uint8, # TODO + shape=shape, + matmul=lambda x: _embedding_matmul(self, x), + todense=lambda: self._todense(), + transpose=lambda: Selection( + indices=put_indices, shape=(self.shape[1], self.shape[0]) + ), + ) + + def _todense(self): + res = np.eye(self.shape[1], self.shape[1]) + return _embedding_matmul(self, res) + + +def _embedding_matmul(embedding, M): + res_shape = np.array(M.shape) + res_shape[-2] = embedding.shape[0] + res = np.full(shape=tuple(res_shape), fill_value=embedding._fill_value) + # print(f"Take indices {embedding._take_indices}") + # print(f"Put indices {embedding._put_indices}") + # print(f"Res, {res}") + # res[np.array(embedding._put_indices).reshape(-1, 1)] = M[embedding._take_indices] + # TODO + # take_from_M = np.take_along_axis( + # M, indices=np.array(embedding._put_indices), axis=-2 + # ) + take_from_M = M[..., np.array(embedding._put_indices), :] + # np.put_along_axis( + # res, indices=np.array(embedding._put_indices), values=take_from_M, axis=-2 + # ) + res[..., np.array(embedding._put_indices), :] = take_from_M + return res From 3a605bfa62a5c2c14a87550fbc05fb06f675236e Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Fri, 30 Jul 2021 12:12:20 +0200 Subject: [PATCH 23/74] Add NotImplementedType --- src/probnum/typing.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/probnum/typing.py b/src/probnum/typing.py index 51e5362f6..c7592f58a 100644 --- a/src/probnum/typing.py +++ b/src/probnum/typing.py @@ -81,3 +81,5 @@ DenseOutputLocationArgType = Union[FloatArgType, np.ndarray] """TimeSeriesPosteriors and derived classes can be evaluated at a single location 't' or an array of locations.""" + +NotImplementedType = type(NotImplemented) From 4f9ef1521c3e73fadfb32fc112a2ce65754f3005 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Fri, 30 Jul 2021 12:12:46 +0200 Subject: [PATCH 24/74] Advances in arithmetics --- src/probnum/linops/__init__.py | 2 +- src/probnum/linops/_arithmetic.py | 128 +++++++++++++++----- src/probnum/linops/_arithmetic_fallbacks.py | 30 ++++- src/probnum/linops/_kronecker.py | 6 +- 4 files changed, 128 insertions(+), 38 deletions(-) diff --git a/src/probnum/linops/__init__.py b/src/probnum/linops/__init__.py index d14c1bf69..f12ae2fce 100644 --- a/src/probnum/linops/__init__.py +++ b/src/probnum/linops/__init__.py @@ -13,7 +13,7 @@ """ from ._kronecker import Kronecker, SymmetricKronecker, Symmetrize -from ._linear_operator import Identity, LinearOperator, Matrix, Selection +from ._linear_operator import Embedding, Identity, LinearOperator, Matrix, Selection from ._scaling import Scaling, Zero from ._utils import LinearOperatorLike, aslinop diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index aeeec8181..55a7765ce 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -4,7 +4,7 @@ import numpy as np import scipy.sparse -from probnum.typing import ScalarArgType, ShapeArgType +from probnum.typing import NotImplementedType, ScalarArgType, ShapeArgType from ._arithmetic_fallbacks import ( NegatedLinearOperator, @@ -18,6 +18,7 @@ from ._linear_operator import ( AdjointLinearOperator, BinaryOperandType, + Embedding, Identity, LinearOperator, Matrix, @@ -46,6 +47,7 @@ _TypeCastLinearOperator, Scaling, Selection, + Embedding, Zero, Kronecker, ] @@ -81,7 +83,7 @@ def matmul(op1: LinearOperator, op2: LinearOperator) -> LinearOperator: ######################################################################################## _BinaryOperatorType = Callable[ - [LinearOperator, LinearOperator], Union[LinearOperator, type(NotImplemented)] + [LinearOperator, LinearOperator], Union[LinearOperator, NotImplementedType] ] _BinaryOperatorRegistryType = Dict[Tuple[type, type], _BinaryOperatorType] @@ -91,13 +93,11 @@ def matmul(op1: LinearOperator, op2: LinearOperator) -> LinearOperator: _mul_fns: _BinaryOperatorRegistryType = {} _matmul_fns: _BinaryOperatorRegistryType = {} -# Scaling -_add_fns[(Scaling, Scaling)] = Scaling._add_scaling -_sub_fns[(Scaling, Scaling)] = Scaling._sub_scaling -_mul_fns[(Scaling, Scaling)] = Scaling._mul_scaling -_matmul_fns[(Scaling, Scaling)] = Scaling._matmul_scaling - +######################################################################################## +# Fill Arithmetics Registries +######################################################################################## +# Scaling def _mul_scalar_scaling(scalar: ScalarArgType, scaling: Scaling) -> Scaling: if scaling.is_isotropic: return Scaling(scalar * scaling.scalar, shape=scaling.shape) @@ -114,6 +114,10 @@ def _mul_scaling_scalar(scaling: Scaling, scalar: ScalarArgType) -> Scaling: _mul_fns[("scalar", Scaling)] = _mul_scalar_scaling _mul_fns[(Scaling, "scalar")] = _mul_scaling_scalar +_add_fns[(Scaling, Scaling)] = Scaling._add_scaling +_sub_fns[(Scaling, Scaling)] = Scaling._sub_scaling +_mul_fns[(Scaling, Scaling)] = Scaling._mul_scaling +_matmul_fns[(Scaling, Scaling)] = Scaling._matmul_scaling # ScaledLinearOperator def _matmul_scaled_op(scaled, anylinop): @@ -127,40 +131,64 @@ def _matmul_op_scaled(anylinop, scaled): for op_type in _AnyLinOp: _matmul_fns[(ScaledLinearOperator, op_type)] = _matmul_scaled_op _matmul_fns[(op_type, ScaledLinearOperator)] = _matmul_op_scaled + _matmul_fns[(NegatedLinearOperator, op_type)] = _matmul_scaled_op + _matmul_fns[(op_type, NegatedLinearOperator)] = _matmul_op_scaled # Kronecker -_add_fns[(Kronecker, Kronecker)] = Kronecker._add_kronecker -_sub_fns[(Kronecker, Kronecker)] = Kronecker._sub_kronecker - def _matmul_scaling_kronecker(scaling: Scaling, kronecker: Kronecker) -> Kronecker: if scaling.is_isotropic: - if ("scalar", type(kronecker.A)) in _mul_fns: - return Kronecker(A=scaling.scalar * kronecker.A, B=kronecker.B) - if ("scalar", type(kronecker.B)) in _mul_fns: - return Kronecker(A=kronecker.A, B=scaling.scalar * kronecker.B) + return scaling.scalar * kronecker return NotImplemented def _matmul_kronecker_scaling(kronecker: Kronecker, scaling: Scaling) -> Kronecker: if scaling.is_isotropic: - if (type(kronecker.B), "scalar") in _mul_fns: - return Kronecker(A=kronecker.A, B=kronecker.B * scaling.scalar) - if (type(kronecker.A), "scalar") in _mul_fns: - return Kronecker(A=kronecker.A, B=kronecker.B * scaling.scalar) + return kronecker * scaling.scalar return NotImplemented +def _mul_scalar_kronecker(scalar: ScalarArgType, kronecker: Kronecker) -> Kronecker: + + prefer_A = ("scalar", type(kronecker.A)) in _mul_fns + prefer_B = ("scalar", type(kronecker.B)) in _mul_fns + + if prefer_A and not prefer_B: + return Kronecker(A=scalar * kronecker.A, B=kronecker.B) + + if prefer_B and not prefer_A: + return Kronecker(A=kronecker.A, B=scalar * kronecker.B) + + return Kronecker(A=scalar * kronecker.A, B=kronecker.B) + + +def _mul_kronecker_scalar(kronecker: Kronecker, scalar: ScalarArgType) -> Kronecker: + + prefer_A = (type(kronecker.A), "scalar") in _mul_fns + prefer_B = (type(kronecker.B), "scalar") in _mul_fns + + if prefer_A and not prefer_B: + return Kronecker(A=kronecker.A * scalar, B=kronecker.B) + + if prefer_B and not prefer_A: + return Kronecker(A=kronecker.A, B=kronecker.B * scalar) + + return Kronecker(A=kronecker.A, B=kronecker.B * scalar) + + _matmul_fns[(Kronecker, Kronecker)] = Kronecker._matmul_kronecker _mul_fns[(Kronecker, Kronecker)] = Kronecker._mul_kronecker +_add_fns[(Kronecker, Kronecker)] = Kronecker._add_kronecker +_sub_fns[(Kronecker, Kronecker)] = Kronecker._sub_kronecker -_mul_fns[("scalar", Kronecker)] = lambda sc, kr: Kronecker(A=sc * kr.A, B=kr.B) -_mul_fns[(Kronecker, "scalar")] = lambda kr, sc: Kronecker(A=kr.A, B=kr.B * sc) +_mul_fns[("scalar", Kronecker)] = _mul_scalar_kronecker +_mul_fns[(Kronecker, "scalar")] = _mul_kronecker_scalar _matmul_fns[(Kronecker, Scaling)] = _matmul_kronecker_scaling _matmul_fns[(Scaling, Kronecker)] = _matmul_scaling_kronecker + # Matrix def _matmul_scaling_matrix(scaling: Scaling, matrix: Matrix) -> Matrix: if scaling.shape[1] != matrix.shape[0]: @@ -176,14 +204,14 @@ def _matmul_matrix_scaling(matrix: Matrix, scaling: Scaling) -> Matrix: return Matrix(A=np.multiply(matrix.A, scaling.factors)) -def _mul_matrix_scalar(mat, scalar) -> Union[type(NotImplemented), Matrix]: +def _mul_matrix_scalar(mat, scalar) -> Union[NotImplementedType, Matrix]: if np.isscalar(scalar): return Matrix(A=scalar * mat.A) return NotImplemented -def _mul_scalar_matrix(scalar, mat) -> Union[type(NotImplemented), Matrix]: +def _mul_scalar_matrix(scalar, mat) -> Union[NotImplementedType, Matrix]: if np.isscalar(scalar): return Matrix(A=scalar * mat.A) @@ -192,7 +220,6 @@ def _mul_scalar_matrix(scalar, mat) -> Union[type(NotImplemented), Matrix]: _mul_fns[(Matrix, "scalar")] = _mul_matrix_scalar _mul_fns[("scalar", Matrix)] = _mul_scalar_matrix - _matmul_fns[(Matrix, Matrix)] = Matrix._matmul_matrix _matmul_fns[(Scaling, Matrix)] = _matmul_scaling_matrix _matmul_fns[(Matrix, Scaling)] = _matmul_matrix_scaling @@ -203,12 +230,31 @@ def _mul_scalar_matrix(scalar, mat) -> Union[type(NotImplemented), Matrix]: _matmul_fns[(Identity, op_type)] = lambda idty, other: other _matmul_fns[(op_type, Identity)] = lambda other, idty: other -_mul_fns[(Identity, "scalar")] = lambda idty, sc: Scaling(sc, shape=idty.shape) -_mul_fns[("scalar", Identity)] = lambda sc, idty: Scaling(sc, shape=idty.shape) +# Selection / Embedding +def _matmul_selection_embedding( + selection: Selection, embedding: Embedding +) -> Union[NotImplementedType, Identity]: + + if (embedding.shape[-1] == selection.shape[-2]) and np.all( + selection.indices == embedding._put_indices + ): + return Identity(shape=(selection.shape[-2], embedding.shape[-1])) -# Selection -_matmul_fns[(Selection, Selection)] = Selection._matmul_selection + return NotImplemented + + +def _matmul_embedding_selection( + embedding: Embedding, selection: Selection +) -> Union[NotImplementedType, Selection]: + if (selection.shape[-1] == embedding.shape[-2]) and np.all( + embedding.put_indices == selection.indices + ): + return Selection() + + +_matmul_fns[(Selection, Embedding)] = _matmul_selection_embedding +# Embedding @ Selection would be Projection # Zero def _matmul_zero_anylinop(z: Zero, op: LinearOperator) -> Zero: @@ -239,12 +285,32 @@ def _add_anylinop_zero(op: LinearOperator, z: Zero) -> Zero: return op +def _sub_zero_anylinop(z: Zero, op: LinearOperator) -> Zero: + if z.shape != op.shape: + raise ValueError(f"shape mismatch") # TODO + + return -op + + +def _sub_anylinop_zero(op: LinearOperator, z: Zero) -> Zero: + if z.shape != op.shape: + raise ValueError(f"shape mismatch") # TODO + + return op + + for op_type in _AnyLinOp: _matmul_fns[(Zero, op_type)] = _matmul_zero_anylinop _matmul_fns[(op_type, Zero)] = _matmul_anylinop_zero _add_fns[(Zero, op_type)] = _add_zero_anylinop _add_fns[(op_type, Zero)] = _add_anylinop_zero - # TODO scalar mult + _sub_fns[(Zero, op_type)] = _sub_zero_anylinop + _sub_fns[(op_type, Zero)] = _sub_anylinop_zero + + +######################################################################################## +# Apply +######################################################################################## def _apply( @@ -254,10 +320,10 @@ def _apply( fallback_operator: Optional[ Callable[ [LinearOperator, LinearOperator], - Union[LinearOperator, type(NotImplemented)], + Union[LinearOperator, NotImplementedType], ] ] = None, -) -> Union[LinearOperator, type(NotImplemented)]: +) -> Union[LinearOperator, NotImplementedType]: if np.isscalar(op1): key1 = "scalar" else: diff --git a/src/probnum/linops/_arithmetic_fallbacks.py b/src/probnum/linops/_arithmetic_fallbacks.py index 7833b486c..1f06f40b9 100644 --- a/src/probnum/linops/_arithmetic_fallbacks.py +++ b/src/probnum/linops/_arithmetic_fallbacks.py @@ -5,7 +5,7 @@ import numpy as np import probnum.utils -from probnum.typing import ScalarArgType +from probnum.typing import NotImplementedType, ScalarArgType from ._linear_operator import BinaryOperandType, LinearOperator @@ -47,6 +47,15 @@ def _inv(self) -> "ScaledLinearOperator": return ScaledLinearOperator(self._linop.inv(), 1.0 / self._scalar) + def __mul__(self, other: BinaryOperandType) -> "LinearOperator": + if np.isscalar(other): + return ScaledLinearOperator(linop=self._linop, scalar=self._scalar * other) + + return super().__mul__(other) + + def __repr__(self) -> str: + return f"{self._scalar} * {self._linop}" + class NegatedLinearOperator(ScaledLinearOperator): def __init__(self, linop: LinearOperator): @@ -55,6 +64,9 @@ def __init__(self, linop: LinearOperator): def __neg__(self) -> "LinearOperator": return self._linop + def __repr__(self) -> str: + return f"-{self._linop}" + class SumLinearOperator(LinearOperator): """Sum of two linear operators.""" @@ -100,6 +112,12 @@ def __init__(self, *summands: LinearOperator): def __neg__(self): return SumLinearOperator(*(-summand for summand in self._summands)) + def __repr__(self): + res = "SumLinearOperator [\n" + for s in self._summands: + res += f"\t{s}, \n" + return res + "]" + @staticmethod def _expand_sum_ops(*summands: LinearOperator) -> Tuple[LinearOperator, ...]: expanded_summands = [] @@ -115,7 +133,7 @@ def _expand_sum_ops(*summands: LinearOperator) -> Tuple[LinearOperator, ...]: def _mul_fallback( op1: BinaryOperandType, op2: BinaryOperandType -) -> Union[LinearOperator, type(NotImplemented)]: +) -> Union[LinearOperator, NotImplementedType]: res = NotImplemented if isinstance(op1, LinearOperator) and isinstance(op2, LinearOperator): @@ -195,10 +213,16 @@ def _expand_prod_ops(*factors: LinearOperator) -> Tuple[LinearOperator, ...]: return tuple(expanded_factors) + def __repr__(self): + res = "ProductLinearOperator [\n" + for s in self._factors: + res += f"\t{s}, \n" + return res + "]" + def _matmul_fallback( op1: BinaryOperandType, op2: BinaryOperandType -) -> Union[LinearOperator, type(NotImplemented)]: +) -> Union[LinearOperator, NotImplementedType]: res = NotImplemented if isinstance(op1, LinearOperator) and isinstance(op2, LinearOperator): diff --git a/src/probnum/linops/_kronecker.py b/src/probnum/linops/_kronecker.py index 6a1dfa149..138690a0a 100644 --- a/src/probnum/linops/_kronecker.py +++ b/src/probnum/linops/_kronecker.py @@ -3,7 +3,7 @@ import numpy as np -from probnum.typing import DTypeArgType +from probnum.typing import DTypeArgType, NotImplementedType from . import _linear_operator, _utils from ._scaling import Scaling @@ -188,7 +188,7 @@ def _mul_kronecker(self, other: "Kronecker") -> "Kronecker": def _add_kronecker( self, other: "Kronecker" - ) -> Union[type(NotImplemented), "Kronecker"]: + ) -> Union[NotImplementedType, "Kronecker"]: if self.A == other.A: return Kronecker(A=self.A, B=self.B + other.B) @@ -200,7 +200,7 @@ def _add_kronecker( def _sub_kronecker( self, other: "Kronecker" - ) -> Union[type(NotImplemented), "Kronecker"]: + ) -> Union[NotImplementedType, "Kronecker"]: if self.A == other.A and self.B == other.B: return Kronecker( From 2ee287f791649bc77cbca12d6714c6d18e1ab5f1 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Mon, 2 Aug 2021 12:51:57 +0200 Subject: [PATCH 25/74] Matrix arithmetics --- src/probnum/linops/_arithmetic.py | 12 ++++++++++++ src/probnum/linops/_linear_operator.py | 3 +++ 2 files changed, 15 insertions(+) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index 55a7765ce..1e734d358 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -225,6 +225,18 @@ def _mul_scalar_matrix(scalar, mat) -> Union[NotImplementedType, Matrix]: _matmul_fns[(Matrix, Scaling)] = _matmul_matrix_scaling +_matmul_fns[(Selection, Matrix)] = lambda sel, mat: Matrix(sel @ mat.A) +_matmul_fns[(Embedding, Matrix)] = lambda emb, mat: Matrix(emb @ mat.A) +_matmul_fns[(Matrix, Selection)] = lambda mat, sel: Matrix(mat.A @ sel) +_matmul_fns[(Matrix, Embedding)] = lambda mat, emb: Matrix(mat.A @ emb) + +_add_fns[(Matrix, Matrix)] = lambda mat1, mat2: Matrix(mat1.A + mat2.A) +_sub_fns[(Matrix, Matrix)] = lambda mat1, mat2: Matrix(mat1.A - mat2.A) + +_matmul_fns[(Matrix, _InverseLinearOperator)] = lambda mat, inv: Matrix(mat.A @ inv) +_matmul_fns[(_InverseLinearOperator, Matrix)] = lambda inv, mat: Matrix(inv @ mat.A) + + # Identity for op_type in _AnyLinOp: _matmul_fns[(Identity, op_type)] = lambda idty, other: other diff --git a/src/probnum/linops/_linear_operator.py b/src/probnum/linops/_linear_operator.py index b0d470b7e..908589ca2 100644 --- a/src/probnum/linops/_linear_operator.py +++ b/src/probnum/linops/_linear_operator.py @@ -1041,6 +1041,9 @@ def __eq__(self, other: LinearOperator) -> bool: return np.all(self.A == other.A) + def __neg__(self) -> "Matrix": + return Matrix(-self.A) + class Identity(LinearOperator): """The identity operator. From 760160afb59e6b9fecfe132f78386d7f8eab3b97 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Mon, 2 Aug 2021 14:05:31 +0200 Subject: [PATCH 26/74] Remove obsolete function stub --- src/probnum/linops/_arithmetic.py | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index 1e734d358..9d0b8d2ed 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -256,15 +256,6 @@ def _matmul_selection_embedding( return NotImplemented -def _matmul_embedding_selection( - embedding: Embedding, selection: Selection -) -> Union[NotImplementedType, Selection]: - if (selection.shape[-1] == embedding.shape[-2]) and np.all( - embedding.put_indices == selection.indices - ): - return Selection() - - _matmul_fns[(Selection, Embedding)] = _matmul_selection_embedding # Embedding @ Selection would be Projection From ab246875f9b4824456d3447027ce930f16313f38 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Mon, 2 Aug 2021 14:15:52 +0200 Subject: [PATCH 27/74] Fix pylint issues --- src/probnum/linops/_arithmetic.py | 12 ++++++------ src/probnum/linops/_linear_operator.py | 8 ++++---- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index 9d0b8d2ed..fe607f1fd 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -262,42 +262,42 @@ def _matmul_selection_embedding( # Zero def _matmul_zero_anylinop(z: Zero, op: LinearOperator) -> Zero: if z.shape[1] != op.shape[0]: - raise ValueError(f"shape mismatch") # TODO + raise ValueError("shape mismatch") # TODO return Zero(shape=(z.shape[0], op.shape[1])) def _matmul_anylinop_zero(op: LinearOperator, z: Zero) -> Zero: if z.shape[0] != op.shape[1]: - raise ValueError(f"shape mismatch") # TODO + raise ValueError("shape mismatch") # TODO return Zero(shape=(op.shape[0], z.shape[1])) def _add_zero_anylinop(z: Zero, op: LinearOperator) -> Zero: if z.shape != op.shape: - raise ValueError(f"shape mismatch") # TODO + raise ValueError("shape mismatch") # TODO return op def _add_anylinop_zero(op: LinearOperator, z: Zero) -> Zero: if z.shape != op.shape: - raise ValueError(f"shape mismatch") # TODO + raise ValueError("shape mismatch") # TODO return op def _sub_zero_anylinop(z: Zero, op: LinearOperator) -> Zero: if z.shape != op.shape: - raise ValueError(f"shape mismatch") # TODO + raise ValueError("shape mismatch") # TODO return -op def _sub_anylinop_zero(op: LinearOperator, z: Zero) -> Zero: if z.shape != op.shape: - raise ValueError(f"shape mismatch") # TODO + raise ValueError("shape mismatch") # TODO return op diff --git a/src/probnum/linops/_linear_operator.py b/src/probnum/linops/_linear_operator.py index 908589ca2..c17b3521d 100644 --- a/src/probnum/linops/_linear_operator.py +++ b/src/probnum/linops/_linear_operator.py @@ -7,7 +7,7 @@ import scipy.sparse.linalg import probnum.utils -from probnum.typing import DTypeArgType, NotImplementedType, ScalarArgType, ShapeArgType +from probnum.typing import DTypeArgType, ScalarArgType, ShapeArgType BinaryOperandType = Union[ "LinearOperator", ScalarArgType, np.ndarray, scipy.sparse.spmatrix @@ -1119,7 +1119,7 @@ def __eq__(self, other: LinearOperator) -> bool: class Selection(LinearOperator): def __init__(self, indices, shape): if np.ndim(indices) > 1: - raise ValueError(f"bla") # TODO + raise ValueError("bla") # TODO self._indices = probnum.utils.as_shape(indices) # TODO assert shape[0] <= shape[1] assert len(self._indices) == shape[0] @@ -1128,7 +1128,7 @@ def __init__(self, indices, shape): dtype=np.uint8, # TODO shape=shape, matmul=lambda x: _selection_matmul(self.indices, x), - todense=lambda: self._todense(), + todense=self._todense, transpose=lambda: Embedding( take_indices=np.arange(len(self._indices)), put_indices=self._indices, @@ -1159,7 +1159,7 @@ def __init__(self, take_indices, put_indices, shape, fill_value=0.0): dtype=np.uint8, # TODO shape=shape, matmul=lambda x: _embedding_matmul(self, x), - todense=lambda: self._todense(), + todense=self._todense, transpose=lambda: Selection( indices=put_indices, shape=(self.shape[1], self.shape[0]) ), From ab2d9ad49b1be72abc22d753392e1be8c7f49601 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Mon, 2 Aug 2021 14:32:42 +0200 Subject: [PATCH 28/74] Fix some ToDo's --- src/probnum/linops/_arithmetic.py | 14 +++--- src/probnum/linops/_linear_operator.py | 61 +++++++++++++++++--------- 2 files changed, 47 insertions(+), 28 deletions(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index fe607f1fd..3ce4ba99e 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -262,42 +262,42 @@ def _matmul_selection_embedding( # Zero def _matmul_zero_anylinop(z: Zero, op: LinearOperator) -> Zero: if z.shape[1] != op.shape[0]: - raise ValueError("shape mismatch") # TODO + raise ValueError(f"matmul received invalid shapes {z.shape} @ {op.shape}" return Zero(shape=(z.shape[0], op.shape[1])) def _matmul_anylinop_zero(op: LinearOperator, z: Zero) -> Zero: - if z.shape[0] != op.shape[1]: - raise ValueError("shape mismatch") # TODO + if op.shape[1] != z.shape[0]: + raise ValueError(f"matmul received invalid shapes {op.shape} @ {z.shape}" return Zero(shape=(op.shape[0], z.shape[1])) def _add_zero_anylinop(z: Zero, op: LinearOperator) -> Zero: if z.shape != op.shape: - raise ValueError("shape mismatch") # TODO + raise ValueError(f"add received invalid shapes {z.shape} + {op.shape}") return op def _add_anylinop_zero(op: LinearOperator, z: Zero) -> Zero: if z.shape != op.shape: - raise ValueError("shape mismatch") # TODO + raise ValueError(f"add received invalid shapes {op.shape} + {z.shape}") return op def _sub_zero_anylinop(z: Zero, op: LinearOperator) -> Zero: if z.shape != op.shape: - raise ValueError("shape mismatch") # TODO + raise ValueError(f"sub received invalid shapes {op.shape} - {z.shape}") return -op def _sub_anylinop_zero(op: LinearOperator, z: Zero) -> Zero: if z.shape != op.shape: - raise ValueError("shape mismatch") # TODO + raise ValueError(f"sub received invalid shapes {op.shape} - {z.shape}") return op diff --git a/src/probnum/linops/_linear_operator.py b/src/probnum/linops/_linear_operator.py index c17b3521d..e337f9c6b 100644 --- a/src/probnum/linops/_linear_operator.py +++ b/src/probnum/linops/_linear_operator.py @@ -1117,15 +1117,23 @@ def __eq__(self, other: LinearOperator) -> bool: class Selection(LinearOperator): - def __init__(self, indices, shape): + def __init__(self, indices, shape, dtype=np.double): if np.ndim(indices) > 1: - raise ValueError("bla") # TODO - self._indices = probnum.utils.as_shape(indices) # TODO - assert shape[0] <= shape[1] + raise ValueError( + "Selection LinOp expects an integer or (1D) iterable of " + f"integers. Received {type(indices)} with shape {np.shape(indices)}." + ) + if shape[0] > shape[1]: + raise ValueError( + f"Invalid shape {shape} for Selection LinOp. If the " + "output-dimension (shape[0]) is larger than the input-dimension " + "(shape[1]), consider using `Embedding`." + ) + self._indices = probnum.utils.as_shape(indices) assert len(self._indices) == shape[0] super().__init__( - dtype=np.uint8, # TODO + dtype=dtype, shape=shape, matmul=lambda x: _selection_matmul(self.indices, x), todense=self._todense, @@ -1150,13 +1158,35 @@ def _selection_matmul(indices, M): class Embedding(LinearOperator): - def __init__(self, take_indices, put_indices, shape, fill_value=0.0): - self._take_indices = probnum.utils.as_shape(take_indices) # TODO - self._put_indices = probnum.utils.as_shape(put_indices) # TODO - self._fill_value = fill_value # TODO + def __init__( + self, take_indices, put_indices, shape, fill_value=0.0, dtype=np.double + ): + if np.ndim(take_indices) > 1: + raise ValueError( + "Embedding LinOp expects an integer or (1D) iterable of " + f"integers. Received {type(take_indices)} with shape " + f"{np.shape(take_indices)}." + ) + if np.ndim(put_indices) > 1: + raise ValueError( + "Embedding LinOp expects an integer or (1D) iterable of " + f"integers. Received {type(put_indices)} with shape " + f"{np.shape(put_indices)}." + ) + + if shape[0] < shape[1]: + raise ValueError( + f"Invalid shape {shape} for Embedding LinOp. If the " + "output-dimension (shape[0]) is smaller than the input-dimension " + "(shape[1]), consider using `Selection`." + ) + + self._take_indices = probnum.utils.as_shape(take_indices) + self._put_indices = probnum.utils.as_shape(put_indices) + self._fill_value = fill_value super().__init__( - dtype=np.uint8, # TODO + dtype=dtype, shape=shape, matmul=lambda x: _embedding_matmul(self, x), todense=self._todense, @@ -1174,17 +1204,6 @@ def _embedding_matmul(embedding, M): res_shape = np.array(M.shape) res_shape[-2] = embedding.shape[0] res = np.full(shape=tuple(res_shape), fill_value=embedding._fill_value) - # print(f"Take indices {embedding._take_indices}") - # print(f"Put indices {embedding._put_indices}") - # print(f"Res, {res}") - # res[np.array(embedding._put_indices).reshape(-1, 1)] = M[embedding._take_indices] - # TODO - # take_from_M = np.take_along_axis( - # M, indices=np.array(embedding._put_indices), axis=-2 - # ) take_from_M = M[..., np.array(embedding._put_indices), :] - # np.put_along_axis( - # res, indices=np.array(embedding._put_indices), values=take_from_M, axis=-2 - # ) res[..., np.array(embedding._put_indices), :] = take_from_M return res From 9da5bbb0358f317a0a7021262e04c83b147caaea Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Mon, 2 Aug 2021 14:38:45 +0200 Subject: [PATCH 29/74] Add switch for matrix-matrix products in LinOps --- src/probnum/_config.py | 7 +++++++ src/probnum/linops/_arithmetic.py | 10 ++++++---- src/probnum/linops/_linear_operator.py | 11 +++++++---- 3 files changed, 20 insertions(+), 8 deletions(-) diff --git a/src/probnum/_config.py b/src/probnum/_config.py index 0cacc412a..7c5f07bcd 100644 --- a/src/probnum/_config.py +++ b/src/probnum/_config.py @@ -134,6 +134,13 @@ def register(self, key: str, default_value: Any, description: str) -> None: "of arrays. LinearOperators define a matrix-vector product implicitly without instantiating the full matrix in memory. This makes them memory- and runtime-efficient for linear algebra operations." ), ), + ( + "collapse_dense_matrix_linop_products", + True, + "If True, the matmul operation applied to two Matrix-LinearOperators will " + "again yield a Matrix-LinearOperator with the computed matrix-matrix product," + " instead of a ProductLinearOperator.", + ), ] # ... and register the default configuration options. diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index 3ce4ba99e..35dbd94b1 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -4,6 +4,7 @@ import numpy as np import scipy.sparse +from probnum import config from probnum.typing import NotImplementedType, ScalarArgType, ShapeArgType from ._arithmetic_fallbacks import ( @@ -233,8 +234,9 @@ def _mul_scalar_matrix(scalar, mat) -> Union[NotImplementedType, Matrix]: _add_fns[(Matrix, Matrix)] = lambda mat1, mat2: Matrix(mat1.A + mat2.A) _sub_fns[(Matrix, Matrix)] = lambda mat1, mat2: Matrix(mat1.A - mat2.A) -_matmul_fns[(Matrix, _InverseLinearOperator)] = lambda mat, inv: Matrix(mat.A @ inv) -_matmul_fns[(_InverseLinearOperator, Matrix)] = lambda inv, mat: Matrix(inv @ mat.A) +if config.collapse_dense_matrix_linop_products: + _matmul_fns[(Matrix, _InverseLinearOperator)] = lambda mat, inv: Matrix(mat.A @ inv) + _matmul_fns[(_InverseLinearOperator, Matrix)] = lambda inv, mat: Matrix(inv @ mat.A) # Identity @@ -262,14 +264,14 @@ def _matmul_selection_embedding( # Zero def _matmul_zero_anylinop(z: Zero, op: LinearOperator) -> Zero: if z.shape[1] != op.shape[0]: - raise ValueError(f"matmul received invalid shapes {z.shape} @ {op.shape}" + raise ValueError(f"matmul received invalid shapes {z.shape} @ {op.shape}") return Zero(shape=(z.shape[0], op.shape[1])) def _matmul_anylinop_zero(op: LinearOperator, z: Zero) -> Zero: if op.shape[1] != z.shape[0]: - raise ValueError(f"matmul received invalid shapes {op.shape} @ {z.shape}" + raise ValueError(f"matmul received invalid shapes {op.shape} @ {z.shape}") return Zero(shape=(op.shape[0], z.shape[1])) diff --git a/src/probnum/linops/_linear_operator.py b/src/probnum/linops/_linear_operator.py index e337f9c6b..f877c3da7 100644 --- a/src/probnum/linops/_linear_operator.py +++ b/src/probnum/linops/_linear_operator.py @@ -7,6 +7,7 @@ import scipy.sparse.linalg import probnum.utils +from probnum import config from probnum.typing import DTypeArgType, ScalarArgType, ShapeArgType BinaryOperandType = Union[ @@ -1029,11 +1030,13 @@ def _sparse_inv(self) -> "Matrix": raise np.linalg.LinAlgError(str(err)) from err def _matmul_matrix(self, other: "Matrix") -> "Matrix": - # TODO: Switch via config option - if not self.shape[1] == other.shape[0]: - raise ValueError(f"Matmul shape mismatch {self.shape} x {other.shape}") + if config.collapse_dense_matrix_linop_products: + if not self.shape[1] == other.shape[0]: + raise ValueError(f"Matmul shape mismatch {self.shape} x {other.shape}") - return Matrix(A=self.A @ other.A) + return Matrix(A=self.A @ other.A) + + return NotImplemented def __eq__(self, other: LinearOperator) -> bool: if not self._is_type_shape_dtype_equal(other): From 31ebf335d9e38b9eb875d3d36172adda66d05144 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Mon, 2 Aug 2021 21:34:42 +0200 Subject: [PATCH 30/74] Add IdentityKronecker operator --- src/probnum/linops/__init__.py | 2 +- src/probnum/linops/_arithmetic.py | 47 ++++++++++++++++- src/probnum/linops/_kronecker.py | 87 +++++++++++++++++++++++++++++++ 3 files changed, 134 insertions(+), 2 deletions(-) diff --git a/src/probnum/linops/__init__.py b/src/probnum/linops/__init__.py index f12ae2fce..0aba7d0ab 100644 --- a/src/probnum/linops/__init__.py +++ b/src/probnum/linops/__init__.py @@ -12,7 +12,7 @@ :class:`~probnum.linops.LinearOperator` instances. """ -from ._kronecker import Kronecker, SymmetricKronecker, Symmetrize +from ._kronecker import IdentityKronecker, Kronecker, SymmetricKronecker, Symmetrize from ._linear_operator import Embedding, Identity, LinearOperator, Matrix, Selection from ._scaling import Scaling, Zero from ._utils import LinearOperatorLike, aslinop diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index 35dbd94b1..d30b3ce62 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -15,7 +15,7 @@ _matmul_fallback, _mul_fallback, ) -from ._kronecker import Kronecker, SymmetricKronecker, Symmetrize +from ._kronecker import IdentityKronecker, Kronecker, SymmetricKronecker, Symmetrize from ._linear_operator import ( AdjointLinearOperator, BinaryOperandType, @@ -190,6 +190,51 @@ def _mul_kronecker_scalar(kronecker: Kronecker, scalar: ScalarArgType) -> Kronec _matmul_fns[(Scaling, Kronecker)] = _matmul_scaling_kronecker +# IdentityKronecker + + +def _matmul_scaling_idkronecker( + scaling: Scaling, idkronecker: IdentityKronecker +) -> IdentityKronecker: + if scaling.is_isotropic: + return scaling.scalar * idkronecker + return NotImplemented + + +def _matmul_idkronecker_scaling( + idkronecker: IdentityKronecker, scaling: Scaling +) -> IdentityKronecker: + if scaling.is_isotropic: + return idkronecker * scaling.scalar + return NotImplemented + + +def _mul_scalar_idkronecker( + scalar: ScalarArgType, idkronecker: IdentityKronecker +) -> IdentityKronecker: + + return IdentityKronecker(A=idkronecker.A, B=scalar * idkronecker.B) + + +def _mul_idkronecker_scalar( + idkronecker: IdentityKronecker, scalar: ScalarArgType +) -> IdentityKronecker: + + return IdentityKronecker(A=idkronecker.A, B=idkronecker.B * scalar) + + +_matmul_fns[ + (IdentityKronecker, IdentityKronecker) +] = IdentityKronecker._matmul_idkronecker +_mul_fns[(IdentityKronecker, IdentityKronecker)] = IdentityKronecker._mul_idkronecker +_add_fns[(IdentityKronecker, IdentityKronecker)] = IdentityKronecker._add_idkronecker +_sub_fns[(IdentityKronecker, IdentityKronecker)] = IdentityKronecker._sub_idkronecker + +_mul_fns[("scalar", IdentityKronecker)] = _mul_scalar_idkronecker +_mul_fns[(IdentityKronecker, "scalar")] = _mul_idkronecker_scalar +_matmul_fns[(IdentityKronecker, Scaling)] = _matmul_idkronecker_scaling +_matmul_fns[(Scaling, IdentityKronecker)] = _matmul_scaling_idkronecker + # Matrix def _matmul_scaling_matrix(scaling: Scaling, matrix: Matrix) -> Matrix: if scaling.shape[1] != matrix.shape[0]: diff --git a/src/probnum/linops/_kronecker.py b/src/probnum/linops/_kronecker.py index 138690a0a..12c5d1c2d 100644 --- a/src/probnum/linops/_kronecker.py +++ b/src/probnum/linops/_kronecker.py @@ -471,3 +471,90 @@ def _cond_identical_factors(self, p) -> np.inexact: return self.A.cond(p=p) * self.B.cond(p=p) return np.linalg.cond(self.todense(cache=False), p=p) + + +class IdentityKronecker(_linear_operator.LinearOperator): + def __init__(self, num_blocks: int, B: _utils.LinearOperatorLike): + self.A = _linear_operator.Identity(num_blocks) + self.B = _utils.aslinop(B) + + if self.B.is_square: + # det(A (x) B) = det(A)^n * det(B)^m + det = lambda: (self.B.det() ** self.A.shape[0]) + logabsdet = lambda: (self.A.shape[0] * self.B.logabsdet()) + trace = lambda: self.A.trace() * self.B.trace() + else: + det = None + logabsdet = None + trace = None + + super().__init__( + dtype=np.result_type(self.A.dtype, self.B.dtype), + shape=( + num_blocks * self.B.shape[0], + num_blocks * self.B.shape[1], + ), + matmul=lambda x: _kronecker_matmul( + self.A, self.B, x + ), # TODO: is there a more efficient way? + rmatmul=lambda x: _kronecker_rmatmul( + self.A, self.B, x + ), # TODO: is there a more efficient way? + todense=lambda: np.kron( + self.A.todense(cache=False), self.B.todense(cache=False) + ), + conjugate=lambda: Kronecker(A=self.A.conj(), B=self.B.conj()), + # (A (x) B)^T = A^T (x) B^T + transpose=lambda: Kronecker(A=self.A.T, B=self.B.T), + # (A (x) B)^H = A^H (x) B^H + adjoint=lambda: Kronecker(A=self.A.H, B=self.B.H), + # (A (x) B)^-1 = A^-1 (x) B^-1 + inverse=lambda: IdentityKronecker(num_blocks, B=self.B.inv()), + rank=lambda: self.A.rank() * self.B.rank(), + cond=self._cond, + det=det, + logabsdet=logabsdet, + trace=trace, + ) + + def _matmul_idkronecker(self, other: "IdentityKronecker") -> "IdentityKronecker": + _A, _B = self.A, self.B + _C, _D = other.A, other.B + if not (_A.shape == _C.shape and _B.shape[1] == _D.shape[0]): + raise ValueError( + f"Matmul shape mismatch {_A.shape} x {_C.shape} " + f"or {_B.shape} x {_D.shape}" + ) + + # Using (A (x) B) @ (C (x) D) = (A @ C) (x) (B @ D) + return IdentityKronecker(A=_A, B=_B @ _D) + + def _mul_idkronecker(self, other: "IdentityKronecker") -> "IdentityKronecker": + _A, _B = self.A, self.B + _C, _D = other.A, other.B + if not (_A.shape == _C.shape and _B.shape == _D.shape): + raise ValueError( + f"Matmul shape mismatch {_A.shape} x {_C.shape} " + f"or {_B.shape} x {_D.shape}" + ) + + # Using (A (x) B) o (C (x) D) = (A o C) (x) (B o D) + return IdentityKronecker(A=_A, B=_B * _D) + + def _add_idkronecker( + self, other: "IdentityKronecker" + ) -> Union[NotImplementedType, "IdentityKronecker"]: + + if self.A.shape == other.A.shape: + return IdentityKronecker(A=self.A, B=self.B + other.B) + + return NotImplemented + + def _sub_idkronecker( + self, other: "IdentityKronecker" + ) -> Union[NotImplementedType, "IdentityKronecker"]: + + if self.A.shape == other.A.shape: + return IdentityKronecker(A=self.A, B=self.B - other.B) + + return NotImplemented From 97597bb8a3e38f36dceef1bd0bad4bc23db11db1 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Tue, 3 Aug 2021 11:50:15 +0200 Subject: [PATCH 31/74] Fix normal rv addition test --- tests/test_randvars/test_normal.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_randvars/test_normal.py b/tests/test_randvars/test_normal.py index 9139eec7f..650905eac 100644 --- a/tests/test_randvars/test_normal.py +++ b/tests/test_randvars/test_normal.py @@ -57,7 +57,7 @@ def setUp(self): ), ( linops.Matrix(A=sparsemat.todense()), - linops.Kronecker(0.1 * linops.Identity(m), linops.Identity(n)), + linops.Kronecker(linops.Identity(m), linops.Identity(n)), ), ( linops.Matrix(A=np.random.uniform(size=(2, 2))), From d63eeb194834b311437e73f4eb58a30b116276dd Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Tue, 10 Aug 2021 11:14:06 +0200 Subject: [PATCH 32/74] Rename config options according to decision --- src/probnum/_config.py | 4 ++-- src/probnum/linops/_linear_operator.py | 2 +- .../markov/discrete/_linear_gaussian.py | 3 ++- .../randprocs/markov/integrator/_integrator.py | 11 +++++++++-- .../randprocs/markov/integrator/_iwp.py | 18 +++++++++--------- 5 files changed, 23 insertions(+), 15 deletions(-) diff --git a/src/probnum/_config.py b/src/probnum/_config.py index 7c5f07bcd..e80608f71 100644 --- a/src/probnum/_config.py +++ b/src/probnum/_config.py @@ -135,8 +135,8 @@ def register(self, key: str, default_value: Any, description: str) -> None: ), ), ( - "collapse_dense_matrix_linop_products", - True, + "lazy_matrix_matrix_matmul", + False, "If True, the matmul operation applied to two Matrix-LinearOperators will " "again yield a Matrix-LinearOperator with the computed matrix-matrix product," " instead of a ProductLinearOperator.", diff --git a/src/probnum/linops/_linear_operator.py b/src/probnum/linops/_linear_operator.py index f877c3da7..c83b7b896 100644 --- a/src/probnum/linops/_linear_operator.py +++ b/src/probnum/linops/_linear_operator.py @@ -1030,7 +1030,7 @@ def _sparse_inv(self) -> "Matrix": raise np.linalg.LinAlgError(str(err)) from err def _matmul_matrix(self, other: "Matrix") -> "Matrix": - if config.collapse_dense_matrix_linop_products: + if config.lazy_matrix_matrix_matmul: if not self.shape[1] == other.shape[0]: raise ValueError(f"Matmul shape mismatch {self.shape} x {other.shape}") diff --git a/src/probnum/randprocs/markov/discrete/_linear_gaussian.py b/src/probnum/randprocs/markov/discrete/_linear_gaussian.py index dcd21a9dd..d6c3a3830 100644 --- a/src/probnum/randprocs/markov/discrete/_linear_gaussian.py +++ b/src/probnum/randprocs/markov/discrete/_linear_gaussian.py @@ -173,7 +173,8 @@ def _forward_rv_classic( info = {"crosscov": crosscov} if compute_gain: if config.matrix_free: - gain = (new_cov.T.inv() @ crosscov.T).T + # gain = (new_cov.T.inv() @ crosscov.T).T + gain = crosscov @ new_cov.inv() else: gain = scipy.linalg.solve(new_cov.T, crosscov.T, assume_a="sym").T info["gain"] = gain diff --git a/src/probnum/randprocs/markov/integrator/_integrator.py b/src/probnum/randprocs/markov/integrator/_integrator.py index b10e46b72..79923e050 100644 --- a/src/probnum/randprocs/markov/integrator/_integrator.py +++ b/src/probnum/randprocs/markov/integrator/_integrator.py @@ -2,6 +2,7 @@ import numpy as np +from probnum import config, linops from probnum.randprocs.markov.integrator import _preconditioner __all__ = ["IntegratorTransition"] @@ -63,8 +64,12 @@ def proj2coord(self, coord: int) -> np.ndarray: """ projvec1d = np.eye(self.num_derivatives + 1)[:, coord] projmat1d = projvec1d.reshape((1, self.num_derivatives + 1)) - projmat = np.kron(np.eye(self.wiener_process_dimension), projmat1d) - return projmat + if config.matrix_free: + return linops.Kronecker( + linops.Identity(self.wiener_process_dimension), projmat1d + ) + + return np.kron(np.eye(self.wiener_process_dimension), projmat1d) @property def _derivwise2coordwise_projmat(self) -> np.ndarray: @@ -96,6 +101,8 @@ def _derivwise2coordwise_projmat(self) -> np.ndarray: * self.wiener_process_dimension : (q + 1) * self.wiener_process_dimension ] + if config.matrix_free: + projmat = linops.aslinop(projmat) return projmat @property diff --git a/src/probnum/randprocs/markov/integrator/_iwp.py b/src/probnum/randprocs/markov/integrator/_iwp.py index 1c3860194..fc7c09b66 100644 --- a/src/probnum/randprocs/markov/integrator/_iwp.py +++ b/src/probnum/randprocs/markov/integrator/_iwp.py @@ -138,9 +138,9 @@ def __init__( def _drift_matrix(self): drift_matrix_1d = np.diag(np.ones(self.num_derivatives), 1) if config.matrix_free: - return linops.Kronecker( - A=linops.Identity(self.wiener_process_dimension), - B=linops.Matrix(A=drift_matrix_1d), + return linops.IdentityKronecker( + num_blocks=self.wiener_process_dimension, + B=drift_matrix_1d, ) return np.kron(np.eye(self.wiener_process_dimension), drift_matrix_1d) @@ -154,9 +154,9 @@ def _dispersion_matrix(self): dispersion_matrix_1d[-1] = 1.0 # Unit diffusion if config.matrix_free: - return linops.Kronecker( - A=linops.Identity(self.wiener_process_dimension), - B=linops.Matrix(A=dispersion_matrix_1d.reshape(-1, 1)), + return linops.IdentityKronecker( + num_blocks=self.wiener_process_dimension, + B=dispersion_matrix_1d.reshape(-1, 1), ) return np.kron(np.eye(self.wiener_process_dimension), dispersion_matrix_1d).T @@ -177,7 +177,7 @@ def equivalent_discretisation_preconditioned(self): if config.matrix_free: state_transition = linops.Kronecker( A=linops.Identity(self.wiener_process_dimension), - B=linops.aslinop(state_transition_1d), + B=state_transition_1d, ) else: state_transition = np.kron( @@ -187,7 +187,7 @@ def equivalent_discretisation_preconditioned(self): if config.matrix_free: process_noise = linops.Kronecker( A=linops.Identity(self.wiener_process_dimension), - B=linops.aslinop(process_noise_1d), + B=process_noise_1d, ) else: process_noise = np.kron( @@ -201,7 +201,7 @@ def equivalent_discretisation_preconditioned(self): if config.matrix_free: process_noise_cholesky = linops.Kronecker( A=linops.Identity(self.wiener_process_dimension), - B=linops.aslinop(process_noise_cholesky_1d), + B=process_noise_cholesky_1d, ) else: process_noise_cholesky = np.kron( From 5b7aad2acc890260b922033a3cca6c99769e03c8 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Tue, 10 Aug 2021 11:57:43 +0200 Subject: [PATCH 33/74] Fix IdentityKronecker arithmetics --- src/probnum/linops/__init__.py | 8 +++++ src/probnum/linops/_arithmetic.py | 30 +++++++++++++++---- src/probnum/linops/_kronecker.py | 27 ++++++++++++----- .../randprocs/markov/integrator/_iwp.py | 15 ++++------ .../markov/integrator/_preconditioner.py | 5 ++-- 5 files changed, 60 insertions(+), 25 deletions(-) diff --git a/src/probnum/linops/__init__.py b/src/probnum/linops/__init__.py index 0aba7d0ab..33a11cdd0 100644 --- a/src/probnum/linops/__init__.py +++ b/src/probnum/linops/__init__.py @@ -20,22 +20,30 @@ # Public classes and functions. Order is reflected in documentation. __all__ = [ "aslinop", + "Embedding", "LinearOperator", "Matrix", "Identity", + "IdentityKronecker", "Scaling", "Kronecker", + "Selection", "SymmetricKronecker", "Symmetrize", + "Zero", ] # Set correct module paths. Corrects links and module paths in documentation. LinearOperator.__module__ = "probnum.linops" +Embedding.__module__ = "probnum.linops" Matrix.__module__ = "probnum.linops" Identity.__module__ = "probnum.linops" +IdentityKronecker.__module__ = "probnum.linops" Scaling.__module__ = "probnum.linops" Kronecker.__module__ = "probnum.linops" +Selection.__module__ = "probnum.linops" SymmetricKronecker.__module__ = "probnum.linops" Symmetrize.__module__ = "probnum.linops" +Zero.__module__ = "probnum.linops" aslinop.__module__ = "probnum.linops" diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index d30b3ce62..dbe155436 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -38,6 +38,7 @@ SumLinearOperator, AdjointLinearOperator, Identity, + IdentityKronecker, LinearOperator, Matrix, TransposedLinearOperator, @@ -213,14 +214,18 @@ def _mul_scalar_idkronecker( scalar: ScalarArgType, idkronecker: IdentityKronecker ) -> IdentityKronecker: - return IdentityKronecker(A=idkronecker.A, B=scalar * idkronecker.B) + return IdentityKronecker( + num_blocks=idkronecker.num_blocks, B=scalar * idkronecker.B + ) def _mul_idkronecker_scalar( idkronecker: IdentityKronecker, scalar: ScalarArgType ) -> IdentityKronecker: - return IdentityKronecker(A=idkronecker.A, B=idkronecker.B * scalar) + return IdentityKronecker( + num_blocks=idkronecker.num_blocks, B=idkronecker.B * scalar + ) _matmul_fns[ @@ -235,6 +240,9 @@ def _mul_idkronecker_scalar( _matmul_fns[(IdentityKronecker, Scaling)] = _matmul_idkronecker_scaling _matmul_fns[(Scaling, IdentityKronecker)] = _matmul_scaling_idkronecker +_matmul_fns[(Kronecker, IdentityKronecker)] = Kronecker._matmul_kronecker +_matmul_fns[(IdentityKronecker, Kronecker)] = Kronecker._matmul_kronecker + # Matrix def _matmul_scaling_matrix(scaling: Scaling, matrix: Matrix) -> Matrix: if scaling.shape[1] != matrix.shape[0]: @@ -279,9 +287,21 @@ def _mul_scalar_matrix(scalar, mat) -> Union[NotImplementedType, Matrix]: _add_fns[(Matrix, Matrix)] = lambda mat1, mat2: Matrix(mat1.A + mat2.A) _sub_fns[(Matrix, Matrix)] = lambda mat1, mat2: Matrix(mat1.A - mat2.A) -if config.collapse_dense_matrix_linop_products: - _matmul_fns[(Matrix, _InverseLinearOperator)] = lambda mat, inv: Matrix(mat.A @ inv) - _matmul_fns[(_InverseLinearOperator, Matrix)] = lambda inv, mat: Matrix(inv @ mat.A) + +def _matmul_matrix_inverse(mat: Matrix, inverse: _InverseLinearOperator) -> Matrix: + if config.lazy_matrix_matrix_matmul: + return Matrix(mat.A @ inverse) + return NotImplemented + + +def _matmul_inverse_matrix(inverse: _InverseLinearOperator, mat: Matrix) -> Matrix: + if config.lazy_matrix_matrix_matmul: + return Matrix(mat.A @ inverse) + return NotImplemented + + +_matmul_fns[(Matrix, _InverseLinearOperator)] = _matmul_matrix_inverse +_matmul_fns[(_InverseLinearOperator, Matrix)] = _matmul_inverse_matrix # Identity diff --git a/src/probnum/linops/_kronecker.py b/src/probnum/linops/_kronecker.py index 12c5d1c2d..caca6ae23 100644 --- a/src/probnum/linops/_kronecker.py +++ b/src/probnum/linops/_kronecker.py @@ -475,6 +475,7 @@ def _cond_identical_factors(self, p) -> np.inexact: class IdentityKronecker(_linear_operator.LinearOperator): def __init__(self, num_blocks: int, B: _utils.LinearOperatorLike): + self._num_blocks = num_blocks self.A = _linear_operator.Identity(num_blocks) self.B = _utils.aslinop(B) @@ -503,11 +504,11 @@ def __init__(self, num_blocks: int, B: _utils.LinearOperatorLike): todense=lambda: np.kron( self.A.todense(cache=False), self.B.todense(cache=False) ), - conjugate=lambda: Kronecker(A=self.A.conj(), B=self.B.conj()), - # (A (x) B)^T = A^T (x) B^T - transpose=lambda: Kronecker(A=self.A.T, B=self.B.T), + conjugate=lambda: IdentityKronecker(self._num_blocks, B=self.B.conj()), + # (A (x) B)^T = A (x) B^T + transpose=lambda: IdentityKronecker(self._num_blocks, B=self.B.T), # (A (x) B)^H = A^H (x) B^H - adjoint=lambda: Kronecker(A=self.A.H, B=self.B.H), + adjoint=lambda: IdentityKronecker(self._num_blocks, B=self.B.H), # (A (x) B)^-1 = A^-1 (x) B^-1 inverse=lambda: IdentityKronecker(num_blocks, B=self.B.inv()), rank=lambda: self.A.rank() * self.B.rank(), @@ -517,6 +518,10 @@ def __init__(self, num_blocks: int, B: _utils.LinearOperatorLike): trace=trace, ) + @property + def num_blocks(self): + return self._num_blocks + def _matmul_idkronecker(self, other: "IdentityKronecker") -> "IdentityKronecker": _A, _B = self.A, self.B _C, _D = other.A, other.B @@ -527,7 +532,7 @@ def _matmul_idkronecker(self, other: "IdentityKronecker") -> "IdentityKronecker" ) # Using (A (x) B) @ (C (x) D) = (A @ C) (x) (B @ D) - return IdentityKronecker(A=_A, B=_B @ _D) + return IdentityKronecker(num_blocks=self._num_blocks, B=_B @ _D) def _mul_idkronecker(self, other: "IdentityKronecker") -> "IdentityKronecker": _A, _B = self.A, self.B @@ -539,14 +544,14 @@ def _mul_idkronecker(self, other: "IdentityKronecker") -> "IdentityKronecker": ) # Using (A (x) B) o (C (x) D) = (A o C) (x) (B o D) - return IdentityKronecker(A=_A, B=_B * _D) + return IdentityKronecker(num_blocks=self._num_blocks, B=_B * _D) def _add_idkronecker( self, other: "IdentityKronecker" ) -> Union[NotImplementedType, "IdentityKronecker"]: if self.A.shape == other.A.shape: - return IdentityKronecker(A=self.A, B=self.B + other.B) + return IdentityKronecker(num_blocks=self._num_blocks, B=self.B + other.B) return NotImplemented @@ -555,6 +560,12 @@ def _sub_idkronecker( ) -> Union[NotImplementedType, "IdentityKronecker"]: if self.A.shape == other.A.shape: - return IdentityKronecker(A=self.A, B=self.B - other.B) + return IdentityKronecker(num_blocks=self._num_blocks, B=self.B - other.B) return NotImplemented + + def _cond(self, p) -> np.inexact: + if p is None or p in (2, 1, np.inf, "fro", -2, -1, -np.inf): + return self.A.cond(p=p) * self.B.cond(p=p) + + return np.linalg.cond(self.todense(cache=False), p=p) diff --git a/src/probnum/randprocs/markov/integrator/_iwp.py b/src/probnum/randprocs/markov/integrator/_iwp.py index fc7c09b66..76e3d2ca7 100644 --- a/src/probnum/randprocs/markov/integrator/_iwp.py +++ b/src/probnum/randprocs/markov/integrator/_iwp.py @@ -175,9 +175,8 @@ def equivalent_discretisation_preconditioned(self): scipy.linalg.pascal(self.num_derivatives + 1, kind="lower", exact=False) ) if config.matrix_free: - state_transition = linops.Kronecker( - A=linops.Identity(self.wiener_process_dimension), - B=state_transition_1d, + state_transition = linops.IdentityKronecker( + num_blocks=self.wiener_process_dimension, B=state_transition_1d ) else: state_transition = np.kron( @@ -185,9 +184,8 @@ def equivalent_discretisation_preconditioned(self): ) process_noise_1d = np.flip(scipy.linalg.hilbert(self.num_derivatives + 1)) if config.matrix_free: - process_noise = linops.Kronecker( - A=linops.Identity(self.wiener_process_dimension), - B=process_noise_1d, + process_noise = linops.IdentityKronecker( + num_blocks=self.wiener_process_dimension, B=process_noise_1d ) else: process_noise = np.kron( @@ -199,9 +197,8 @@ def equivalent_discretisation_preconditioned(self): process_noise_cholesky_1d = np.linalg.cholesky(process_noise_1d) if config.matrix_free: - process_noise_cholesky = linops.Kronecker( - A=linops.Identity(self.wiener_process_dimension), - B=process_noise_cholesky_1d, + process_noise_cholesky = linops.IdentityKronecker( + num_blocks=self.wiener_process_dimension, B=process_noise_cholesky_1d ) else: process_noise_cholesky = np.kron( diff --git a/src/probnum/randprocs/markov/integrator/_preconditioner.py b/src/probnum/randprocs/markov/integrator/_preconditioner.py index e6cb633a5..eb782af76 100644 --- a/src/probnum/randprocs/markov/integrator/_preconditioner.py +++ b/src/probnum/randprocs/markov/integrator/_preconditioner.py @@ -80,9 +80,8 @@ def from_order(cls, order, dimension): def __call__(self, step): scaling_vector = np.abs(step) ** self.powers / self.scales if config.matrix_free: - return linops.Kronecker( - A=linops.Identity(self.dimension), - B=linops.Scaling(factors=scaling_vector), + return linops.IdentityKronecker( + num_blocks=self.dimension, B=linops.Scaling(factors=scaling_vector) ) return np.kron(np.eye(self.dimension), np.diag(scaling_vector)) From 0338560f0762c2dcdb77be55bb703b83076a3772 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Tue, 10 Aug 2021 17:56:49 +0200 Subject: [PATCH 34/74] Half-baked arithmetics tests --- tests/test_linops/test_arithmetics.py | 112 ++++++++++++++++++++++ tests/test_linops/test_linearoperators.py | 38 -------- 2 files changed, 112 insertions(+), 38 deletions(-) create mode 100644 tests/test_linops/test_arithmetics.py delete mode 100644 tests/test_linops/test_linearoperators.py diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py new file mode 100644 index 000000000..fb96e75b7 --- /dev/null +++ b/tests/test_linops/test_arithmetics.py @@ -0,0 +1,112 @@ +"""Tests for linear operator arithmetics.""" +import numpy as np +import pytest + +from probnum import linops + +matrix_shape = (4, 4) + +kron_factor_shape = (2, 2) + + +def get_kronecker(): + return linops.Kronecker( + np.random.rand(*kron_factor_shape), np.random.rand(*kron_factor_shape) + ) + + +def get_zero(): + return linops.Zero(matrix_shape) + + +def get_scaling(): + return linops.Scaling(factors=np.random.rand(matrix_shape[0])) + + +def get_inverse_linop(): + rand_matrix = np.random.rand(*matrix_shape) + _mat = linops.Matrix(rand_matrix.T @ rand_matrix) + return _mat.inv() + + +def get_matrix_linop(): + return linops.Matrix(np.random.rand(*matrix_shape)) + + +def get_idkron_linop(): + return linops.IdentityKronecker( + num_blocks=kron_factor_shape[0], B=np.random.rand(*kron_factor_shape) + ) + + +def get_id_linop(): + return linops.Identity(matrix_shape[0]) + + +def get_sum_linop(): + summand1 = linops.Matrix(np.random.rand(*matrix_shape)) + summand2 = linops.Matrix(np.random.rand(*matrix_shape)) + return summand1 + summand2 + + +def get_product_linop(): + factor1 = linops.Matrix(np.random.rand(*matrix_shape)) + factor2 = linops.Matrix(np.random.rand(*matrix_shape)) + return factor1 @ factor2 + + +def get_negated_linop(): + op = linops.Matrix(np.random.rand(*matrix_shape)) + return -op + + +def get_scaled_linop(): + op = linops.Matrix(np.random.rand(*matrix_shape)) + return 3.14 * op + + +def get_transposed_linop(): + op = linops.Matrix(np.random.rand(*matrix_shape)) + return op.T + + +@pytest.mark.parametrize( + "linop1", + [ + get_kronecker(), + get_zero(), + get_scaling(), + get_inverse_linop(), + get_matrix_linop(), + get_idkron_linop(), + get_id_linop(), + get_sum_linop(), + get_product_linop(), + get_negated_linop(), + get_scaled_linop(), + get_transposed_linop(), + ], +) +@pytest.mark.parametrize( + "linop2", + [ + get_kronecker(), + get_zero(), + get_scaling(), + get_inverse_linop(), + get_matrix_linop(), + get_idkron_linop(), + get_id_linop(), + get_sum_linop(), + get_product_linop(), + get_negated_linop(), + get_scaled_linop(), + get_transposed_linop(), + ], +) +def test_arithmetics(linop1, linop2): + + res_linop = linop1 @ linop2 + assert res_linop.ndim == 2 + assert res_linop.shape[0] == linop1.shape[0] + assert res_linop.shape[1] == linop2.shape[1] diff --git a/tests/test_linops/test_linearoperators.py b/tests/test_linops/test_linearoperators.py deleted file mode 100644 index 63ca3a83c..000000000 --- a/tests/test_linops/test_linearoperators.py +++ /dev/null @@ -1,38 +0,0 @@ -"""Tests for linear operators.""" -import itertools -import unittest - -import numpy as np - -from probnum import linops -from tests.testing import NumpyAssertions - - -class LinearOperatorArithmeticTestCase(unittest.TestCase, NumpyAssertions): - """Test linear operator arithmetic.""" - - def setUp(self): - """Resources for tests.""" - # Random Seed - rng = np.random.default_rng(42) - - # Scalars and arrays - self.scalars = [0, int(1), 0.1, -4.2, np.nan, np.inf] - self.arrays = [rng.normal(size=[5, 4]), np.array([[3, 4], [1, 5]])] - - def test_scalar_mult(self): - """Matrix linear operator multiplication with scalars.""" - for A, alpha in list(itertools.product(self.arrays, self.scalars)): - with self.subTest(): - Aop = linops.Matrix(A) - - self.assertAllClose((alpha * Aop).todense(), alpha * A) - - def test_addition(self): - """Linear operator addition.""" - for A, B in list(zip(self.arrays, self.arrays)): - with self.subTest(): - Aop = linops.Matrix(A) - Bop = linops.Matrix(B) - - self.assertAllClose((Aop + Bop).todense(), A + B) From ec7f7ae4e7ef41a1140901287c032a4114b07de4 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Tue, 10 Aug 2021 18:03:40 +0200 Subject: [PATCH 35/74] Trigger build From e3977c57d0a64c98ba5406eb39bc2006a145a672 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Wed, 11 Aug 2021 16:55:14 +0200 Subject: [PATCH 36/74] Remove mul for Kronecker --- src/probnum/linops/_arithmetic.py | 3 --- src/probnum/linops/_kronecker.py | 24 ------------------------ 2 files changed, 27 deletions(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index dbe155436..e8a0d1227 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -39,7 +39,6 @@ AdjointLinearOperator, Identity, IdentityKronecker, - LinearOperator, Matrix, TransposedLinearOperator, SymmetricKronecker, @@ -181,7 +180,6 @@ def _mul_kronecker_scalar(kronecker: Kronecker, scalar: ScalarArgType) -> Kronec _matmul_fns[(Kronecker, Kronecker)] = Kronecker._matmul_kronecker -_mul_fns[(Kronecker, Kronecker)] = Kronecker._mul_kronecker _add_fns[(Kronecker, Kronecker)] = Kronecker._add_kronecker _sub_fns[(Kronecker, Kronecker)] = Kronecker._sub_kronecker @@ -231,7 +229,6 @@ def _mul_idkronecker_scalar( _matmul_fns[ (IdentityKronecker, IdentityKronecker) ] = IdentityKronecker._matmul_idkronecker -_mul_fns[(IdentityKronecker, IdentityKronecker)] = IdentityKronecker._mul_idkronecker _add_fns[(IdentityKronecker, IdentityKronecker)] = IdentityKronecker._add_idkronecker _sub_fns[(IdentityKronecker, IdentityKronecker)] = IdentityKronecker._sub_idkronecker diff --git a/src/probnum/linops/_kronecker.py b/src/probnum/linops/_kronecker.py index caca6ae23..748f774fa 100644 --- a/src/probnum/linops/_kronecker.py +++ b/src/probnum/linops/_kronecker.py @@ -174,18 +174,6 @@ def _matmul_kronecker(self, other: "Kronecker") -> "Kronecker": # Using (A (x) B) @ (C (x) D) = (A @ C) (x) (B @ D) return Kronecker(A=_A @ _C, B=_B @ _D) - def _mul_kronecker(self, other: "Kronecker") -> "Kronecker": - _A, _B = self.A, self.B - _C, _D = other.A, other.B - if not (_A.shape == _C.shape and _B.shape == _D.shape): - raise ValueError( - f"Matmul shape mismatch {_A.shape} x {_C.shape} " - f"or {_B.shape} x {_D.shape}" - ) - - # Using (A (x) B) o (C (x) D) = (A o C) (x) (B o D) - return Kronecker(A=_A * _C, B=_B * _D) - def _add_kronecker( self, other: "Kronecker" ) -> Union[NotImplementedType, "Kronecker"]: @@ -534,18 +522,6 @@ def _matmul_idkronecker(self, other: "IdentityKronecker") -> "IdentityKronecker" # Using (A (x) B) @ (C (x) D) = (A @ C) (x) (B @ D) return IdentityKronecker(num_blocks=self._num_blocks, B=_B @ _D) - def _mul_idkronecker(self, other: "IdentityKronecker") -> "IdentityKronecker": - _A, _B = self.A, self.B - _C, _D = other.A, other.B - if not (_A.shape == _C.shape and _B.shape == _D.shape): - raise ValueError( - f"Matmul shape mismatch {_A.shape} x {_C.shape} " - f"or {_B.shape} x {_D.shape}" - ) - - # Using (A (x) B) o (C (x) D) = (A o C) (x) (B o D) - return IdentityKronecker(num_blocks=self._num_blocks, B=_B * _D) - def _add_idkronecker( self, other: "IdentityKronecker" ) -> Union[NotImplementedType, "IdentityKronecker"]: From b0e2d96e10ae5482bbd01886b5404579ad6495e6 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Wed, 11 Aug 2021 16:55:24 +0200 Subject: [PATCH 37/74] Improve tests for linop arithmetics --- tests/test_linops/test_arithmetics.py | 218 ++++++++++++++------------ 1 file changed, 114 insertions(+), 104 deletions(-) diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py index fb96e75b7..f518d722e 100644 --- a/tests/test_linops/test_arithmetics.py +++ b/tests/test_linops/test_arithmetics.py @@ -1,112 +1,122 @@ """Tests for linear operator arithmetics.""" -import numpy as np -import pytest - -from probnum import linops - -matrix_shape = (4, 4) - -kron_factor_shape = (2, 2) - - -def get_kronecker(): - return linops.Kronecker( - np.random.rand(*kron_factor_shape), np.random.rand(*kron_factor_shape) - ) - - -def get_zero(): - return linops.Zero(matrix_shape) - - -def get_scaling(): - return linops.Scaling(factors=np.random.rand(matrix_shape[0])) - - -def get_inverse_linop(): - rand_matrix = np.random.rand(*matrix_shape) - _mat = linops.Matrix(rand_matrix.T @ rand_matrix) - return _mat.inv() - - -def get_matrix_linop(): - return linops.Matrix(np.random.rand(*matrix_shape)) - - -def get_idkron_linop(): - return linops.IdentityKronecker( - num_blocks=kron_factor_shape[0], B=np.random.rand(*kron_factor_shape) - ) +import operator -def get_id_linop(): - return linops.Identity(matrix_shape[0]) - - -def get_sum_linop(): - summand1 = linops.Matrix(np.random.rand(*matrix_shape)) - summand2 = linops.Matrix(np.random.rand(*matrix_shape)) - return summand1 + summand2 - - -def get_product_linop(): - factor1 = linops.Matrix(np.random.rand(*matrix_shape)) - factor2 = linops.Matrix(np.random.rand(*matrix_shape)) - return factor1 @ factor2 - - -def get_negated_linop(): - op = linops.Matrix(np.random.rand(*matrix_shape)) - return -op - - -def get_scaled_linop(): - op = linops.Matrix(np.random.rand(*matrix_shape)) - return 3.14 * op - +import numpy as np +import pytest -def get_transposed_linop(): - op = linops.Matrix(np.random.rand(*matrix_shape)) - return op.T +from probnum.linops._arithmetic import _add_fns, _matmul_fns, _mul_fns, _sub_fns +from probnum.linops._arithmetic_fallbacks import ( + NegatedLinearOperator, + ProductLinearOperator, + ScaledLinearOperator, + SumLinearOperator, +) +from probnum.linops._kronecker import ( + IdentityKronecker, + Kronecker, + SymmetricKronecker, + Symmetrize, +) +from probnum.linops._linear_operator import ( + AdjointLinearOperator, + BinaryOperandType, + Embedding, + Identity, + LinearOperator, + Matrix, + Selection, + TransposedLinearOperator, + _ConjugateLinearOperator, + _InverseLinearOperator, + _TypeCastLinearOperator, +) +from probnum.linops._scaling import Scaling, Zero +from probnum.problems.zoo.linalg import random_spd_matrix + + +def get_linop(linop_type): + if linop_type is Kronecker: + return linop_type(np.random.rand(2, 2), np.random.rand(2, 2)) + elif linop_type is IdentityKronecker: + return linop_type(2, np.random.rand(2, 2)) + elif linop_type is Zero or linop_type is Identity: + return linop_type(shape=(4, 4)) + elif linop_type is Scaling: + return linop_type(factors=np.random.rand(4)) + elif linop_type is Matrix: + return linop_type(np.random.rand(4, 4)) + elif linop_type is _InverseLinearOperator: + _posdef_randmat = random_spd_matrix(rng=np.random.default_rng(123), dim=4) + return Matrix(_posdef_randmat).inv() + elif linop_type is TransposedLinearOperator: + return TransposedLinearOperator(linop=Matrix(np.random.rand(4, 4))) + elif linop_type is Embedding: + return Embedding(take_indices=(0, 1, 2), put_indices=(1, 0, 3), shape=(4, 3)) + elif linop_type is Selection: + return Selection(indices=(1, 0, 3), shape=(3, 4)) + elif linop_type is NegatedLinearOperator: + return NegatedLinearOperator(linop=Matrix(np.random.rand(4, 4))) + elif linop_type is ScaledLinearOperator: + return ScaledLinearOperator(linop=Matrix(np.random.rand(4, 4)), scalar=3.14) + elif linop_type is ProductLinearOperator: + return ProductLinearOperator( + Matrix(np.random.rand(4, 4)), Matrix(np.random.rand(4, 4)) + ) + elif linop_type is SumLinearOperator: + return SumLinearOperator( + Matrix(np.random.rand(4, 4)), Matrix(np.random.rand(4, 4)) + ) + elif linop_type is AdjointLinearOperator: + return AdjointLinearOperator(linop=Identity(4)) + elif linop_type is _ConjugateLinearOperator: + return _ConjugateLinearOperator(linop=Identity(4, dtype=np.complex64)) + elif linop_type is SymmetricKronecker: + return SymmetricKronecker(Identity(2), Identity(2)) + elif linop_type is Symmetrize: + return Symmetrize(2) + elif linop_type is _TypeCastLinearOperator: + return _TypeCastLinearOperator( + linop=Matrix(np.random.rand(4, 4)), dtype=np.float32 + ) + elif isinstance(linop_type, str) and linop_type == "scalar": + return 1.3579 + else: + raise TypeError(f"Don't know what to do with type {linop_type}.") @pytest.mark.parametrize( - "linop1", - [ - get_kronecker(), - get_zero(), - get_scaling(), - get_inverse_linop(), - get_matrix_linop(), - get_idkron_linop(), - get_id_linop(), - get_sum_linop(), - get_product_linop(), - get_negated_linop(), - get_scaled_linop(), - get_transposed_linop(), - ], -) -@pytest.mark.parametrize( - "linop2", - [ - get_kronecker(), - get_zero(), - get_scaling(), - get_inverse_linop(), - get_matrix_linop(), - get_idkron_linop(), - get_id_linop(), - get_sum_linop(), - get_product_linop(), - get_negated_linop(), - get_scaled_linop(), - get_transposed_linop(), - ], + "op", [operator.matmul, operator.mul, operator.add, operator.sub] ) -def test_arithmetics(linop1, linop2): - - res_linop = linop1 @ linop2 - assert res_linop.ndim == 2 - assert res_linop.shape[0] == linop1.shape[0] - assert res_linop.shape[1] == linop2.shape[1] +def test_arithmetics(op): + + registry = { + operator.matmul: _matmul_fns, + operator.mul: _mul_fns, + operator.add: _add_fns, + operator.sub: _sub_fns, + }[op] + + for (l_type, r_type) in registry.keys(): + if ( + l_type is Selection + or l_type is Embedding + or r_type is Selection + or r_type is Embedding + ): + # Checked seperatly + continue + + linop1 = get_linop(l_type) + linop2 = get_linop(r_type) + + res_linop = op(linop1, linop2) + assert res_linop.ndim == 2 + + if isinstance(l_type, str) and l_type == "scalar": + assert res_linop.shape == linop2.shape + elif isinstance(r_type, str) and r_type == "scalar": + assert res_linop.shape == linop1.shape + else: + assert res_linop.shape[0] == linop1.shape[0] + assert res_linop.shape[1] == linop2.shape[1] From 3e9f10df6f0fb8cf0bf545d1c192c95da87644b0 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Wed, 11 Aug 2021 17:00:31 +0200 Subject: [PATCH 38/74] Fix pylint issues --- tests/test_linops/test_arithmetics.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py index f518d722e..d3075e845 100644 --- a/tests/test_linops/test_arithmetics.py +++ b/tests/test_linops/test_arithmetics.py @@ -20,10 +20,8 @@ ) from probnum.linops._linear_operator import ( AdjointLinearOperator, - BinaryOperandType, Embedding, Identity, - LinearOperator, Matrix, Selection, TransposedLinearOperator, @@ -36,6 +34,7 @@ def get_linop(linop_type): + # pylint: disable=too-many-return-statements if linop_type is Kronecker: return linop_type(np.random.rand(2, 2), np.random.rand(2, 2)) elif linop_type is IdentityKronecker: From 817baa07c973b8f47f72781e022412c22a756bf7 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Wed, 11 Aug 2021 17:22:53 +0200 Subject: [PATCH 39/74] Choose neutral ground for scalar times kronecker --- src/probnum/linops/_arithmetic.py | 26 ++++---------------------- 1 file changed, 4 insertions(+), 22 deletions(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index e8a0d1227..23f26bc99 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -152,31 +152,13 @@ def _matmul_kronecker_scaling(kronecker: Kronecker, scaling: Scaling) -> Kroneck def _mul_scalar_kronecker(scalar: ScalarArgType, kronecker: Kronecker) -> Kronecker: - - prefer_A = ("scalar", type(kronecker.A)) in _mul_fns - prefer_B = ("scalar", type(kronecker.B)) in _mul_fns - - if prefer_A and not prefer_B: - return Kronecker(A=scalar * kronecker.A, B=kronecker.B) - - if prefer_B and not prefer_A: - return Kronecker(A=kronecker.A, B=scalar * kronecker.B) - - return Kronecker(A=scalar * kronecker.A, B=kronecker.B) + sqrt_scalar = np.sqrt(scalar) + return Kronecker(A=sqrt_scalar * kronecker.A, B=sqrt_scalar * kronecker.B) def _mul_kronecker_scalar(kronecker: Kronecker, scalar: ScalarArgType) -> Kronecker: - - prefer_A = (type(kronecker.A), "scalar") in _mul_fns - prefer_B = (type(kronecker.B), "scalar") in _mul_fns - - if prefer_A and not prefer_B: - return Kronecker(A=kronecker.A * scalar, B=kronecker.B) - - if prefer_B and not prefer_A: - return Kronecker(A=kronecker.A, B=kronecker.B * scalar) - - return Kronecker(A=kronecker.A, B=kronecker.B * scalar) + sqrt_scalar = np.sqrt(scalar) + return Kronecker(A=sqrt_scalar * kronecker.A, B=sqrt_scalar * kronecker.B) _matmul_fns[(Kronecker, Kronecker)] = Kronecker._matmul_kronecker From dcbe5bf19f73c486f967ec4901222d9b78bf196c Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Wed, 11 Aug 2021 17:23:02 +0200 Subject: [PATCH 40/74] Add coverage --- tests/test_linops/test_arithmetics.py | 43 +++++++++++++++++---------- 1 file changed, 28 insertions(+), 15 deletions(-) diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py index d3075e845..03a8c45ca 100644 --- a/tests/test_linops/test_arithmetics.py +++ b/tests/test_linops/test_arithmetics.py @@ -1,5 +1,6 @@ """Tests for linear operator arithmetics.""" +import itertools import operator import numpy as np @@ -33,18 +34,28 @@ from probnum.problems.zoo.linalg import random_spd_matrix +def _aslist(arg): + try: + return list(arg) + except TypeError: + return [arg] + + def get_linop(linop_type): # pylint: disable=too-many-return-statements if linop_type is Kronecker: - return linop_type(np.random.rand(2, 2), np.random.rand(2, 2)) + return Kronecker(np.random.rand(2, 2), np.random.rand(2, 2)) elif linop_type is IdentityKronecker: - return linop_type(2, np.random.rand(2, 2)) + return IdentityKronecker(2, np.random.rand(2, 2)) elif linop_type is Zero or linop_type is Identity: return linop_type(shape=(4, 4)) elif linop_type is Scaling: - return linop_type(factors=np.random.rand(4)) + return ( + Scaling(factors=np.random.rand(4)), + Scaling(factors=3.14 * np.ones(4)), + ) elif linop_type is Matrix: - return linop_type(np.random.rand(4, 4)) + return Matrix(np.random.rand(4, 4)) elif linop_type is _InverseLinearOperator: _posdef_randmat = random_spd_matrix(rng=np.random.default_rng(123), dim=4) return Matrix(_posdef_randmat).inv() @@ -106,16 +117,18 @@ def test_arithmetics(op): # Checked seperatly continue - linop1 = get_linop(l_type) - linop2 = get_linop(r_type) + linops1 = get_linop(l_type) + linops2 = get_linop(r_type) + + for (linop1, linop2) in itertools.product(_aslist(linops1), _aslist(linops2)): - res_linop = op(linop1, linop2) - assert res_linop.ndim == 2 + res_linop = op(linop1, linop2) + assert res_linop.ndim == 2 - if isinstance(l_type, str) and l_type == "scalar": - assert res_linop.shape == linop2.shape - elif isinstance(r_type, str) and r_type == "scalar": - assert res_linop.shape == linop1.shape - else: - assert res_linop.shape[0] == linop1.shape[0] - assert res_linop.shape[1] == linop2.shape[1] + if isinstance(l_type, str) and l_type == "scalar": + assert res_linop.shape == linop2.shape + elif isinstance(r_type, str) and r_type == "scalar": + assert res_linop.shape == linop1.shape + else: + assert res_linop.shape[0] == linop1.shape[0] + assert res_linop.shape[1] == linop2.shape[1] From 63f614f6f1ede75057e43d29c07c9cad8675be83 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Wed, 11 Aug 2021 17:31:50 +0200 Subject: [PATCH 41/74] Fix wrong is_scalar = False in Scaling --- src/probnum/linops/_scaling.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/probnum/linops/_scaling.py b/src/probnum/linops/_scaling.py index 223dbf9bd..5d1b68ac4 100644 --- a/src/probnum/linops/_scaling.py +++ b/src/probnum/linops/_scaling.py @@ -124,6 +124,8 @@ def __init__( trace = lambda: self.shape[0] * self._scalar elif np.ndim(factors) == 1: + if np.all(factors == factors[0]): + self._scalar = factors[0] # Anisotropic scaling self._factors = np.asarray(factors, dtype=dtype) self._factors.setflags(write=False) From b0c0ba5e61762aecf09014a1ef6559dccf61b9b7 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Wed, 11 Aug 2021 17:32:41 +0200 Subject: [PATCH 42/74] Check isotropic scaling --- tests/test_linops/test_arithmetics.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py index 03a8c45ca..d1d377bdc 100644 --- a/tests/test_linops/test_arithmetics.py +++ b/tests/test_linops/test_arithmetics.py @@ -35,9 +35,10 @@ def _aslist(arg): + """Converts anything to a list. Non-iterables become single-element lists.""" try: return list(arg) - except TypeError: + except TypeError: # excepts TypeError: '' object is not iterable return [arg] @@ -52,7 +53,7 @@ def get_linop(linop_type): elif linop_type is Scaling: return ( Scaling(factors=np.random.rand(4)), - Scaling(factors=3.14 * np.ones(4)), + Scaling(factors=3.14), ) elif linop_type is Matrix: return Matrix(np.random.rand(4, 4)) From 26600276fdc20c5d89087ef8178a634020ed8d90 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 12 Aug 2021 08:26:37 +0200 Subject: [PATCH 43/74] Fix tests --- tests/test_linops/test_arithmetics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py index d1d377bdc..d5740941a 100644 --- a/tests/test_linops/test_arithmetics.py +++ b/tests/test_linops/test_arithmetics.py @@ -53,7 +53,7 @@ def get_linop(linop_type): elif linop_type is Scaling: return ( Scaling(factors=np.random.rand(4)), - Scaling(factors=3.14), + Scaling(factors=3.14, shape=(4, 4)), ) elif linop_type is Matrix: return Matrix(np.random.rand(4, 4)) From 8d076ca9a2665f341173f9b385453f0206d56d8f Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 12 Aug 2021 09:23:56 +0200 Subject: [PATCH 44/74] Check failure cases as well --- src/probnum/linops/_arithmetic.py | 46 ++++++++++++++--- src/probnum/linops/_arithmetic_fallbacks.py | 2 +- tests/test_linops/test_arithmetics.py | 55 +++++++++++---------- 3 files changed, 68 insertions(+), 35 deletions(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index 23f26bc99..01a0e6097 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -140,12 +140,22 @@ def _matmul_op_scaled(anylinop, scaled): def _matmul_scaling_kronecker(scaling: Scaling, kronecker: Kronecker) -> Kronecker: + if scaling.shape[1] != kronecker.shape[0]: + raise ValueError( + f"matmul received invalid shapes {scaling.shape} @ {kronecker.shape}" + ) + if scaling.is_isotropic: return scaling.scalar * kronecker return NotImplemented def _matmul_kronecker_scaling(kronecker: Kronecker, scaling: Scaling) -> Kronecker: + if kronecker.shape[1] != scaling.shape[0]: + raise ValueError( + f"matmul received invalid shapes {kronecker.shape} @ {scaling.shape}" + ) + if scaling.is_isotropic: return kronecker * scaling.scalar return NotImplemented @@ -177,6 +187,12 @@ def _mul_kronecker_scalar(kronecker: Kronecker, scalar: ScalarArgType) -> Kronec def _matmul_scaling_idkronecker( scaling: Scaling, idkronecker: IdentityKronecker ) -> IdentityKronecker: + + if scaling.shape[1] != idkronecker.shape[0]: + raise ValueError( + f"matmul received invalid shapes {scaling.shape} @ {idkronecker.shape}" + ) + if scaling.is_isotropic: return scaling.scalar * idkronecker return NotImplemented @@ -185,6 +201,12 @@ def _matmul_scaling_idkronecker( def _matmul_idkronecker_scaling( idkronecker: IdentityKronecker, scaling: Scaling ) -> IdentityKronecker: + + if idkronecker.shape[1] != scaling.shape[0]: + raise ValueError( + f"matmul received invalid shapes {idkronecker.shape} @ {scaling.shape}" + ) + if scaling.is_isotropic: return idkronecker * scaling.scalar return NotImplemented @@ -224,16 +246,10 @@ def _mul_idkronecker_scalar( # Matrix def _matmul_scaling_matrix(scaling: Scaling, matrix: Matrix) -> Matrix: - if scaling.shape[1] != matrix.shape[0]: - return NotImplemented - return Matrix(A=np.multiply(scaling.factors[:, np.newaxis], matrix.A)) def _matmul_matrix_scaling(matrix: Matrix, scaling: Scaling) -> Matrix: - if matrix.shape[1] != scaling.shape[0]: - return NotImplemented - return Matrix(A=np.multiply(matrix.A, scaling.factors)) @@ -284,9 +300,23 @@ def _matmul_inverse_matrix(inverse: _InverseLinearOperator, mat: Matrix) -> Matr # Identity +def _matmul_id_any(idty: Identity, anyop: LinearOperator) -> LinearOperator: + if idty.shape[1] != anyop.shape[0]: + raise ValueError(f"matmul received invalid shapes {idty.shape} @ {anyop.shape}") + + return anyop + + +def _matmul_any_id(anyop: LinearOperator, idty: Identity) -> LinearOperator: + if anyop.shape[1] != idty.shape[0]: + raise ValueError(f"matmul received invalid shapes {anyop.shape} @ {idty.shape}") + + return anyop + + for op_type in _AnyLinOp: - _matmul_fns[(Identity, op_type)] = lambda idty, other: other - _matmul_fns[(op_type, Identity)] = lambda other, idty: other + _matmul_fns[(Identity, op_type)] = _matmul_id_any + _matmul_fns[(op_type, Identity)] = _matmul_any_id # Selection / Embedding diff --git a/src/probnum/linops/_arithmetic_fallbacks.py b/src/probnum/linops/_arithmetic_fallbacks.py index 1f06f40b9..5159def81 100644 --- a/src/probnum/linops/_arithmetic_fallbacks.py +++ b/src/probnum/linops/_arithmetic_fallbacks.py @@ -166,7 +166,7 @@ def __init__(self, *factors: LinearOperator): ): raise ValueError( f"Shape mismatch: Cannot multiply linear operators with shapes " - f"{' x '.join(factor.shape for factor in factors)}." + f"{' x '.join(str(factor.shape) for factor in factors)}." ) self._factors = ProductLinearOperator._expand_prod_ops(*factors) diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py index d5740941a..025c25b43 100644 --- a/tests/test_linops/test_arithmetics.py +++ b/tests/test_linops/test_arithmetics.py @@ -1,7 +1,6 @@ """Tests for linear operator arithmetics.""" import itertools -import operator import numpy as np import pytest @@ -45,18 +44,26 @@ def _aslist(arg): def get_linop(linop_type): # pylint: disable=too-many-return-statements if linop_type is Kronecker: - return Kronecker(np.random.rand(2, 2), np.random.rand(2, 2)) + return ( + Kronecker(np.random.rand(2, 2), np.random.rand(2, 2)), + Kronecker(np.random.rand(4, 3), np.random.rand(2, 3)), + ) elif linop_type is IdentityKronecker: - return IdentityKronecker(2, np.random.rand(2, 2)) + return ( + IdentityKronecker(2, np.random.rand(2, 2)), + IdentityKronecker(3, np.random.rand(3, 4)), + ) elif linop_type is Zero or linop_type is Identity: - return linop_type(shape=(4, 4)) + return (linop_type(shape=(4, 4)), linop_type(shape=(3, 3))) elif linop_type is Scaling: return ( Scaling(factors=np.random.rand(4)), Scaling(factors=3.14, shape=(4, 4)), + Scaling(factors=np.random.rand(6), shape=(6, 6)), + Scaling(factors=3.14, shape=(3, 3)), ) elif linop_type is Matrix: - return Matrix(np.random.rand(4, 4)) + return (Matrix(np.random.rand(4, 4)), Matrix(np.random.rand(6, 3))) elif linop_type is _InverseLinearOperator: _posdef_randmat = random_spd_matrix(rng=np.random.default_rng(123), dim=4) return Matrix(_posdef_randmat).inv() @@ -96,19 +103,9 @@ def get_linop(linop_type): raise TypeError(f"Don't know what to do with type {linop_type}.") -@pytest.mark.parametrize( - "op", [operator.matmul, operator.mul, operator.add, operator.sub] -) -def test_arithmetics(op): - - registry = { - operator.matmul: _matmul_fns, - operator.mul: _mul_fns, - operator.add: _add_fns, - operator.sub: _sub_fns, - }[op] +def test_matmul(): - for (l_type, r_type) in registry.keys(): + for (l_type, r_type) in _matmul_fns.keys(): if ( l_type is Selection or l_type is Embedding @@ -123,13 +120,19 @@ def test_arithmetics(op): for (linop1, linop2) in itertools.product(_aslist(linops1), _aslist(linops2)): - res_linop = op(linop1, linop2) - assert res_linop.ndim == 2 - - if isinstance(l_type, str) and l_type == "scalar": - assert res_linop.shape == linop2.shape - elif isinstance(r_type, str) and r_type == "scalar": - assert res_linop.shape == linop1.shape + if linop1.shape[1] != linop2.shape[0]: + with pytest.raises(ValueError): + res_linop = linop1 @ linop2 + print(res_linop) + print(linop1, linop2) else: - assert res_linop.shape[0] == linop1.shape[0] - assert res_linop.shape[1] == linop2.shape[1] + res_linop = linop1 @ linop2 + assert res_linop.ndim == 2 + + if isinstance(l_type, str) and l_type == "scalar": + assert res_linop.shape == linop2.shape + elif isinstance(r_type, str) and r_type == "scalar": + assert res_linop.shape == linop1.shape + else: + assert res_linop.shape[0] == linop1.shape[0] + assert res_linop.shape[1] == linop2.shape[1] From 36975523b58fcbd662c87d47a27abdabc221aa3c Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 12 Aug 2021 09:42:58 +0200 Subject: [PATCH 45/74] Improve kronecker matmul --- src/probnum/linops/_kronecker.py | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/src/probnum/linops/_kronecker.py b/src/probnum/linops/_kronecker.py index 748f774fa..2157d699b 100644 --- a/src/probnum/linops/_kronecker.py +++ b/src/probnum/linops/_kronecker.py @@ -163,13 +163,15 @@ def _cond(self, p) -> np.inexact: return np.linalg.cond(self.todense(cache=False), p=p) def _matmul_kronecker(self, other: "Kronecker") -> "Kronecker": + if self.shape[1] != other.shape[0]: + raise ValueError( + f"matmul received invalid shapes {self.shape} @ {other.shape}" + ) + _A, _B = self.A, self.B _C, _D = other.A, other.B if not (_A.shape[1] == _C.shape[0] and _B.shape[1] == _D.shape[0]): - raise ValueError( - f"Matmul shape mismatch {_A.shape} x {_C.shape} " - f"or {_B.shape} x {_D.shape}" - ) + return NotImplemented # Using (A (x) B) @ (C (x) D) = (A @ C) (x) (B @ D) return Kronecker(A=_A @ _C, B=_B @ _D) @@ -511,14 +513,16 @@ def num_blocks(self): return self._num_blocks def _matmul_idkronecker(self, other: "IdentityKronecker") -> "IdentityKronecker": - _A, _B = self.A, self.B - _C, _D = other.A, other.B - if not (_A.shape == _C.shape and _B.shape[1] == _D.shape[0]): + if self.shape[1] != other.shape[0]: raise ValueError( - f"Matmul shape mismatch {_A.shape} x {_C.shape} " - f"or {_B.shape} x {_D.shape}" + f"matmul received invalid shapes {self.shape} @ {other.shape}" ) + _A, _B = self.A, self.B + _C, _D = other.A, other.B + if not (_A.shape[1] == _C.shape[0] and _B.shape[1] == _D.shape[0]): + return NotImplemented + # Using (A (x) B) @ (C (x) D) = (A @ C) (x) (B @ D) return IdentityKronecker(num_blocks=self._num_blocks, B=_B @ _D) From 56990d64d295c0fe26bd484890621970531bdd15 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 12 Aug 2021 09:43:03 +0200 Subject: [PATCH 46/74] Adapt tests --- tests/test_linops/test_arithmetics.py | 99 ++++++++++++++++++++++++++- 1 file changed, 97 insertions(+), 2 deletions(-) diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py index 025c25b43..68b9b98af 100644 --- a/tests/test_linops/test_arithmetics.py +++ b/tests/test_linops/test_arithmetics.py @@ -123,8 +123,6 @@ def test_matmul(): if linop1.shape[1] != linop2.shape[0]: with pytest.raises(ValueError): res_linop = linop1 @ linop2 - print(res_linop) - print(linop1, linop2) else: res_linop = linop1 @ linop2 assert res_linop.ndim == 2 @@ -136,3 +134,100 @@ def test_matmul(): else: assert res_linop.shape[0] == linop1.shape[0] assert res_linop.shape[1] == linop2.shape[1] + + +def test_mul(): + + for (l_type, r_type) in _mul_fns.keys(): + if ( + l_type is Selection + or l_type is Embedding + or r_type is Selection + or r_type is Embedding + ): + # Checked seperatly + continue + + linops1 = get_linop(l_type) + linops2 = get_linop(r_type) + + for (linop1, linop2) in itertools.product(_aslist(linops1), _aslist(linops2)): + + if isinstance(l_type, str) and l_type == "scalar": + res_linop = linop1 * linop2 + assert res_linop.shape == linop2.shape + elif isinstance(r_type, str) and r_type == "scalar": + res_linop = linop1 * linop2 + assert res_linop.shape == linop1.shape + else: + if linop1.shape != linop2.shape: + with pytest.raises(ValueError): + res_linop = linop1 * linop2 + else: + res_linop = linop1 * linop2 + assert res_linop.shape == linop1.shape == linop2.shape + + +def test_add(): + + for (l_type, r_type) in _add_fns.keys(): + if ( + l_type is Selection + or l_type is Embedding + or r_type is Selection + or r_type is Embedding + ): + # Checked seperatly + continue + + linops1 = get_linop(l_type) + linops2 = get_linop(r_type) + + for (linop1, linop2) in itertools.product(_aslist(linops1), _aslist(linops2)): + + if linop1.shape != linop2.shape: + with pytest.raises(ValueError): + res_linop = linop1 + linop2 + else: + res_linop = linop1 + linop2 + assert res_linop.shape == linop1.shape == linop2.shape + + +def test_sub(): + + for (l_type, r_type) in _sub_fns.keys(): + if ( + l_type is Selection + or l_type is Embedding + or r_type is Selection + or r_type is Embedding + ): + # Checked seperatly + continue + + linops1 = get_linop(l_type) + linops2 = get_linop(r_type) + + for (linop1, linop2) in itertools.product(_aslist(linops1), _aslist(linops2)): + + if linop1.shape != linop2.shape: + with pytest.raises(ValueError): + res_linop = linop1 - linop2 + else: + res_linop = linop1 - linop2 + assert res_linop.shape == linop1.shape == linop2.shape + + +def test_kronecker_matmul(): + # Checks the case in which the shapes of the Kronecker-structured matrices + # are valid in itself but the respective Kronecker factors (k1.A @ k2.A and/or + # k1.B @ k2.B) have invalid shapes for matmul. + k1 = Kronecker(np.random.rand(4, 2), np.random.rand(2, 3)) # (8, 6) + k2 = Kronecker(np.random.rand(3, 2), np.random.rand(2, 3)) # (6, 6) + + # Even though the shapes fit, and Kronecker @ Kronecker = Kronecker .... + assert k1.shape[1] == k2.shape[0] + + # The result does not have a Kronecker structure + res = k1 @ k2 + assert not isinstance(res, Kronecker) From 41fe141c6015b12ddfbe664376fdeaff83b472d4 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 12 Aug 2021 09:46:40 +0200 Subject: [PATCH 47/74] Add test for selection and embedding --- tests/test_linops/test_arithmetics.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py index 68b9b98af..3b7bf98bf 100644 --- a/tests/test_linops/test_arithmetics.py +++ b/tests/test_linops/test_arithmetics.py @@ -231,3 +231,12 @@ def test_kronecker_matmul(): # The result does not have a Kronecker structure res = k1 @ k2 assert not isinstance(res, Kronecker) + + +def test_selection_embedding(): + sel = get_linop(Selection) + emb = get_linop(Embedding) + + product = sel @ emb + assert product.shape[0] == sel.shape[0] + assert product.shape[1] == emb.shape[1] From c8f83bd599a0cdbe530989c8366c42e4b277b26c Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 12 Aug 2021 09:58:17 +0200 Subject: [PATCH 48/74] Use lambda expr in scalar matrix mult --- src/probnum/linops/_arithmetic.py | 18 ++---------------- 1 file changed, 2 insertions(+), 16 deletions(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index 01a0e6097..ec1124493 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -253,22 +253,8 @@ def _matmul_matrix_scaling(matrix: Matrix, scaling: Scaling) -> Matrix: return Matrix(A=np.multiply(matrix.A, scaling.factors)) -def _mul_matrix_scalar(mat, scalar) -> Union[NotImplementedType, Matrix]: - if np.isscalar(scalar): - return Matrix(A=scalar * mat.A) - - return NotImplemented - - -def _mul_scalar_matrix(scalar, mat) -> Union[NotImplementedType, Matrix]: - if np.isscalar(scalar): - return Matrix(A=scalar * mat.A) - - return NotImplemented - - -_mul_fns[(Matrix, "scalar")] = _mul_matrix_scalar -_mul_fns[("scalar", Matrix)] = _mul_scalar_matrix +_mul_fns[(Matrix, "scalar")] = lambda mat, scal: Matrix(A=scal * mat.A) +_mul_fns[("scalar", Matrix)] = lambda scal, mat: Matrix(A=scal * mat.A) _matmul_fns[(Matrix, Matrix)] = Matrix._matmul_matrix _matmul_fns[(Scaling, Matrix)] = _matmul_scaling_matrix _matmul_fns[(Matrix, Scaling)] = _matmul_matrix_scaling From e079aea21565c847f6abd40f219192fc95a1f528 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 12 Aug 2021 10:12:04 +0200 Subject: [PATCH 49/74] Extended config option test --- src/probnum/linops/_arithmetic.py | 18 ++++++++++++------ tests/test_linops/test_arithmetics.py | 22 ++++++++++++++++++++++ 2 files changed, 34 insertions(+), 6 deletions(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index ec1124493..8b881b625 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -269,20 +269,26 @@ def _matmul_matrix_scaling(matrix: Matrix, scaling: Scaling) -> Matrix: _sub_fns[(Matrix, Matrix)] = lambda mat1, mat2: Matrix(mat1.A - mat2.A) -def _matmul_matrix_inverse(mat: Matrix, inverse: _InverseLinearOperator) -> Matrix: +def _matmul_matrix_wrapped( + mat: Matrix, wrapped: Union[_InverseLinearOperator, TransposedLinearOperator] +) -> Union[Matrix, NotImplementedType]: if config.lazy_matrix_matrix_matmul: - return Matrix(mat.A @ inverse) + return Matrix(mat.A @ wrapped) return NotImplemented -def _matmul_inverse_matrix(inverse: _InverseLinearOperator, mat: Matrix) -> Matrix: +def _matmul_wrapped_matrix( + wrapped: Union[_InverseLinearOperator, TransposedLinearOperator], mat: Matrix +) -> Union[Matrix, NotImplementedType]: if config.lazy_matrix_matrix_matmul: - return Matrix(mat.A @ inverse) + return Matrix(wrapped @ mat.A) return NotImplemented -_matmul_fns[(Matrix, _InverseLinearOperator)] = _matmul_matrix_inverse -_matmul_fns[(_InverseLinearOperator, Matrix)] = _matmul_inverse_matrix +_matmul_fns[(Matrix, _InverseLinearOperator)] = _matmul_matrix_wrapped +_matmul_fns[(_InverseLinearOperator, Matrix)] = _matmul_wrapped_matrix +_matmul_fns[(Matrix, TransposedLinearOperator)] = _matmul_matrix_wrapped +_matmul_fns[(TransposedLinearOperator, Matrix)] = _matmul_wrapped_matrix # Identity diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py index 3b7bf98bf..d80acb41a 100644 --- a/tests/test_linops/test_arithmetics.py +++ b/tests/test_linops/test_arithmetics.py @@ -5,6 +5,7 @@ import numpy as np import pytest +from probnum import config from probnum.linops._arithmetic import _add_fns, _matmul_fns, _mul_fns, _sub_fns from probnum.linops._arithmetic_fallbacks import ( NegatedLinearOperator, @@ -240,3 +241,24 @@ def test_selection_embedding(): product = sel @ emb assert product.shape[0] == sel.shape[0] assert product.shape[1] == emb.shape[1] + + +def test_lazy_matrix_matrix_matmul_option(): + mat1 = get_linop(Matrix)[0] + mat2 = get_linop(Matrix)[0] + inv = get_linop(_InverseLinearOperator) + transposed = get_linop(TransposedLinearOperator) + + with config(lazy_matrix_matrix_matmul=False): + assert isinstance(mat1 @ mat2, ProductLinearOperator) + assert isinstance(mat1 @ inv, ProductLinearOperator) + assert isinstance(inv @ mat2, ProductLinearOperator) + assert isinstance(mat1 @ transposed, ProductLinearOperator) + assert isinstance(transposed @ mat2, ProductLinearOperator) + + with config(lazy_matrix_matrix_matmul=True): + assert isinstance(mat1 @ mat2, Matrix) + assert isinstance(mat1 @ inv, Matrix) + assert isinstance(inv @ mat2, Matrix) + assert isinstance(mat1 @ transposed, Matrix) + assert isinstance(transposed @ mat2, Matrix) From c5e12723e97d46fa77c805fc3a469f9f981d7d4b Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 12 Aug 2021 10:50:31 +0200 Subject: [PATCH 50/74] Test arithmetic fallbacks --- src/probnum/linops/_arithmetic_fallbacks.py | 20 +-------- .../test_linops/test_arithmetics_fallbacks.py | 45 +++++++++++++++++++ 2 files changed, 46 insertions(+), 19 deletions(-) create mode 100644 tests/test_linops/test_arithmetics_fallbacks.py diff --git a/src/probnum/linops/_arithmetic_fallbacks.py b/src/probnum/linops/_arithmetic_fallbacks.py index 5159def81..3c60390e3 100644 --- a/src/probnum/linops/_arithmetic_fallbacks.py +++ b/src/probnum/linops/_arithmetic_fallbacks.py @@ -47,12 +47,6 @@ def _inv(self) -> "ScaledLinearOperator": return ScaledLinearOperator(self._linop.inv(), 1.0 / self._scalar) - def __mul__(self, other: BinaryOperandType) -> "LinearOperator": - if np.isscalar(other): - return ScaledLinearOperator(linop=self._linop, scalar=self._scalar * other) - - return super().__mul__(other) - def __repr__(self) -> str: return f"{self._scalar} * {self._linop}" @@ -72,11 +66,6 @@ class SumLinearOperator(LinearOperator): """Sum of two linear operators.""" def __init__(self, *summands: LinearOperator): - if not all(isinstance(summand, LinearOperator) for summand in summands): - raise TypeError("All summands must be `LinearOperator`s") - - if len(summands) < 2: - raise ValueError("There must be at least two summands") if not all(summand.shape == summands[0].shape for summand in summands): raise ValueError("All summands must have the same shape") @@ -136,9 +125,7 @@ def _mul_fallback( ) -> Union[LinearOperator, NotImplementedType]: res = NotImplemented - if isinstance(op1, LinearOperator) and isinstance(op2, LinearOperator): - pass # TODO: Implement generic Hadamard product - elif isinstance(op1, LinearOperator): + if isinstance(op1, LinearOperator): if np.ndim(op2) == 0: res = ScaledLinearOperator(op1, op2) elif isinstance(op2, LinearOperator): @@ -154,11 +141,6 @@ class ProductLinearOperator(LinearOperator): """(Operator) Product of two linear operators.""" def __init__(self, *factors: LinearOperator): - if not all(isinstance(factor, LinearOperator) for factor in factors): - raise TypeError("All factors must be `LinearOperator`s") - - if len(factors) < 2: - raise ValueError("There must be at least two factors") if not all( lfactor.shape[1] == rfactor.shape[0] diff --git a/tests/test_linops/test_arithmetics_fallbacks.py b/tests/test_linops/test_arithmetics_fallbacks.py new file mode 100644 index 000000000..ad810f1d5 --- /dev/null +++ b/tests/test_linops/test_arithmetics_fallbacks.py @@ -0,0 +1,45 @@ +"""Tests for linear operator arithmetics fallbacks.""" + +import numpy as np +import pytest + +from probnum.linops._arithmetic_fallbacks import ( + NegatedLinearOperator, + ProductLinearOperator, + ScaledLinearOperator, + SumLinearOperator, +) +from probnum.linops._linear_operator import Matrix +from probnum.problems.zoo.linalg import random_spd_matrix + + +@pytest.fixture +def rng(): + return np.random.default_rng(123) + + +@pytest.fixture +def scalar(): + return 3.14 + + +@pytest.fixture +def rand_spd_mat(rng): + return Matrix(random_spd_matrix(rng, dim=4)) + + +def test_scaled_linop(rand_spd_mat, scalar): + with pytest.raises(TypeError): + ScaledLinearOperator(np.random.rand(4, 4), scalar=scalar) + with pytest.raises(TypeError): + ScaledLinearOperator(rand_spd_mat, scalar=np.ones(4)) + + scaled1 = ScaledLinearOperator(rand_spd_mat, scalar=0.0) + scaled2 = ScaledLinearOperator(rand_spd_mat, scalar=scalar) + + with pytest.raises(np.linalg.LinAlgError): + scaled1.inv() + + assert np.allclose( + scaled2.inv().todense(), (1.0 / scalar) * scaled2._linop.inv().todense() + ) From 34ac9e6a76bdea6a3cd0ba81e74e3892f2c5815a Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 12 Aug 2021 10:57:21 +0200 Subject: [PATCH 51/74] Now build. From f362ba6723770f8aad04ee92b78886151393d373 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 12 Aug 2021 12:47:16 +0200 Subject: [PATCH 52/74] Fix pylint --- tests/test_linops/test_arithmetics.py | 1 + tests/test_linops/test_arithmetics_fallbacks.py | 5 +---- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py index d80acb41a..687dbc897 100644 --- a/tests/test_linops/test_arithmetics.py +++ b/tests/test_linops/test_arithmetics.py @@ -1,4 +1,5 @@ """Tests for linear operator arithmetics.""" +# pylint: disable=consider-iterating-dictionary import itertools diff --git a/tests/test_linops/test_arithmetics_fallbacks.py b/tests/test_linops/test_arithmetics_fallbacks.py index ad810f1d5..f2d9138ab 100644 --- a/tests/test_linops/test_arithmetics_fallbacks.py +++ b/tests/test_linops/test_arithmetics_fallbacks.py @@ -3,11 +3,8 @@ import numpy as np import pytest -from probnum.linops._arithmetic_fallbacks import ( - NegatedLinearOperator, - ProductLinearOperator, +from probnum.linops._arithmetic_fallbacks import ( # NegatedLinearOperator,; ProductLinearOperator,; SumLinearOperator, ScaledLinearOperator, - SumLinearOperator, ) from probnum.linops._linear_operator import Matrix from probnum.problems.zoo.linalg import random_spd_matrix From b9a27df6fa7b3242b9a8707ccdae96193436009b Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 12 Aug 2021 12:51:03 +0200 Subject: [PATCH 53/74] Cover Kronecker distributive --- tests/test_linops/test_arithmetics.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py index 687dbc897..dec3d739c 100644 --- a/tests/test_linops/test_arithmetics.py +++ b/tests/test_linops/test_arithmetics.py @@ -46,9 +46,13 @@ def _aslist(arg): def get_linop(linop_type): # pylint: disable=too-many-return-statements if linop_type is Kronecker: + _A1 = np.random.rand(2, 2) + _B1 = np.random.rand(2, 3) return ( - Kronecker(np.random.rand(2, 2), np.random.rand(2, 2)), - Kronecker(np.random.rand(4, 3), np.random.rand(2, 3)), + Kronecker(_A1, np.random.rand(2, 2)), + Kronecker(np.random.rand(4, 3), _B1), + Kronecker(_A1, np.random.rand(2, 2)), + Kronecker(np.random.rand(2, 2), _B1), ) elif linop_type is IdentityKronecker: return ( From 0447cdb3f6bd99194a853d6fb76f0c7cfccfe497 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 12 Aug 2021 13:14:20 +0200 Subject: [PATCH 54/74] Add matrix free IWP transition tests --- .../test_markov/test_integrator/test_iwp.py | 50 +++++++++++++++++++ 1 file changed, 50 insertions(+) diff --git a/tests/test_randprocs/test_markov/test_integrator/test_iwp.py b/tests/test_randprocs/test_markov/test_integrator/test_iwp.py index 542796a24..7927c959c 100644 --- a/tests/test_randprocs/test_markov/test_integrator/test_iwp.py +++ b/tests/test_randprocs/test_markov/test_integrator/test_iwp.py @@ -277,3 +277,53 @@ def test_forward_realization_values( ) np.testing.assert_allclose(ah_22_ibm @ real, rv.mean) np.testing.assert_allclose(diffusion * qh_22_ibm, rv.cov) + + +class TestIntegratedWienerTransitionValuesLinOps: + + # Replacement for an __init__ in the pytest language. See: + # https://stackoverflow.com/questions/21430900/py-test-skips-test-class-if-constructor-is-defined + @pytest.fixture(autouse=True) + def _setup( + self, + ): + wiener_process_dimension = 1 # make tests compatible with some_normal_rv1, etc. + with config(matrix_free_linalg=True): + self.transition = randprocs.markov.integrator.IntegratedWienerTransition( + num_derivatives=2, + wiener_process_dimension=wiener_process_dimension, + forward_implementation="classic", + backward_implementation="classic", + ) + + def test_discretise_values(self, ah_22_ibm, qh_22_ibm, dt): + with config(matrix_free_linalg=True): + discrete_model = self.transition.discretise(dt=dt) + np.testing.assert_allclose( + discrete_model.state_trans_mat.todense(), ah_22_ibm + ) + np.testing.assert_allclose( + discrete_model.proc_noise_cov_mat.todense(), qh_22_ibm + ) + + def test_forward_rv_values(self, normal_rv3x3, diffusion, ah_22_ibm, qh_22_ibm, dt): + with config(matrix_free_linalg=True): + rv, _ = self.transition.forward_rv( + normal_rv3x3, t=0.0, dt=dt, _diffusion=diffusion + ) + np.testing.assert_allclose(ah_22_ibm @ normal_rv3x3.mean, rv[:3].mean) + np.testing.assert_allclose( + ah_22_ibm @ normal_rv3x3.cov @ ah_22_ibm.T + diffusion * qh_22_ibm, + rv.cov.todense(), + ) + + def test_forward_realization_values( + self, normal_rv3x3, diffusion, ah_22_ibm, qh_22_ibm, dt + ): + with config(matrix_free_linalg=True): + real = normal_rv3x3.mean + rv, _ = self.transition.forward_realization( + real, t=0.0, dt=dt, _diffusion=diffusion + ) + np.testing.assert_allclose(ah_22_ibm @ real, rv.mean) + np.testing.assert_allclose(diffusion * qh_22_ibm, rv.cov.todense()) From 10b4877ed16391a5bc0ab4a8cc228ae23a7d7d1f Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 12 Aug 2021 14:05:44 +0200 Subject: [PATCH 55/74] Fix embedding matmul and add tests --- src/probnum/linops/_linear_operator.py | 5 ++-- .../selectionembedding_cases.py | 28 +++++++++++++++++++ 2 files changed, 30 insertions(+), 3 deletions(-) create mode 100644 tests/test_linops/test_linops_cases/selectionembedding_cases.py diff --git a/src/probnum/linops/_linear_operator.py b/src/probnum/linops/_linear_operator.py index c83b7b896..696962311 100644 --- a/src/probnum/linops/_linear_operator.py +++ b/src/probnum/linops/_linear_operator.py @@ -1199,14 +1199,13 @@ def __init__( ) def _todense(self): - res = np.eye(self.shape[1], self.shape[1]) - return _embedding_matmul(self, res) + return self.T.todense().T def _embedding_matmul(embedding, M): res_shape = np.array(M.shape) res_shape[-2] = embedding.shape[0] res = np.full(shape=tuple(res_shape), fill_value=embedding._fill_value) - take_from_M = M[..., np.array(embedding._put_indices), :] + take_from_M = M[..., np.array(embedding._take_indices), :] res[..., np.array(embedding._put_indices), :] = take_from_M return res diff --git a/tests/test_linops/test_linops_cases/selectionembedding_cases.py b/tests/test_linops/test_linops_cases/selectionembedding_cases.py new file mode 100644 index 000000000..a5e626a45 --- /dev/null +++ b/tests/test_linops/test_linops_cases/selectionembedding_cases.py @@ -0,0 +1,28 @@ +from typing import Tuple + +import numpy as np +import pytest + +import probnum as pn + + +def case_embedding(): + return ( + pn.linops.Embedding( + take_indices=(0, 1, 2), put_indices=(1, 0, 3), shape=(4, 3) + ), + np.array( + [ + [0.0, 1.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 0.0, 0.0], + [0.0, 0.0, 1.0], + ] + ), + ) + + +def case_selection(): + return pn.linops.Selection(indices=(1, 0, 3), shape=(3, 4)), np.array( + [[0.0, 1.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0]] + ) From b7f1a7f6eb0ff3c57739c793d88f50d3c8ece2e4 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 12 Aug 2021 14:05:57 +0200 Subject: [PATCH 56/74] Test scaling inequality --- tests/test_linops/test_arithmetics.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py index dec3d739c..8562d0b55 100644 --- a/tests/test_linops/test_arithmetics.py +++ b/tests/test_linops/test_arithmetics.py @@ -267,3 +267,18 @@ def test_lazy_matrix_matrix_matmul_option(): assert isinstance(inv @ mat2, Matrix) assert isinstance(mat1 @ transposed, Matrix) assert isinstance(transposed @ mat2, Matrix) + + +def test_equality(): + scalings = get_linop(Scaling) + int_scaling = Scaling(2, shape=(4, 4)) + for s1, s2 in itertools.product(scalings, _aslist(scalings) + [int_scaling]): + if ( + s1.shape == s2.shape + and s1.dtype == s2.dtype + and np.all(s1.todense() == s2.todense()) + ): + assert s1 == s2 + + else: + assert s1 != s2 From 373dd2ffaecd7d4871eb8c1bb295cfccc951d53a Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 12 Aug 2021 14:19:25 +0200 Subject: [PATCH 57/74] Extend tests --- src/probnum/linops/_arithmetic_fallbacks.py | 3 --- tests/test_linops/test_arithmetics.py | 13 ++++++++++--- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/probnum/linops/_arithmetic_fallbacks.py b/src/probnum/linops/_arithmetic_fallbacks.py index 3c60390e3..1ecdb289c 100644 --- a/src/probnum/linops/_arithmetic_fallbacks.py +++ b/src/probnum/linops/_arithmetic_fallbacks.py @@ -131,9 +131,6 @@ def _mul_fallback( elif isinstance(op2, LinearOperator): if np.ndim(op1) == 0: res = ScaledLinearOperator(op2, op1) - else: - raise TypeError("At least one of the two operands must be a `LinearOperator`.") - return res diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py index 8562d0b55..d5c4f44a4 100644 --- a/tests/test_linops/test_arithmetics.py +++ b/tests/test_linops/test_arithmetics.py @@ -242,10 +242,17 @@ def test_kronecker_matmul(): def test_selection_embedding(): sel = get_linop(Selection) emb = get_linop(Embedding) + emb2 = Embedding( + take_indices=emb._take_indices, put_indices=emb._put_indices, shape=(5, 3) + ) - product = sel @ emb - assert product.shape[0] == sel.shape[0] - assert product.shape[1] == emb.shape[1] + product1 = sel @ emb + assert product1.shape[0] == sel.shape[0] + assert product1.shape[1] == emb.shape[1] + + product2 = sel @ emb2 + assert product2.shape[0] == sel.shape[0] + assert product2.shape[1] == emb2.shape[1] def test_lazy_matrix_matrix_matmul_option(): From 78d36983daf77d9b462e7e99d1673401f1f4a2f6 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 12 Aug 2021 14:20:43 +0200 Subject: [PATCH 58/74] Fix pylint --- .../test_linops/test_linops_cases/selectionembedding_cases.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/tests/test_linops/test_linops_cases/selectionembedding_cases.py b/tests/test_linops/test_linops_cases/selectionembedding_cases.py index a5e626a45..ac5a54103 100644 --- a/tests/test_linops/test_linops_cases/selectionembedding_cases.py +++ b/tests/test_linops/test_linops_cases/selectionembedding_cases.py @@ -1,7 +1,4 @@ -from typing import Tuple - import numpy as np -import pytest import probnum as pn From 5963565ed24b7e7bd3b6bc6f0ca661812c99c654 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 12 Aug 2021 14:40:50 +0200 Subject: [PATCH 59/74] Stabilize tests --- src/probnum/linops/_kronecker.py | 5 ----- tests/test_linops/test_arithmetics.py | 4 ++-- 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/src/probnum/linops/_kronecker.py b/src/probnum/linops/_kronecker.py index 2157d699b..99d30f162 100644 --- a/src/probnum/linops/_kronecker.py +++ b/src/probnum/linops/_kronecker.py @@ -192,11 +192,6 @@ def _sub_kronecker( self, other: "Kronecker" ) -> Union[NotImplementedType, "Kronecker"]: - if self.A == other.A and self.B == other.B: - return Kronecker( - A=Scaling(0.0, shape=self.A.shape), B=Scaling(0.0, shape=self.B.shape) - ) - if self.A == other.A: return Kronecker(A=self.A, B=self.B - other.B) diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py index d5c4f44a4..bf8e8b2d1 100644 --- a/tests/test_linops/test_arithmetics.py +++ b/tests/test_linops/test_arithmetics.py @@ -46,8 +46,8 @@ def _aslist(arg): def get_linop(linop_type): # pylint: disable=too-many-return-statements if linop_type is Kronecker: - _A1 = np.random.rand(2, 2) - _B1 = np.random.rand(2, 3) + _A1 = np.ones((2, 2)) + _B1 = np.ones((2, 3)) return ( Kronecker(_A1, np.random.rand(2, 2)), Kronecker(np.random.rand(4, 3), _B1), From 2d8024c244bd9d8ed6c50019c21013b5e1df2978 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 12 Aug 2021 14:50:58 +0200 Subject: [PATCH 60/74] Test IdentityKronecker --- src/probnum/linops/_kronecker.py | 1 - tests/test_linops/test_arithmetics.py | 15 +++++++++++++++ 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/src/probnum/linops/_kronecker.py b/src/probnum/linops/_kronecker.py index 99d30f162..b14e8775a 100644 --- a/src/probnum/linops/_kronecker.py +++ b/src/probnum/linops/_kronecker.py @@ -6,7 +6,6 @@ from probnum.typing import DTypeArgType, NotImplementedType from . import _linear_operator, _utils -from ._scaling import Scaling class Symmetrize(_linear_operator.LinearOperator): diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py index bf8e8b2d1..f6d0b2541 100644 --- a/tests/test_linops/test_arithmetics.py +++ b/tests/test_linops/test_arithmetics.py @@ -239,6 +239,21 @@ def test_kronecker_matmul(): assert not isinstance(res, Kronecker) +def test_idkronecker_matmul(): + # Checks the case in which the shapes of the Kronecker-structured matrices + # are valid in itself but the respective Kronecker factors (k1.A @ k2.A and/or + # k1.B @ k2.B) have invalid shapes for matmul. + k1 = IdentityKronecker(4, np.random.rand(2, 3)) # (8, 12) + k2 = IdentityKronecker(2, np.random.rand(6, 2)) # (12, 4) + + # Even though the shapes fit, and IdentityKronecker @ IdentityKronecker = IdentityKronecker .... + assert k1.shape[1] == k2.shape[0] + + # The result does not have a IdentityKronecker structure + res = k1 @ k2 + assert not isinstance(res, IdentityKronecker) + + def test_selection_embedding(): sel = get_linop(Selection) emb = get_linop(Embedding) From 316533bd5f7b22a9f4fce1310f0b97e7f2d431e9 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 18 Nov 2021 16:23:38 +0100 Subject: [PATCH 61/74] Add restructured text to config docs --- src/probnum/_config.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/probnum/_config.py b/src/probnum/_config.py index e80608f71..70f5fd585 100644 --- a/src/probnum/_config.py +++ b/src/probnum/_config.py @@ -130,16 +130,18 @@ def register(self, key: str, default_value: Any, description: str) -> None: "matrix_free", False, ( - "If True, wherever possible, LinearOperators are used instead " - "of arrays. LinearOperators define a matrix-vector product implicitly without instantiating the full matrix in memory. This makes them memory- and runtime-efficient for linear algebra operations." + "If ``True``, wherever possible, ``LinearOperator``s are used instead " + "of arrays. ``LinearOperator``s define a matrix-vector product implicitly without " + "instantiating the full matrix in memory. This makes them memory- and runtime-efficient " + "for linear algebra operations." ), ), ( "lazy_matrix_matrix_matmul", False, - "If True, the matmul operation applied to two Matrix-LinearOperators will " - "again yield a Matrix-LinearOperator with the computed matrix-matrix product," - " instead of a ProductLinearOperator.", + "If ``True``, the matmul operation applied to two ``Matrix``-``LinearOperator``s will " + "again yield a ``Matrix``-``LinearOperator`` with the computed matrix-matrix product," + " instead of a ``ProductLinearOperator``.", ), ] From f01c8c802b007e2d487b8e2c09e1ed1c832593f0 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 18 Nov 2021 16:28:45 +0100 Subject: [PATCH 62/74] Handle mult with negative scalars --- src/probnum/linops/_arithmetic.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index 8b881b625..3af966ab2 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -162,11 +162,15 @@ def _matmul_kronecker_scaling(kronecker: Kronecker, scaling: Scaling) -> Kroneck def _mul_scalar_kronecker(scalar: ScalarArgType, kronecker: Kronecker) -> Kronecker: + if scalar < 0.0: + return NotImplemented sqrt_scalar = np.sqrt(scalar) return Kronecker(A=sqrt_scalar * kronecker.A, B=sqrt_scalar * kronecker.B) def _mul_kronecker_scalar(kronecker: Kronecker, scalar: ScalarArgType) -> Kronecker: + if scalar < 0.0: + return NotImplemented sqrt_scalar = np.sqrt(scalar) return Kronecker(A=sqrt_scalar * kronecker.A, B=sqrt_scalar * kronecker.B) From dfdc9d0809f007c5835a0fedd0021ed189fdbfc4 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 18 Nov 2021 16:31:24 +0100 Subject: [PATCH 63/74] Replace isscalar by ndim == 0 --- src/probnum/linops/_arithmetic.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index 3af966ab2..54ba02a64 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -399,12 +399,12 @@ def _apply( ] ] = None, ) -> Union[LinearOperator, NotImplementedType]: - if np.isscalar(op1): + if np.ndim(op1) == 0: key1 = "scalar" else: key1 = type(op1) - if np.isscalar(op2): + if np.ndim(op2) == 0: key2 = "scalar" else: key2 = type(op2) From 8495c9cc7dfa6665f648ea63074e183759452104 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 18 Nov 2021 16:32:18 +0100 Subject: [PATCH 64/74] Replace "scalar" by numbers.Number --- src/probnum/linops/_arithmetic.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index 54ba02a64..1173b663d 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -1,4 +1,5 @@ """Linear operator arithmetic.""" +import numbers from typing import Any, Callable, Dict, Optional, Tuple, Union import numpy as np @@ -113,8 +114,8 @@ def _mul_scaling_scalar(scaling: Scaling, scalar: ScalarArgType) -> Scaling: return Scaling(scalar * scaling.factors, shape=scaling.shape) -_mul_fns[("scalar", Scaling)] = _mul_scalar_scaling -_mul_fns[(Scaling, "scalar")] = _mul_scaling_scalar +_mul_fns[(numbers.Number, Scaling)] = _mul_scalar_scaling +_mul_fns[(Scaling, numbers.Number)] = _mul_scaling_scalar _add_fns[(Scaling, Scaling)] = Scaling._add_scaling _sub_fns[(Scaling, Scaling)] = Scaling._sub_scaling _mul_fns[(Scaling, Scaling)] = Scaling._mul_scaling @@ -179,8 +180,8 @@ def _mul_kronecker_scalar(kronecker: Kronecker, scalar: ScalarArgType) -> Kronec _add_fns[(Kronecker, Kronecker)] = Kronecker._add_kronecker _sub_fns[(Kronecker, Kronecker)] = Kronecker._sub_kronecker -_mul_fns[("scalar", Kronecker)] = _mul_scalar_kronecker -_mul_fns[(Kronecker, "scalar")] = _mul_kronecker_scalar +_mul_fns[(numbers.Number, Kronecker)] = _mul_scalar_kronecker +_mul_fns[(Kronecker, numbers.Number)] = _mul_kronecker_scalar _matmul_fns[(Kronecker, Scaling)] = _matmul_kronecker_scaling _matmul_fns[(Scaling, Kronecker)] = _matmul_scaling_kronecker @@ -240,8 +241,8 @@ def _mul_idkronecker_scalar( _add_fns[(IdentityKronecker, IdentityKronecker)] = IdentityKronecker._add_idkronecker _sub_fns[(IdentityKronecker, IdentityKronecker)] = IdentityKronecker._sub_idkronecker -_mul_fns[("scalar", IdentityKronecker)] = _mul_scalar_idkronecker -_mul_fns[(IdentityKronecker, "scalar")] = _mul_idkronecker_scalar +_mul_fns[(numbers.Number, IdentityKronecker)] = _mul_scalar_idkronecker +_mul_fns[(IdentityKronecker, numbers.Number)] = _mul_idkronecker_scalar _matmul_fns[(IdentityKronecker, Scaling)] = _matmul_idkronecker_scaling _matmul_fns[(Scaling, IdentityKronecker)] = _matmul_scaling_idkronecker @@ -257,8 +258,8 @@ def _matmul_matrix_scaling(matrix: Matrix, scaling: Scaling) -> Matrix: return Matrix(A=np.multiply(matrix.A, scaling.factors)) -_mul_fns[(Matrix, "scalar")] = lambda mat, scal: Matrix(A=scal * mat.A) -_mul_fns[("scalar", Matrix)] = lambda scal, mat: Matrix(A=scal * mat.A) +_mul_fns[(Matrix, numbers.Number)] = lambda mat, scal: Matrix(A=scal * mat.A) +_mul_fns[(numbers.Number, Matrix)] = lambda scal, mat: Matrix(A=scal * mat.A) _matmul_fns[(Matrix, Matrix)] = Matrix._matmul_matrix _matmul_fns[(Scaling, Matrix)] = _matmul_scaling_matrix _matmul_fns[(Matrix, Scaling)] = _matmul_matrix_scaling @@ -400,12 +401,12 @@ def _apply( ] = None, ) -> Union[LinearOperator, NotImplementedType]: if np.ndim(op1) == 0: - key1 = "scalar" + key1 = numbers.Number else: key1 = type(op1) if np.ndim(op2) == 0: - key2 = "scalar" + key2 = numbers.Number else: key2 = type(op2) From 580ee8cdbebf00ed987c93d8a773ed296beccecb Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 18 Nov 2021 16:34:08 +0100 Subject: [PATCH 65/74] Nevermind... --- src/probnum/linops/_arithmetic.py | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index 1173b663d..b1803166d 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -1,11 +1,10 @@ """Linear operator arithmetic.""" -import numbers from typing import Any, Callable, Dict, Optional, Tuple, Union import numpy as np import scipy.sparse -from probnum import config +from probnum import config, utils from probnum.typing import NotImplementedType, ScalarArgType, ShapeArgType from ._arithmetic_fallbacks import ( @@ -114,8 +113,8 @@ def _mul_scaling_scalar(scaling: Scaling, scalar: ScalarArgType) -> Scaling: return Scaling(scalar * scaling.factors, shape=scaling.shape) -_mul_fns[(numbers.Number, Scaling)] = _mul_scalar_scaling -_mul_fns[(Scaling, numbers.Number)] = _mul_scaling_scalar +_mul_fns[(np.number, Scaling)] = _mul_scalar_scaling +_mul_fns[(Scaling, np.number)] = _mul_scaling_scalar _add_fns[(Scaling, Scaling)] = Scaling._add_scaling _sub_fns[(Scaling, Scaling)] = Scaling._sub_scaling _mul_fns[(Scaling, Scaling)] = Scaling._mul_scaling @@ -180,8 +179,8 @@ def _mul_kronecker_scalar(kronecker: Kronecker, scalar: ScalarArgType) -> Kronec _add_fns[(Kronecker, Kronecker)] = Kronecker._add_kronecker _sub_fns[(Kronecker, Kronecker)] = Kronecker._sub_kronecker -_mul_fns[(numbers.Number, Kronecker)] = _mul_scalar_kronecker -_mul_fns[(Kronecker, numbers.Number)] = _mul_kronecker_scalar +_mul_fns[(np.number, Kronecker)] = _mul_scalar_kronecker +_mul_fns[(Kronecker, np.number)] = _mul_kronecker_scalar _matmul_fns[(Kronecker, Scaling)] = _matmul_kronecker_scaling _matmul_fns[(Scaling, Kronecker)] = _matmul_scaling_kronecker @@ -241,8 +240,8 @@ def _mul_idkronecker_scalar( _add_fns[(IdentityKronecker, IdentityKronecker)] = IdentityKronecker._add_idkronecker _sub_fns[(IdentityKronecker, IdentityKronecker)] = IdentityKronecker._sub_idkronecker -_mul_fns[(numbers.Number, IdentityKronecker)] = _mul_scalar_idkronecker -_mul_fns[(IdentityKronecker, numbers.Number)] = _mul_idkronecker_scalar +_mul_fns[(np.number, IdentityKronecker)] = _mul_scalar_idkronecker +_mul_fns[(IdentityKronecker, np.number)] = _mul_idkronecker_scalar _matmul_fns[(IdentityKronecker, Scaling)] = _matmul_idkronecker_scaling _matmul_fns[(Scaling, IdentityKronecker)] = _matmul_scaling_idkronecker @@ -258,8 +257,8 @@ def _matmul_matrix_scaling(matrix: Matrix, scaling: Scaling) -> Matrix: return Matrix(A=np.multiply(matrix.A, scaling.factors)) -_mul_fns[(Matrix, numbers.Number)] = lambda mat, scal: Matrix(A=scal * mat.A) -_mul_fns[(numbers.Number, Matrix)] = lambda scal, mat: Matrix(A=scal * mat.A) +_mul_fns[(Matrix, np.number)] = lambda mat, scal: Matrix(A=scal * mat.A) +_mul_fns[(np.number, Matrix)] = lambda scal, mat: Matrix(A=scal * mat.A) _matmul_fns[(Matrix, Matrix)] = Matrix._matmul_matrix _matmul_fns[(Scaling, Matrix)] = _matmul_scaling_matrix _matmul_fns[(Matrix, Scaling)] = _matmul_matrix_scaling @@ -401,12 +400,14 @@ def _apply( ] = None, ) -> Union[LinearOperator, NotImplementedType]: if np.ndim(op1) == 0: - key1 = numbers.Number + key1 = np.number + op1 = utils.as_numpy_scalar(op1) else: key1 = type(op1) if np.ndim(op2) == 0: - key2 = numbers.Number + key2 = np.number + op2 = utils.as_numpy_scalar(op1) else: key2 = type(op2) From a6542fb7e9e70414cd93dfcd5cdac7944904abd3 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 18 Nov 2021 16:54:43 +0100 Subject: [PATCH 66/74] Incoroporate PR comments --- src/probnum/linops/_kronecker.py | 12 ++++++------ src/probnum/linops/_scaling.py | 17 ++++++++--------- 2 files changed, 14 insertions(+), 15 deletions(-) diff --git a/src/probnum/linops/_kronecker.py b/src/probnum/linops/_kronecker.py index b14e8775a..f16305d9a 100644 --- a/src/probnum/linops/_kronecker.py +++ b/src/probnum/linops/_kronecker.py @@ -179,10 +179,10 @@ def _add_kronecker( self, other: "Kronecker" ) -> Union[NotImplementedType, "Kronecker"]: - if self.A == other.A: + if self.A is other.A or self.A == other.A: return Kronecker(A=self.A, B=self.B + other.B) - if self.B == other.B: + if self.B is other.B or self.B == other.B: return Kronecker(A=self.A + other.A, B=other.B) return NotImplemented @@ -191,10 +191,10 @@ def _sub_kronecker( self, other: "Kronecker" ) -> Union[NotImplementedType, "Kronecker"]: - if self.A == other.A: + if self.A is other.A or self.A == other.A: return Kronecker(A=self.A, B=self.B - other.B) - if self.B == other.B: + if self.B is other.B or self.B == other.B: return Kronecker(A=self.A - other.A, B=other.B) return NotImplemented @@ -481,10 +481,10 @@ def __init__(self, num_blocks: int, B: _utils.LinearOperatorLike): ), matmul=lambda x: _kronecker_matmul( self.A, self.B, x - ), # TODO: is there a more efficient way? + ), # TODO: can be implemented more efficiently rmatmul=lambda x: _kronecker_rmatmul( self.A, self.B, x - ), # TODO: is there a more efficient way? + ), # TODO: can be implemented more efficiently todense=lambda: np.kron( self.A.todense(cache=False), self.B.todense(cache=False) ), diff --git a/src/probnum/linops/_scaling.py b/src/probnum/linops/_scaling.py index 5d1b68ac4..c72ccfc66 100644 --- a/src/probnum/linops/_scaling.py +++ b/src/probnum/linops/_scaling.py @@ -124,8 +124,6 @@ def __init__( trace = lambda: self.shape[0] * self._scalar elif np.ndim(factors) == 1: - if np.all(factors == factors[0]): - self._scalar = factors[0] # Anisotropic scaling self._factors = np.asarray(factors, dtype=dtype) self._factors.setflags(write=False) @@ -201,6 +199,8 @@ def scalar(self) -> Optional[np.number]: @property def is_isotropic(self) -> bool: """Whether scaling is uniform / isotropic.""" + if self._scalar is None and np.all(self.factors == self.factors[0]): + self._scalar = self.factors[0] return self._scalar is not None def __eq__(self, other: _linear_operator.LinearOperator) -> bool: @@ -324,16 +324,15 @@ def _cond_isotropic(self, p: Union[None, int, float, str]) -> np.inexact: class Zero(_linear_operator.LinearOperator): def __init__(self, shape, dtype=np.float64): - zero = np.zeros(shape=(), dtype=dtype) - matmul = lambda x: zero * x - rmatmul = lambda x: zero * x - apply = lambda x, axis: zero * x + matmul = lambda x: np.zeros(x.shape, np.result_type(x, self.dtype)) + rmatmul = lambda x: np.zeros(x.shape, np.result_type(x, self.dtype)) + apply = lambda x, axis: np.zeros(x.shape, np.result_type(x, self.dtype)) todense = lambda: np.zeros(shape=shape, dtype=dtype) rank = lambda: np.intp(0) - eigvals = lambda: np.full(shape[0], zero, dtype=dtype) - det = lambda: (zero.astype(dtype, copy=False) ** shape[0]) + eigvals = lambda: np.zeros(shape=(shape[0],), dtype=dtype) + det = lambda: np.zeros(shape=(), dtype=dtype) - trace = lambda: zero + trace = lambda: np.zeros(shape=(), dtype=dtype) super().__init__( shape, From 9d69385a9dc0a57b543a505ae2619e332fdde941 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 18 Nov 2021 16:58:18 +0100 Subject: [PATCH 67/74] Remove remaining matrix_free_linalg --- .../test_markov/test_integrator/test_iwp.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_randprocs/test_markov/test_integrator/test_iwp.py b/tests/test_randprocs/test_markov/test_integrator/test_iwp.py index 7927c959c..d9099149c 100644 --- a/tests/test_randprocs/test_markov/test_integrator/test_iwp.py +++ b/tests/test_randprocs/test_markov/test_integrator/test_iwp.py @@ -288,7 +288,7 @@ def _setup( self, ): wiener_process_dimension = 1 # make tests compatible with some_normal_rv1, etc. - with config(matrix_free_linalg=True): + with config(matrix_free=True): self.transition = randprocs.markov.integrator.IntegratedWienerTransition( num_derivatives=2, wiener_process_dimension=wiener_process_dimension, @@ -297,7 +297,7 @@ def _setup( ) def test_discretise_values(self, ah_22_ibm, qh_22_ibm, dt): - with config(matrix_free_linalg=True): + with config(matrix_free=True): discrete_model = self.transition.discretise(dt=dt) np.testing.assert_allclose( discrete_model.state_trans_mat.todense(), ah_22_ibm @@ -307,7 +307,7 @@ def test_discretise_values(self, ah_22_ibm, qh_22_ibm, dt): ) def test_forward_rv_values(self, normal_rv3x3, diffusion, ah_22_ibm, qh_22_ibm, dt): - with config(matrix_free_linalg=True): + with config(matrix_free=True): rv, _ = self.transition.forward_rv( normal_rv3x3, t=0.0, dt=dt, _diffusion=diffusion ) @@ -320,7 +320,7 @@ def test_forward_rv_values(self, normal_rv3x3, diffusion, ah_22_ibm, qh_22_ibm, def test_forward_realization_values( self, normal_rv3x3, diffusion, ah_22_ibm, qh_22_ibm, dt ): - with config(matrix_free_linalg=True): + with config(matrix_free=True): real = normal_rv3x3.mean rv, _ = self.transition.forward_realization( real, t=0.0, dt=dt, _diffusion=diffusion From ce071e934b0d93066c5d02fbea52778d2722e9bf Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 18 Nov 2021 17:03:09 +0100 Subject: [PATCH 68/74] Fix docs --- src/probnum/_config.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/probnum/_config.py b/src/probnum/_config.py index 70f5fd585..b28814961 100644 --- a/src/probnum/_config.py +++ b/src/probnum/_config.py @@ -130,8 +130,8 @@ def register(self, key: str, default_value: Any, description: str) -> None: "matrix_free", False, ( - "If ``True``, wherever possible, ``LinearOperator``s are used instead " - "of arrays. ``LinearOperator``s define a matrix-vector product implicitly without " + "If ``True``, wherever possible, ``LinearOperator`` s are used instead " + "of arrays. ``LinearOperator`` s define a matrix-vector product implicitly without " "instantiating the full matrix in memory. This makes them memory- and runtime-efficient " "for linear algebra operations." ), @@ -139,7 +139,7 @@ def register(self, key: str, default_value: Any, description: str) -> None: ( "lazy_matrix_matrix_matmul", False, - "If ``True``, the matmul operation applied to two ``Matrix``-``LinearOperator``s will " + "If ``True``, the matmul operation applied to two ``Matrix``-``LinearOperator`` s will " "again yield a ``Matrix``-``LinearOperator`` with the computed matrix-matrix product," " instead of a ``ProductLinearOperator``.", ), From 0391203d7ae1059c76285ca2cfb6ace949dd7735 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 18 Nov 2021 17:08:37 +0100 Subject: [PATCH 69/74] pylint --- tests/test_linops/test_arithmetics.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py index f6d0b2541..2ea41b8a6 100644 --- a/tests/test_linops/test_arithmetics.py +++ b/tests/test_linops/test_arithmetics.py @@ -44,7 +44,7 @@ def _aslist(arg): def get_linop(linop_type): - # pylint: disable=too-many-return-statements + # pylint: disable=too-many-return-statements,too-complex,too-many-branches if linop_type is Kronecker: _A1 = np.ones((2, 2)) _B1 = np.ones((2, 3)) @@ -143,7 +143,7 @@ def test_matmul(): def test_mul(): - + # pylint: disable=else-if-used for (l_type, r_type) in _mul_fns.keys(): if ( l_type is Selection From 92fbda29d94e79ed6998c298ed9b5ed68e5b3e34 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 18 Nov 2021 17:10:49 +0100 Subject: [PATCH 70/74] Fix tests ("scalar") --- tests/test_linops/test_arithmetics.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py index 2ea41b8a6..ae8b5a948 100644 --- a/tests/test_linops/test_arithmetics.py +++ b/tests/test_linops/test_arithmetics.py @@ -103,7 +103,7 @@ def get_linop(linop_type): return _TypeCastLinearOperator( linop=Matrix(np.random.rand(4, 4)), dtype=np.float32 ) - elif isinstance(linop_type, str) and linop_type == "scalar": + elif linop_type is np.number: return 1.3579 else: raise TypeError(f"Don't know what to do with type {linop_type}.") @@ -133,9 +133,9 @@ def test_matmul(): res_linop = linop1 @ linop2 assert res_linop.ndim == 2 - if isinstance(l_type, str) and l_type == "scalar": + if l_type is np.number: assert res_linop.shape == linop2.shape - elif isinstance(r_type, str) and r_type == "scalar": + elif r_type is np.number: assert res_linop.shape == linop1.shape else: assert res_linop.shape[0] == linop1.shape[0] @@ -159,10 +159,10 @@ def test_mul(): for (linop1, linop2) in itertools.product(_aslist(linops1), _aslist(linops2)): - if isinstance(l_type, str) and l_type == "scalar": + if l_type is np.number: res_linop = linop1 * linop2 assert res_linop.shape == linop2.shape - elif isinstance(r_type, str) and r_type == "scalar": + elif r_type is np.number: res_linop = linop1 * linop2 assert res_linop.shape == linop1.shape else: From a220f9b637e9aa953508526c38eb554e8f9de826 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 18 Nov 2021 17:13:11 +0100 Subject: [PATCH 71/74] Fix small error --- src/probnum/linops/_arithmetic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/probnum/linops/_arithmetic.py b/src/probnum/linops/_arithmetic.py index b1803166d..6d993cf4d 100644 --- a/src/probnum/linops/_arithmetic.py +++ b/src/probnum/linops/_arithmetic.py @@ -407,7 +407,7 @@ def _apply( if np.ndim(op2) == 0: key2 = np.number - op2 = utils.as_numpy_scalar(op1) + op2 = utils.as_numpy_scalar(op2) else: key2 = type(op2) From 7cb002d86d62a5baeda2e60d57b15d5209ea9209 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 18 Nov 2021 17:15:06 +0100 Subject: [PATCH 72/74] Fix pylint part 2 --- tests/test_linops/test_arithmetics.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/tests/test_linops/test_arithmetics.py b/tests/test_linops/test_arithmetics.py index ae8b5a948..734927858 100644 --- a/tests/test_linops/test_arithmetics.py +++ b/tests/test_linops/test_arithmetics.py @@ -143,7 +143,6 @@ def test_matmul(): def test_mul(): - # pylint: disable=else-if-used for (l_type, r_type) in _mul_fns.keys(): if ( l_type is Selection @@ -165,13 +164,12 @@ def test_mul(): elif r_type is np.number: res_linop = linop1 * linop2 assert res_linop.shape == linop1.shape - else: - if linop1.shape != linop2.shape: - with pytest.raises(ValueError): - res_linop = linop1 * linop2 - else: + elif linop1.shape != linop2.shape: + with pytest.raises(ValueError): res_linop = linop1 * linop2 - assert res_linop.shape == linop1.shape == linop2.shape + else: + res_linop = linop1 * linop2 + assert res_linop.shape == linop1.shape == linop2.shape def test_add(): From 50fd28e9fd5f5b1c6e96fadc27cbad8de3eb3e33 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 18 Nov 2021 17:19:57 +0100 Subject: [PATCH 73/74] pylint part 3 --- src/probnum/linops/_linear_operator.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/probnum/linops/_linear_operator.py b/src/probnum/linops/_linear_operator.py index 696962311..86330e0d9 100644 --- a/src/probnum/linops/_linear_operator.py +++ b/src/probnum/linops/_linear_operator.py @@ -571,9 +571,6 @@ def _is_type_shape_dtype_equal(self, other: "LinearOperator") -> bool: and self.dtype == other.dtype ) - def __eq__(self, other: "LinearOperator") -> bool: - raise NotImplementedError - def __matmul__( self, other: BinaryOperandType ) -> Union["LinearOperator", np.ndarray]: From ac2c5e0e47421423cf21abb2e0de05d573c790f0 Mon Sep 17 00:00:00 2001 From: Jonathan Schmidt Date: Thu, 18 Nov 2021 17:27:41 +0100 Subject: [PATCH 74/74] PYLIIIIINT part 4 --- src/probnum/linops/_arithmetic_fallbacks.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/probnum/linops/_arithmetic_fallbacks.py b/src/probnum/linops/_arithmetic_fallbacks.py index 1ecdb289c..768e0f083 100644 --- a/src/probnum/linops/_arithmetic_fallbacks.py +++ b/src/probnum/linops/_arithmetic_fallbacks.py @@ -1,3 +1,4 @@ +"""Fallback-implementations of LinearOperator arithmetic.""" import functools import operator from typing import Tuple, Union