@@ -1,135 +1,359 @@
Loading
1 1
"""Kernel / covariance function."""
2 2
3 3
import abc
4 -
from typing import Generic, Optional, Tuple, TypeVar, Union
4 +
from typing import Optional, Union
5 5
6 6
import numpy as np
7 7
8 -
import probnum.utils as _utils
9 -
from probnum.typing import IntArgType, ShapeArgType, ShapeType
8 +
from probnum import utils as _pn_utils
9 +
from probnum.typing import ArrayLike, IntArgType, ShapeArgType, ShapeType
10 10
11 -
_InputType = TypeVar("InputType")
12 11
12 +
class Kernel(abc.ABC):
13 +
    r"""(Cross-)covariance function(s)
13 14
14 -
class Kernel(Generic[_InputType], abc.ABC):
15 -
    """Kernel / covariance function.
15 +
    Abstract base class representing one or multiple (cross-)covariance function(s),
16 +
    also known as kernels.
17 +
    A cross-covariance function
16 18
17 -
    Abstract base class for kernels / covariance functions. Kernels are a
18 -
    generalization of a positive-definite function or matrix. They
19 -
    typically define the covariance function of a random process and thus describe
20 -
    its spatial or temporal variation. If evaluated at two sets of points a kernel
21 -
    gives the covariance of the random process at these locations.
19 +
    .. math::
20 +
        :nowrap:
21 +
22 +
        \begin{equation}
23 +
            k_{fg} \colon
24 +
            \mathcal{X}^{d_\text{in}} \times \mathcal{X}^{d_\text{in}}
25 +
            \to \mathbb{R}
26 +
        \end{equation}
27 +
28 +
    is a function of two arguments :math:`x_0` and :math:`x_1`, which represents the
29 +
    covariance between two evaluations :math:`f(x_0)` and :math:`g(x_1)` of two
30 +
    scalar-valued random processes :math:`f` and :math:`g` (or, equivalently, two
31 +
    outputs :math:`h_i(x_0)` and :math:`h_j(x_1)` of a vector-valued random process
32 +
    :math:`h`).
33 +
    If :math:`f = g`, then the cross-covariance function is also referred to as a
34 +
    covariance function, in which case it must be symmetric and positive
35 +
    (semi-)definite.
36 +
37 +
    An instance of a :class:`Kernel` can compute multiple different (cross-)covariance
38 +
    functions on the same pair of inputs simultaneously. For instance, it can be used to
39 +
    compute the full covariance matrix
40 +
41 +
    .. math::
42 +
        :nowrap:
43 +
44 +
        \begin{equation}
45 +
            C^f \colon
46 +
            \mathcal{X}^{d_\text{in}} \times \mathcal{X}^{d_\text{in}}
47 +
            \to \mathbb{R}^{d_\text{out} \times d_\text{out}},
48 +
            C^f_{i j}(x_0, x_1) := k_{f_i f_j}(x_0, x_1)
49 +
        \end{equation}
50 +
51 +
    of the vector-valued random process :math:`f`. To this end, we understand any
52 +
    :class:`Kernel` with a non-empty :attr:`shape` as a tensor with the given
53 +
    :attr:`shape`, which contains different (cross-)covariance functions as its entries.
22 54
23 55
    Parameters
24 56
    ----------
25 57
    input_dim :
26 58
        Input dimension of the kernel.
27 -
    output_dim :
28 -
        Output dimension of the kernel.
59 +
    shape :
60 +
        If ``shape`` is set to ``()``, the :class:`Kernel` instance represents a
61 +
        single (cross-)covariance function. Otherwise, i.e. if ``shape`` is non-empty,
62 +
        the :class:`Kernel` instance represents a tensor of (cross-)covariance
63 +
        functions whose shape is given by ``shape``.
29 64
30 65
    Examples
31 66
    --------
32 -
    Kernels are implemented by subclassing this abstract base class.
33 -
34 -
    >>> from probnum.kernels import Kernel
35 -
    ...
36 -
    >>> class CustomLinearKernel(Kernel):
37 -
    ...
38 -
    ...     def __init__(self, constant=0.0):
39 -
    ...         self.constant = constant
40 -
    ...         super().__init__(input_dim=1, output_dim=1)
41 -
    ...
42 -
    ...     def __call__(self, x0, x1=None):
43 -
    ...         # Check and reshape inputs
44 -
    ...         x0, x1, kernshape = self._check_and_reshape_inputs(x0, x1)
45 -
    ...
46 -
    ...         # Compute kernel matrix
47 -
    ...         if x1 is None:
48 -
    ...             x1 = x0
49 -
    ...         kernmat = x0 @ x1.T + self.constant
50 -
    ...
51 -
    ...         return Kernel._reshape_kernelmatrix(kernmat, newshape=kernshape)
52 -
53 -
    We can now evaluate the kernel like so.
54 -
55 -
    >>> import numpy as np
56 -
    >>> k = CustomLinearKernel(constant=1.0)
57 -
    >>> k(np.linspace(0, 1, 4)[:, None])
58 -
    array([[1.        , 1.        , 1.        , 1.        ],
59 -
           [1.        , 1.11111111, 1.22222222, 1.33333333],
60 -
           [1.        , 1.22222222, 1.44444444, 1.66666667],
61 -
           [1.        , 1.33333333, 1.66666667, 2.        ]])
67 +
68 +
    >>> D = 3
69 +
    >>> k = pn.kernels.Linear(input_dim=D)
70 +
71 +
    Generate some input data.
72 +
73 +
    >>> xs = np.repeat(np.linspace(0, 1, 4)[:, None], D, axis=-1)
74 +
    >>> xs.shape
75 +
    (4, 3)
76 +
    >>> xs
77 +
    array([[0.        , 0.        , 0.        ],
78 +
           [0.33333333, 0.33333333, 0.33333333],
79 +
           [0.66666667, 0.66666667, 0.66666667],
80 +
           [1.        , 1.        , 1.        ]])
81 +
82 +
    We can compute kernel matrices like so.
83 +
84 +
    >>> k.matrix(xs)
85 +
    array([[0.        , 0.        , 0.        , 0.        ],
86 +
           [0.        , 0.33333333, 0.66666667, 1.        ],
87 +
           [0.        , 0.66666667, 1.33333333, 2.        ],
88 +
           [0.        , 1.        , 2.        , 3.        ]])
89 +
90 +
    Inputs to :meth:`Kernel.__call__` are broadcast according to the "kernel
91 +
    broadcasting" rules detailed in the "Notes" section of the :meth:`Kernel._call__`
92 +
    documentation.
93 +
94 +
    >>> k(xs[:, None, :], xs[None, :, :])  # same as `.matrix`
95 +
    array([[0.        , 0.        , 0.        , 0.        ],
96 +
           [0.        , 0.33333333, 0.66666667, 1.        ],
97 +
           [0.        , 0.66666667, 1.33333333, 2.        ],
98 +
           [0.        , 1.        , 2.        , 3.        ]])
99 +
100 +
    A shape of ``1`` along the last axis is broadcast to :attr:`input_dim`.
101 +
102 +
    >>> xs_d1 = xs[:, [0]]
103 +
    >>> xs_d1.shape
104 +
    (4, 1)
105 +
    >>> xs_d1
106 +
    array([[0.        ],
107 +
           [0.33333333],
108 +
           [0.66666667],
109 +
           [1.        ]])
110 +
    >>> k(xs_d1[:, None, :], xs_d1[None, :, :])  # same as `.matrix`
111 +
    array([[0.        , 0.        , 0.        , 0.        ],
112 +
           [0.        , 0.33333333, 0.66666667, 1.        ],
113 +
           [0.        , 0.66666667, 1.33333333, 2.        ],
114 +
           [0.        , 1.        , 2.        , 3.        ]])
115 +
    >>> k(xs[:, None, :], xs_d1[None, :, :])  # same as `.matrix`
116 +
    array([[0.        , 0.        , 0.        , 0.        ],
117 +
           [0.        , 0.33333333, 0.66666667, 1.        ],
118 +
           [0.        , 0.66666667, 1.33333333, 2.        ],
119 +
           [0.        , 1.        , 2.        , 3.        ]])
120 +
121 +
    No broadcasting is applied if both inputs have the same shape. For instance, one can
122 +
    efficiently compute just the diagonal of the kernel matrix via
123 +
124 +
    >>> k(xs, xs)
125 +
    array([0.        , 0.33333333, 1.33333333, 3.        ])
126 +
    >>> k(xs, None)  # x1 = None is an efficient way to set x1 == x0
127 +
    array([0.        , 0.33333333, 1.33333333, 3.        ])
128 +
129 +
    and the diagonal above the main diagonal of the kernel matrix is retrieved through
130 +
131 +
    >>> k(xs[:-1, :], xs[1:, :])
132 +
    array([0.        , 0.66666667, 2.        ])
62 133
    """
63 134
64 -
    # pylint: disable="invalid-name"
65 135
    def __init__(
66 136
        self,
67 137
        input_dim: IntArgType,
68 -
        output_dim: IntArgType = 1,
138 +
        shape: ShapeArgType = (),
69 139
    ):
70 -
        self._input_dim = np.int_(_utils.as_numpy_scalar(input_dim))
71 -
        self._output_dim = np.int_(_utils.as_numpy_scalar(output_dim))
140 +
        self._input_dim = int(input_dim)
141 +
142 +
        self._shape = _pn_utils.as_shape(shape)
72 143
73 144
    def __repr__(self) -> str:
74 145
        return f"<{self.__class__.__name__}>"
75 146
76 -
    @abc.abstractmethod
77 147
    def __call__(
78 -
        self, x0: _InputType, x1: Optional[_InputType] = None
79 -
    ) -> Union[np.ndarray, np.float_]:
80 -
        """Evaluate the kernel.
148 +
        self,
149 +
        x0: ArrayLike,
150 +
        x1: Optional[ArrayLike],
151 +
    ) -> Union[np.ndarray, np.floating]:
152 +
        """Evaluate the (cross-)covariance function(s).
81 153
82 -
        Computes the covariance function at ``x0`` and ``x1``. If the inputs have
83 -
        more than one dimension the covariance function is evaluated pairwise for all
84 -
        observations determined by the first dimension of ``x0`` and ``x1``. If
85 -
        only ``x0`` is given the kernel matrix :math:`K=k(X_0, X_0)` is computed.
154 +
        The inputs are broadcast to a common shape following the "kernel broadcasting"
155 +
        rules outlined in the "Notes" section.
86 156
87 157
        Parameters
88 158
        ----------
89 -
        x0 :
90 -
            *shape=(input_dim,) or (n0, input_dim)* -- First input.
91 -
        x1 :
92 -
            *shape=(input_dim,) or (n1, input_dim)* -- Second input.
159 +
        x0 : array-like
160 +
            An array of shape ``()`` or ``(Nn, ..., N2, N1, D_in)``, where ``D_in`` is
161 +
            either ``1`` or :attr:`input_dim`, whose entries will be passed to the first
162 +
            argument of the kernel.
163 +
        x1 : array-like
164 +
            An array of shape ``()`` or ``(Mm, ..., M2, M1, D_in)``, where ``D_in`` is
165 +
            either ``1`` or :attr:`input_dim`, whose entries will be
166 +
            passed to the second argument of the kernel. Can also be set to ``None``,
167 +
            in which case the function will behave as if ``x1 = x0``.
93 168
94 169
        Returns
95 170
        -------
96 -
        cov :
97 -
            *shape=(), (output_dim, output_dim) or (n0, n1) or (n0, n1, output_dim,
98 -
            output_dim)* -- Kernel evaluated at ``x0`` and ``x1`` or kernel matrix
99 -
            containing pairwise evaluations for all observations in ``x0`` (and ``x1``).
171 +
        k_x0_x1 : numpy.ndarray
172 +
            The (cross-)covariance function(s) evaluated at ``x0`` and ``x1``.
173 +
            If :attr:`shape` is ``()``, this method returns an array of shape
174 +
            ``(Lk, ..., L2, L1)`` whose entry at index ``(ik, ..., i2, i1)`` contains
175 +
            the evaluation of the (cross-)covariance function at the inputs
176 +
            ``x0[ik, ..., i2, i1, :] and ``x1[il, ..., i2, i1, :])``.
177 +
            For any non-empty :attr:`shape`, it returns an array of shape
178 +
            ``(Sl, ..., S2, S1, Lk, ..., L2, L1)``, where ``S`` is :attr:`shape`, whose
179 +
            entry at index ``(sl, ..., s2, s1, ik, ..., i2, i1)`` contains
180 +
            evaluation of the (cross-)covariance function at index ``(sl, ..., s2, s1)``
181 +
            at the inputs ``x0[ik, ..., i2, i1, :]`` and ``x1[ik, ..., i2, i1, :]``.
182 +
            Above, we assume that ``x0`` and ``x1`` have been broadcast according to the
183 +
            rules described in the "Notes" section.
184 +
185 +
        Raises
186 +
        ------
187 +
        ValueError
188 +
            If the inputs can not be "kernel broadcast" to a common shape.
189 +
190 +
        See Also
191 +
        --------
192 +
        matrix: Convenience function to compute a kernel matrix, i.e. a matrix of
193 +
            pairwise evaluations of the kernel on two sets of points.
194 +
195 +
        Notes
196 +
        -----
197 +
        A :class:`Kernel` operates on its two inputs by a slightly modified version of
198 +
        Numpy's broadcasting rules. First of all, the operation of the kernel is
199 +
        vectorized over all but the last dimension, applying standard broadcasting
200 +
        rules. An input with shape ``()`` is promoted to an input with shape ``(1,)``.
201 +
        Additionally, a ``1`` along the last axis of an input is interpreted as a (set
202 +
        of) point(s) with equal coordinates in all input dimensions, i.e. the inputs are
203 +
        broadcast to :attr:`input_dim` dimensions along the last axis. We refer to this
204 +
        modified set of broadcasting rules as "kernel broadcasting".
205 +
206 +
        Examples
207 +
        --------
208 +
        See documentation of class :class:`Kernel`.
100 209
        """
101 -
        raise NotImplementedError
102 210
103 -
    @property
104 -
    def input_dim(self) -> int:
105 -
        """Dimension of arguments of the covariance function.
211 +
        x0 = np.atleast_1d(x0)
212 +
213 +
        if x1 is not None:
214 +
            x1 = np.atleast_1d(x1)
215 +
216 +
        # Shape checking
217 +
        broadcast_input_shape = self._kernel_broadcast_shapes(x0, x1)
218 +
219 +
        # Evaluate the kernel
220 +
        k_x0_x1 = self._evaluate(x0, x1)
221 +
222 +
        assert k_x0_x1.shape == self._shape + broadcast_input_shape[:-1]
223 +
224 +
        return k_x0_x1
106 225
107 -
        The dimension of inputs to the covariance function :math:`k : \\mathbb{R}^{
108 -
        d_{in}} \\times \\mathbb{R}^{d_{in}} \\rightarrow
109 -
        \\mathbb{R}^{d_{out} \\times d_{out}}`.
226 +
    def matrix(
227 +
        self,
228 +
        x0: ArrayLike,
229 +
        x1: Optional[ArrayLike] = None,
230 +
    ) -> np.ndarray:
231 +
        """A convenience function for computing a kernel matrix for two sets of inputs.
232 +
233 +
        This is syntactic sugar for ``k(x0[:, None, :], x1[None, :, :])``. Hence, it
234 +
        computes the matrix of pairwise covariances between two sets of input points.
235 +
        If ``k`` represents a covariance function, then the resulting matrix will be
236 +
        symmetric positive (semi-)definite for ``x0 == x1``.
237 +
238 +
        Parameters
239 +
        ----------
240 +
        x0 : array-like
241 +
            First set of inputs to the (cross-)covariance function as an array of shape
242 +
            ``(M, D)``, where ``D`` is either 1 or :attr:`input_dim`.
243 +
        x1 : array-like
244 +
            Optional second set of inputs to the (cross-)covariance function as an array
245 +
            of shape ``(N, D)``, where ``D`` is either 1 or :attr:`input_dim`.
246 +
            If ``x1`` is not specified, the function behaves as if ``x1 = x0``.
247 +
248 +
        Returns
249 +
        -------
250 +
        kernmat : numpy.ndarray
251 +
            The matrix / stack of matrices containing the pairwise evaluations of the
252 +
            (cross-)covariance function(s) on ``x0`` and ``x1`` as an array of shape
253 +
            ``(M, N)`` if :attr:`shape` is ``()`` or
254 +
            ``(S[l - 1], ..., S[1], S[0], M, N)``, where ``S`` is :attr:`shape` if
255 +
            :attr:`shape` is non-empty.
256 +
257 +
        Raises
258 +
        ------
259 +
        ValueError
260 +
            If the shapes of the inputs don't match the specification.
261 +
262 +
        See Also
263 +
        --------
264 +
        __call__: Evaluate the kernel more flexibly.
265 +
266 +
        Examples
267 +
        --------
268 +
        See documentation of class :class:`Kernel`.
110 269
        """
270 +
271 +
        x0 = np.array(x0)
272 +
        x1 = x0 if x1 is None else np.array(x1)
273 +
274 +
        # Shape checking
275 +
        errmsg = (
276 +
            "`{argname}` must have shape `(N, D)` or `(D,)`, where `D` is the input "
277 +
            f"dimension of the kernel (D = {self.input_dim}), but an array with shape "
278 +
            "`{shape}` was given."
279 +
        )
280 +
281 +
        if not (1 <= x0.ndim <= 2 and x0.shape[-1] in (self.input_dim, 1)):
282 +
            raise ValueError(errmsg.format(argname="x0", shape=x0.shape))
283 +
284 +
        if not (1 <= x1.ndim <= 2 and x1.shape[-1] in (self.input_dim, 1)):
285 +
            raise ValueError(errmsg.format(argname="x1", shape=x1.shape))
286 +
287 +
        # Pairwise kernel evaluation
288 +
        return self(x0[:, None, :], x1[None, :, :])
289 +
290 +
    @property
291 +
    def input_dim(self) -> int:
292 +
        """Dimension of arguments of the covariance function."""
111 293
        return self._input_dim
112 294
113 295
    @property
114 -
    def output_dim(self) -> int:
115 -
        """Dimension of the evaluated covariance function.
296 +
    def shape(self) -> ShapeType:
297 +
        """If :attr:`shape` is ``()``, the :class:`Kernel` instance represents a
298 +
        single (cross-)covariance function. Otherwise, i.e. if :attr:`shape` is
299 +
        non-empty, the :class:`Kernel` instance represents a tensor of
300 +
        (cross-)covariance functions whose shape is given by ``shape``."""
301 +
        return self._shape
116 302
117 -
        The resulting evaluated kernel :math:`k(x_0, x_1) \\in
118 -
        \\mathbb{R}^{d_{out} \\times d_{out}}` has *shape=(output_dim,
119 -
        output_dim)*.
303 +
    @abc.abstractmethod
304 +
    def _evaluate(
305 +
        self,
306 +
        x0: ArrayLike,
307 +
        x1: Optional[ArrayLike],
308 +
    ) -> Union[np.ndarray, np.float_]:
309 +
        """Implementation of the kernel evaluation which is called after input checking.
310 +
311 +
        When implementing a particular kernel, the subclass should implement the kernel
312 +
        computation by overwriting this method. It is called by the
313 +
        :meth:`Kernel.__call__` method after applying input checking and promoting
314 +
        scalar inputs, i.e. inputs with shape ``()``, to arrays with shape ``(1,)``.
315 +
        The implementation of the kernel must implement the rules of kernel broadcasting
316 +
        which are outlined in the "Notes" section of the :meth:`Kernel.__call__`
317 +
        docstring. Note that the inputs are not automatically broadcast to a common
318 +
        shape, but it is guaranteed that the inputs can be broadcast to a common shape
319 +
        when applying the rules of "kernel broadcasting".
320 +
321 +
        Parameters
322 +
        ----------
323 +
        x0 : array-like
324 +
            An array of shape ``(Nn, ..., N2, N1, D_in)``, where ``D_in`` is either
325 +
            ``1`` or :attr:`input_dim`, whose entries will be passed to the first
326 +
            argument of the kernel.
327 +
        x1 : array-like
328 +
            An array of shape ``(Mm, ..., M2, M1, D_in)``, where ``D_in`` is either
329 +
            ``1`` or :attr:`input_dim`, whose entries will be passed to the second
330 +
            argument of the kernel. Can also be set to ``None``, in which case the
331 +
            method must behave as if ``x1 == x0``.
332 +
333 +
        Returns
334 +
        -------
335 +
        k_x0_x1 : numpy.ndarray
336 +
            The (cross-)covariance function(s) evaluated at ``x0`` and ``x1`` as an
337 +
            array of shape ``(S[l - 1], ..., S[1], S[0], Lk, ..., L2, L1)``,
338 +
            where ``S`` is :attr:`shape`, whose entry at index
339 +
            ``(sl, ..., s2, s1, ik, ..., i2, i1)`` contains the
340 +
            evaluation of the (cross-)covariance function at index ``(sl, ..., s2, s1)``
341 +
            at the inputs ``x0[ik, ..., i2, i1, :]`` and ``x1[ik, ..., i2, i1, :]``,
342 +
            where we assume that ``x0`` and ``x1`` have been broadcast according to the
343 +
            rules described in the "Notes" section of :meth:`Kernel.__call__`.
120 344
        """
121 -
        return self._output_dim
122 345
123 -
    def _check_and_reshape_inputs(
346 +
    def _kernel_broadcast_shapes(
124 347
        self,
125 -
        x0: _InputType,
126 -
        x1: Optional[_InputType] = None,
127 -
    ) -> Tuple[np.ndarray, Optional[np.ndarray], ShapeType]:
128 -
        """Check and transform inputs of the covariance function.
348 +
        x0: np.ndarray,
349 +
        x1: Optional[np.ndarray] = None,
350 +
    ) -> ShapeType:
351 +
        """Applies the "kernel broadcasting" rules to the input shapes.
129 352
130 -
        Checks the shape of the inputs to the covariance function and
131 -
        transforms the inputs into two-dimensional :class:`numpy.ndarray`s such that
132 -
        inputs are stacked row-wise.
353 +
        A description of the "kernel broadcasting" rules can be found in the docstring
354 +
        of :meth:`Kernel.__call__`.
355 +
        Appropriate exceptions are raised if the shapes can not be "kernel broadcast" to
356 +
        a common shape.
133 357
134 358
        Parameters
135 359
        ----------
@@ -140,135 +364,100 @@
Loading
140 364
141 365
        Returns
142 366
        -------
143 -
        x0 :
144 -
            First input to the covariance function.
145 -
        x1 :
146 -
            Second input to the covariance function.
147 -
        kernshape :
148 -
            Shape of the evaluation of the covariance function.
367 +
        broadcast_input_shape : tuple of int
368 +
            Shape of the inputs after applying "kernel broadcasting".
149 369
150 370
        Raises
151 371
        -------
152 -
        ValueError :
153 -
            If input shapes of x0 and x1 do not match the kernel input dimension or
154 -
            each other.
372 +
        ValueError
373 +
            If the inputs can not be "kernel broadcast" to a common shape.
155 374
        """
156 -
        # pylint: disable="too-many-boolean-expressions"
157 375
158 -
        # Check and promote shapes
159 -
        x0 = np.asarray(x0)
160 -
        if x1 is None:
161 -
            if (
162 -
                (x0.ndim == 0 and self.input_dim > 1)  # Scalar input
163 -
                or (x0.ndim == 1 and x0.shape[0] != self.input_dim)  # Vector input
164 -
                or (x0.ndim >= 2 and x0.shape[1] != self.input_dim)  # Matrix input
165 -
            ):
376 +
        err_msg = (
377 +
            "`{argname}` can not be broadcast to the kernel's input dimension "
378 +
            f"{self.input_dim} along the last axis, since "
379 +
            "`{argname}.shape = {shape}`."
380 +
        )
381 +
382 +
        if x0.shape[-1] not in (self.input_dim, 1):
383 +
            # This will never be called if the original input was a scalar
384 +
            raise ValueError(err_msg.format(argname="x0", shape=x0.shape))
385 +
386 +
        broadcast_input_shape = x0.shape
387 +
388 +
        if x1 is not None:
389 +
            if x1.shape[-1] not in (self.input_dim, 1):
390 +
                # This will never be called if the original input was a scalar
391 +
                raise ValueError(err_msg.format(argname="x1", shape=x1.shape))
392 +
393 +
            try:
394 +
                # Ironically, `np.broadcast_arrays` seems to be more efficient than
395 +
                # `np.broadcast_shapes`
396 +
                broadcast_input_shape = np.broadcast_arrays(x0, x1)[0].shape
397 +
            except ValueError as v:
166 398
                raise ValueError(
167 -
                    f"Argument shape x0.shape={x0.shape} does not match "
168 -
                    "kernel input dimension."
169 -
                )
399 +
                    f"The input arrays `x0` and `x1` with shapes {x0.shape} and "
400 +
                    f"{x1.shape} can not be broadcast to a common shape."
401 +
                ) from v
170 402
171 -
            # Determine correct shape for the kernel matrix as the output of __call__
172 -
            kernshape = self._get_shape_kernelmatrix(
173 -
                x0_shape=x0.shape, x1_shape=x0.shape
174 -
            )
403 +
        if broadcast_input_shape[-1] == 1:
404 +
            broadcast_input_shape = broadcast_input_shape[:-1] + (self.input_dim,)
175 405
176 -
            return np.atleast_2d(x0), None, kernshape
177 -
        else:
178 -
            x1 = np.asarray(x1)
179 -
            err_msg = (
180 -
                f"Argument shapes x0.shape={x0.shape} and x1.shape="
181 -
                f"{x1.shape} do not match kernel input dimension "
182 -
                f"{self.input_dim}. Try passing either two vectors or two "
183 -
                "matrices with the second dimension equal to the kernel input "
184 -
                "dimension."
185 -
            )
406 +
        return broadcast_input_shape
186 407
187 -
            # Promote unequal shapes
188 -
            if x0.ndim < 2 and x1.ndim == 2:
189 -
                x0 = np.atleast_2d(x0)
190 -
            if x1.ndim < 2 and x0.ndim == 2:
191 -
                x1 = np.atleast_2d(x1)
192 -
            if x0.ndim != x1.ndim:  # Shape mismatch
193 -
                raise ValueError(err_msg)
194 -
195 -
            # Check shapes
196 -
            if (
197 -
                (x0.ndim == 0 and self.input_dim > 1)  # Scalar input
198 -
                or (
199 -
                    x0.ndim == 1  # Vector input
200 -
                    and not (x0.shape[0] == x1.shape[0] == self.input_dim)
201 -
                )
202 -
                or (
203 -
                    x0.ndim == 2  # Matrix input
204 -
                    and not (x0.shape[1] == x1.shape[1] == self.input_dim)
205 -
                )
206 -
            ):
207 -
                raise ValueError(err_msg)
208 -
209 -
            # Determine correct shape for the kernel matrix as the output of __call__
210 -
            kernshape = self._get_shape_kernelmatrix(
211 -
                x0_shape=x0.shape, x1_shape=x1.shape
212 -
            )
408 +
    def _euclidean_inner_products(
409 +
        self, x0: np.ndarray, x1: Optional[np.ndarray]
410 +
    ) -> np.ndarray:
411 +
        """Implementation of the Euclidean inner product, which supports kernel
412 +
        broadcasting semantics."""
413 +
        prods = x0 ** 2 if x1 is None else x0 * x1
213 414
214 -
            return np.atleast_2d(x0), np.atleast_2d(x1), kernshape
415 +
        if prods.shape[-1] == 1:
416 +
            return self.input_dim * prods[..., 0]
215 417
216 -
    def _get_shape_kernelmatrix(
217 -
        self,
218 -
        x0_shape: ShapeArgType,
219 -
        x1_shape: ShapeArgType,
220 -
    ) -> ShapeType:
221 -
        """Determine the shape of the kernel matrix based on the given arguments.
418 +
        return np.sum(prods, axis=-1)
222 419
223 -
        Determine the correct shape of the covariance function evaluated at the given
224 -
        input arguments. If inputs are vectors the output is a numpy scalar if the
225 -
        output dimension of the kernel is 1, otherwise *shape=(output_dim,
226 -
        output_dim)*. If inputs represent multiple observations, then the resulting
227 -
        matrix has *shape=(n0, n1) or (n0, n1, output_dim, output_dim)*.
228 420
229 -
        Parameters
230 -
        ----------
231 -
        x0_shape :
232 -
            Shape of the first input to the covariance function.
233 -
        x1_shape :
234 -
            Shape of the second input to the covariance function.
235 -
        """
236 -
        if len(x0_shape) <= 1 and len(x1_shape) <= 1:
237 -
            if self.output_dim == 1:
238 -
                kern_shape = 0
239 -
            else:
240 -
                kern_shape = ()
241 -
        else:
242 -
            kern_shape = (x0_shape[0], x1_shape[0])
243 -
244 -
        if self.output_dim > 1:
245 -
            kern_shape += (
246 -
                self.output_dim,
247 -
                self.output_dim,
421 +
class IsotropicMixin(abc.ABC):  # pylint: disable=too-few-public-methods
422 +
    r"""Mixin for isotropic kernels.
423 +
424 +
    An isotropic kernel is a kernel which only depends on the Euclidean norm of the
425 +
    distance between the arguments, i.e.
426 +
427 +
    .. math ::
428 +
429 +
        k(x_0, x_1) = k(\lVert x_0 - x_1 \rVert_2).
430 +
431 +
    Hence, all isotropic kernels are stationary.
432 +
    """
433 +
434 +
    def _squared_euclidean_distances(
435 +
        self, x0: np.ndarray, x1: Optional[np.ndarray]
436 +
    ) -> np.ndarray:
437 +
        """Implementation of the squared Euclidean distance, which supports kernel
438 +
        broadcasting semantics."""
439 +
        if x1 is None:
440 +
            return np.zeros_like(  # pylint: disable=unexpected-keyword-arg
441 +
                x0,
442 +
                shape=x0.shape[:-1],
248 443
            )
249 444
250 -
        return _utils.as_shape(kern_shape)
445 +
        sqdiffs = (x0 - x1) ** 2
446 +
447 +
        if sqdiffs.shape[-1] == 1:
448 +
            return self.input_dim * sqdiffs[..., 0]
251 449
252 -
    @staticmethod
253 -
    def _reshape_kernelmatrix(
254 -
        kerneval: np.ndarray, newshape: ShapeArgType
450 +
        return np.sum(sqdiffs, axis=-1)
451 +
452 +
    def _euclidean_distances(
453 +
        self, x0: np.ndarray, x1: Optional[np.ndarray]
255 454
    ) -> np.ndarray:
256 -
        """Reshape the evaluation of the covariance function.
257 -
258 -
        Reshape the given evaluation of the covariance function to the correct shape,
259 -
        determined by the inputs x0 and x1. This method is designed to be called by
260 -
        subclasses of :class:`Kernel` in their :meth:`__call__` function to ensure
261 -
        the returned quantity has the correct shape independent of the implementation of
262 -
        the kernel.
263 -
264 -
        Parameters:
265 -
        -----------
266 -
        kerneval
267 -
            Covariance function evaluated at ``x0`` and ``x1``.
268 -
        newshape :
269 -
            New shape of the evaluation of the covariance function.
270 -
        """
271 -
        if newshape[0] == 0:
272 -
            return _utils.as_numpy_scalar(kerneval.squeeze())
273 -
        else:
274 -
            return kerneval.reshape(newshape)
455 +
        """Implementation of the Euclidean distance, which supports kernel
456 +
        broadcasting semantics."""
457 +
        if x1 is None:
458 +
            return np.zeros_like(  # pylint: disable=unexpected-keyword-arg
459 +
                x0,
460 +
                shape=x0.shape[:-1],
461 +
            )
462 +
463 +
        return np.sqrt(self._squared_euclidean_distances(x0, x1))

@@ -14,7 +14,7 @@
Loading
14 14
"""
15 15
16 16
import numbers
17 -
from typing import Iterable, Tuple, Union
17 +
from typing import Iterable, Sequence, Tuple, Union
18 18
19 19
import numpy as np
20 20
import scipy.sparse
@@ -29,6 +29,21 @@
Loading
29 29
# Argument Types
30 30
########################################################################################
31 31
32 +
try:
33 +
    from numpy.typing import ArrayLike  # pylint: disable=unused-import
34 +
except ImportError:
35 +
    # This is a very basic implementation which may be altered to align more with
36 +
    # NumPy's definition of an array-like
37 +
    _ArrayLikeScalar = Union[np.generic, bool, int, float, complex, str, bytes]
38 +
    _ArrayLikeNestedSequence = Union[  # Dirty hack to avoid recursive types
39 +
        Sequence[_ArrayLikeScalar],
40 +
        Sequence[Sequence[_ArrayLikeScalar]],
41 +
        Sequence[Sequence[Sequence[_ArrayLikeScalar]]],
42 +
        Sequence[Sequence[Sequence[Sequence[_ArrayLikeScalar]]]],
43 +
        Sequence[Sequence[Sequence[Sequence[Sequence[_ArrayLikeScalar]]]]],
44 +
    ]
45 +
    ArrayLike = Union[_ArrayLikeScalar, _ArrayLikeNestedSequence]
46 +
32 47
IntArgType = Union[int, numbers.Integral, np.integer]
33 48
FloatArgType = Union[float, numbers.Real, np.floating]
34 49

@@ -48,11 +48,16 @@
Loading
48 48
    def __init__(
49 49
        self,
50 50
        input_dim: IntArgType,
51 -
        output_dim: IntArgType,
51 +
        output_dim: Optional[IntArgType],
52 52
        dtype: DTypeArgType,
53 53
    ):
54 54
        self._input_dim = np.int_(_utils.as_numpy_scalar(input_dim))
55 -
        self._output_dim = np.int_(_utils.as_numpy_scalar(output_dim))
55 +
56 +
        self._output_dim = None
57 +
58 +
        if output_dim is not None:
59 +
            self._output_dim = np.int_(_utils.as_numpy_scalar(output_dim))
60 +
56 61
        self._dtype = np.dtype(dtype)
57 62
58 63
    def __repr__(self) -> str:
@@ -155,6 +160,65 @@
Loading
155 160
        """  # pylint: disable=trailing-whitespace
156 161
        raise NotImplementedError
157 162
163 +
    def covmatrix(
164 +
        self, args0: _InputType, args1: Optional[_InputType] = None
165 +
    ) -> _OutputType:
166 +
        """A convenience function for the covariance matrix of two sets of inputs.
167 +
168 +
        This is syntactic sugar for ``proc.cov(x0[:, None, :], x1[None, :, :])``. Hence,
169 +
        it computes the matrix of pairwise covariances between two sets of input points.
170 +
171 +
        Parameters
172 +
        ----------
173 +
        x0 : array-like
174 +
            First set of inputs to the covariance function as an array of shape
175 +
            ``(M, D)``, where ``D`` is either 1 or :attr:`input_dim`.
176 +
        x1 : array-like
177 +
            Optional second set of inputs to the covariance function as an array
178 +
            of shape ``(N, D)``, where ``D`` is either 1 or :attr:`input_dim`.
179 +
            If ``x1`` is not specified, the function behaves as if ``x1 = x0``.
180 +
181 +
        Returns
182 +
        -------
183 +
        kernmat : numpy.ndarray
184 +
            The matrix / stack of matrices containing the pairwise evaluations of the
185 +
            covariance function(s) on ``x0`` and ``x1`` as an array of shape
186 +
            ``(M, N)`` if :attr:`shape` is ``()`` or
187 +
            ``(S[l - 1], ..., S[1], S[0], M, N)``, where ``S`` is :attr:`shape` if
188 +
            :attr:`shape` is non-empty.
189 +
190 +
        Raises
191 +
        ------
192 +
        ValueError
193 +
            If the shapes of the inputs don't match the specification.
194 +
195 +
        See Also
196 +
        --------
197 +
        RandomProcess.cov: Evaluate the kernel more flexibly.
198 +
199 +
        Examples
200 +
        --------
201 +
        See documentation of class :class:`Kernel`.
202 +
        """
203 +
        args0 = np.array(args0)
204 +
        args1 = args0 if args1 is None else np.array(args1)
205 +
206 +
        # Shape checking
207 +
        errmsg = (
208 +
            "`{argname}` must have shape `(N, D)` or `(D,)`, where `D` is the input "
209 +
            f"dimension of the random process (D = {self.input_dim}), but an array "
210 +
            "with shape `{shape}` was given."
211 +
        )
212 +
213 +
        if not (1 <= args0.ndim <= 2 and args0.shape[-1] == self.input_dim):
214 +
            raise ValueError(errmsg.format(argname="args0", shape=args0.shape))
215 +
216 +
        if not (1 <= args1.ndim <= 2 and args1.shape[-1] == self.input_dim):
217 +
            raise ValueError(errmsg.format(argname="args1", shape=args1.shape))
218 +
219 +
        # Pairwise kernel evaluation
220 +
        return self.cov(args0[:, None, :], args1[None, :, :])
221 +
158 222
    def var(self, args: _InputType) -> _OutputType:
159 223
        """Variance function.
160 224
@@ -173,18 +237,17 @@
Loading
173 237
            process at ``args``.
174 238
        """
175 239
        try:
176 -
            cov = self.cov(args0=args)
177 -
            if cov.ndim < 2:
178 -
                return cov
179 -
            elif cov.ndim == 2:
180 -
                return np.diag(cov)
181 -
            else:
182 -
                return np.vstack(
183 -
                    [np.diagonal(cov[:, :, i, i]) for i in range(self.output_dim)]
184 -
                ).T
240 +
            var = self.cov(args0=args)
185 241
        except NotImplementedError as exc:
186 242
            raise NotImplementedError from exc
187 243
244 +
        if var.ndim == args.ndim - 1:
245 +
            return var
246 +
247 +
        assert var.ndim == args.ndim + 1 and var.shape[-2:] == 2 * (self.output_dim,)
248 +
249 +
        return np.diagonal(var, axis1=-2, axis2=-1)
250 +
188 251
    def std(self, args: _InputType) -> _OutputType:
189 252
        """Standard deviation function.
190 253

@@ -3,21 +3,22 @@
Loading
3 3
from typing import Optional
4 4
5 5
import numpy as np
6 -
import scipy.spatial.distance
7 6
8 7
import probnum.utils as _utils
9 8
from probnum.typing import IntArgType, ScalarArgType
10 9
11 -
from ._kernel import Kernel
10 +
from ._kernel import IsotropicMixin, Kernel
12 11
13 -
_InputType = np.ndarray
14 12
13 +
class ExpQuad(Kernel, IsotropicMixin):
14 +
    r"""Exponentiated quadratic / RBF kernel.
15 15
16 -
class ExpQuad(Kernel[_InputType]):
17 -
    """Exponentiated quadratic / RBF kernel.
16 +
    Covariance function defined by
18 17
19 -
    Covariance function defined by :math:`k(x_0, x_1) = \\exp \\big(-\\frac{\\lVert
20 -
    x_0 - x_1 \\rVert^2}{2l^2}\\big)`. This kernel is also known as the squared
18 +
    .. math ::
19 +
        k(x_0, x_1) = \exp \left( -\frac{\lVert x_0 - x_1 \rVert_2^2}{2 l^2} \right).
20 +
21 +
    This kernel is also known as the squared
21 22
    exponential or radial basis function kernel.
22 23
23 24
    Parameters
@@ -25,8 +26,8 @@
Loading
25 26
    input_dim :
26 27
        Input dimension of the kernel.
27 28
    lengthscale
28 -
        Lengthscale of the kernel. Describes the input scale on which the process
29 -
        varies.
29 +
        Lengthscale :math:`l` of the kernel. Describes the input scale on which the
30 +
        process varies.
30 31
31 32
    See Also
32 33
    --------
@@ -38,7 +39,8 @@
Loading
38 39
    >>> import numpy as np
39 40
    >>> from probnum.kernels import ExpQuad
40 41
    >>> K = ExpQuad(input_dim=1, lengthscale=0.1)
41 -
    >>> K(np.linspace(0, 1, 3)[:, None])
42 +
    >>> xs = np.linspace(0, 1, 3)[:, None]
43 +
    >>> K.matrix(xs)
42 44
    array([[1.00000000e+00, 3.72665317e-06, 1.92874985e-22],
43 45
           [3.72665317e-06, 1.00000000e+00, 3.72665317e-06],
44 46
           [1.92874985e-22, 3.72665317e-06, 1.00000000e+00]])
@@ -46,23 +48,12 @@
Loading
46 48
47 49
    def __init__(self, input_dim: IntArgType, lengthscale: ScalarArgType = 1.0):
48 50
        self.lengthscale = _utils.as_numpy_scalar(lengthscale)
49 -
        super().__init__(input_dim=input_dim, output_dim=1)
50 -
51 -
    def __call__(self, x0: _InputType, x1: Optional[_InputType] = None) -> np.ndarray:
51 +
        super().__init__(input_dim=input_dim)
52 52
53 -
        x0, x1, kernshape = self._check_and_reshape_inputs(x0, x1)
54 -
55 -
        # Compute pairwise euclidean distances ||x0 - x1|| / l
53 +
    def _evaluate(self, x0: np.ndarray, x1: Optional[np.ndarray] = None) -> np.ndarray:
56 54
        if x1 is None:
57 -
            pdists = scipy.spatial.distance.squareform(
58 -
                scipy.spatial.distance.pdist(x0 / self.lengthscale, metric="euclidean")
59 -
            )
60 -
        else:
61 -
            pdists = scipy.spatial.distance.cdist(
62 -
                x0 / self.lengthscale, x1 / self.lengthscale, metric="euclidean"
63 -
            )
64 -
65 -
        # Compute kernel matrix
66 -
        kernmat = np.exp(-(pdists ** 2) / 2.0)
55 +
            return np.ones_like(x0[..., 0])
67 56
68 -
        return Kernel._reshape_kernelmatrix(kernmat, newshape=kernshape)
57 +
        return np.exp(
58 +
            -self._squared_euclidean_distances(x0, x1) / (2.0 * self.lengthscale ** 2)
59 +
        )

@@ -9,34 +9,41 @@
Loading
9 9
import probnum.utils as _utils
10 10
from probnum.typing import IntArgType, ScalarArgType
11 11
12 -
from ._kernel import Kernel
12 +
from ._kernel import IsotropicMixin, Kernel
13 13
14 -
_InputType = np.ndarray
15 14
15 +
class Matern(Kernel, IsotropicMixin):
16 +
    r"""Matern kernel.
16 17
17 -
class Matern(Kernel[_InputType]):
18 -
    """Matern kernel.
18 +
    Covariance function defined by
19 19
20 -
    Covariance function defined by :math:`k(x_0, x_1) = \\frac{1}{\\Gamma(\\nu)2^{
21 -
    \\nu-1}}\\big(\\frac{\\sqrt{2\\nu}}{l} \\lVert x_0 , x_1\\rVert \\big)^\\nu
22 -
    K_\\nu\\big(\\frac{\\sqrt{2\\nu}}{l} \\lVert x_0 , x_1 \\rVert \\big)`, where
23 -
    :math:`K_\\nu` is a modified Bessel function. The Matern
24 -
    kernel generalizes the :class:`~probnum.kernels.ExpQuad` kernel
25 -
    via its additional parameter :math:`\\nu` controlling the smoothness of the
26 -
    function. For :math:`\\nu \\rightarrow \\infty` the Matern kernel converges to
27 -
    the :class:`~probnum.kernels.ExpQuad` kernel. A Gaussian process
28 -
    with Matern covariance function is :math:`\\lceil \\nu \\rceil - 1` times
29 -
    differentiable.
20 +
    .. math::
21 +
        :nowrap:
22 +
23 +
        \begin{equation}
24 +
            k(x_0, x_1)
25 +
            =
26 +
            \frac{1}{\Gamma(\nu) 2^{\nu - 1}}
27 +
            \left( \frac{\sqrt{2 \nu}}{l} \lVert x_0 - x_1 \rVert_2 \right)^\nu
28 +
            K_\nu \left( \frac{\sqrt{2 \nu}}{l} \lVert x_0 - x_1 \rVert_2 \right),
29 +
        \end{equation}
30 +
31 +
    where :math:`K_\nu` is a modified Bessel function. The Matern kernel generalizes the
32 +
    :class:`~probnum.kernels.ExpQuad` kernel via its additional parameter :math:`\nu`
33 +
    controlling the smoothness of the function. For :math:`\nu \rightarrow \infty`
34 +
    the Matern kernel converges to the :class:`~probnum.kernels.ExpQuad` kernel. A
35 +
    Gaussian process with Matern covariance function is :math:`\lceil \nu \rceil - 1`
36 +
    times differentiable.
30 37
31 38
    Parameters
32 39
    ----------
33 40
    input_dim :
34 41
        Input dimension of the kernel.
35 42
    lengthscale :
36 -
        Lengthscale of the kernel. Describes the input scale on which the process
37 -
        varies.
43 +
        Lengthscale :math:`l` of the kernel. Describes the input scale on which the
44 +
        process varies.
38 45
    nu :
39 -
        Hyperparameter controlling differentiability.
46 +
        Hyperparameter :math:`\nu` controlling differentiability.
40 47
41 48
    See Also
42 49
    --------
@@ -47,7 +54,8 @@
Loading
47 54
    >>> import numpy as np
48 55
    >>> from probnum.kernels import Matern
49 56
    >>> K = Matern(input_dim=1, lengthscale=0.1, nu=2.5)
50 -
    >>> K(np.linspace(0, 1, 3)[:, None])
57 +
    >>> xs = np.linspace(0, 1, 3)[:, None]
58 +
    >>> K.matrix(xs)
51 59
    array([[1.00000000e+00, 7.50933789e-04, 3.69569622e-08],
52 60
           [7.50933789e-04, 1.00000000e+00, 7.50933789e-04],
53 61
           [3.69569622e-08, 7.50933789e-04, 1.00000000e+00]])
@@ -59,7 +67,6 @@
Loading
59 67
        lengthscale: ScalarArgType = 1.0,
60 68
        nu: ScalarArgType = 1.5,
61 69
    ):
62 -
        # pylint: disable="invalid-name"
63 70
        self.lengthscale = _utils.as_numpy_scalar(lengthscale)
64 71
        if not self.lengthscale > 0:
65 72
            raise ValueError(f"Lengthscale l={self.lengthscale} must be positive.")
@@ -67,44 +74,35 @@
Loading
67 74
        if not self.nu > 0:
68 75
            raise ValueError(f"Hyperparameter nu={self.nu} must be positive.")
69 76
70 -
        super().__init__(input_dim=input_dim, output_dim=1)
71 -
72 -
    def __call__(self, x0: _InputType, x1: Optional[_InputType] = None) -> np.ndarray:
77 +
        super().__init__(input_dim=input_dim)
73 78
74 -
        x0, x1, kernshape = self._check_and_reshape_inputs(x0, x1)
75 -
76 -
        # Compute pairwise euclidean distances ||x0 - x1|| / l
77 -
        if x1 is None:
78 -
            pdists = scipy.spatial.distance.squareform(
79 -
                scipy.spatial.distance.pdist(x0 / self.lengthscale, metric="euclidean")
80 -
            )
81 -
        else:
82 -
            pdists = scipy.spatial.distance.cdist(
83 -
                x0 / self.lengthscale, x1 / self.lengthscale, metric="euclidean"
84 -
            )
79 +
    def _evaluate(self, x0: np.ndarray, x1: Optional[np.ndarray] = None) -> np.ndarray:
80 +
        distances = self._euclidean_distances(x0, x1)
85 81
86 82
        # Kernel matrix computation dependent on differentiability
87 83
        if self.nu == 0.5:
88 -
            kernmat = np.exp(-pdists)
89 -
        elif self.nu == 1.5:
90 -
            scaled_pdists = np.sqrt(3) * pdists
91 -
            kernmat = (1.0 + scaled_pdists) * np.exp(-scaled_pdists)
92 -
        elif self.nu == 2.5:
93 -
            scaled_pdists = np.sqrt(5) * pdists
94 -
            kernmat = (1.0 + scaled_pdists + scaled_pdists ** 2 / 3.0) * np.exp(
95 -
                -scaled_pdists
96 -
            )
97 -
        elif self.nu == np.inf:
98 -
            kernmat = np.exp(-(pdists ** 2) / 2.0)
99 -
        else:
100 -
            # The modified Bessel function K_nu is not defined for z=0
101 -
            pdists[pdists == 0.0] += np.finfo(float).eps
102 -
            scaled_pdists = np.sqrt(2 * self.nu) * pdists
103 -
            kernmat = (
104 -
                2 ** (1.0 - self.nu)
105 -
                / scipy.special.gamma(self.nu)
106 -
                * scaled_pdists ** self.nu
107 -
                * scipy.special.kv(self.nu, scaled_pdists)
84 +
            return np.exp(-1.0 / self.lengthscale * distances)
85 +
86 +
        if self.nu == 1.5:
87 +
            scaled_distances = -np.sqrt(3) / self.lengthscale * distances
88 +
            return (1.0 + scaled_distances) * np.exp(-scaled_distances)
89 +
90 +
        if self.nu == 2.5:
91 +
            scaled_distances = np.sqrt(5) / self.lengthscale * distances
92 +
            return (1.0 + scaled_distances + scaled_distances ** 2 / 3.0) * np.exp(
93 +
                -scaled_distances
108 94
            )
109 95
110 -
        return Kernel._reshape_kernelmatrix(kernmat, newshape=kernshape)
96 +
        if self.nu == np.inf:
97 +
            return np.exp(-1.0 / (2.0 * self.lengthscale ** 2) * distances ** 2)
98 +
99 +
        # The modified Bessel function K_nu is not defined for z=0
100 +
        distances = np.maximum(distances, np.finfo(distances.dtype).eps)
101 +
102 +
        scaled_distances = np.sqrt(2 * self.nu) / self.lengthscale * distances
103 +
        return (
104 +
            2 ** (1.0 - self.nu)
105 +
            / scipy.special.gamma(self.nu)
106 +
            * scaled_distances ** self.nu
107 +
            * scipy.special.kv(self.nu, scaled_distances)
108 +
        )

@@ -3,34 +3,41 @@
Loading
3 3
from typing import Optional
4 4
5 5
import numpy as np
6 -
import scipy.spatial.distance
7 6
8 7
import probnum.utils as _utils
9 8
from probnum.typing import IntArgType, ScalarArgType
10 9
11 -
from ._kernel import Kernel
10 +
from ._kernel import IsotropicMixin, Kernel
12 11
13 -
_InputType = np.ndarray
14 12
13 +
class RatQuad(Kernel, IsotropicMixin):
14 +
    r"""Rational quadratic kernel.
15 15
16 -
class RatQuad(Kernel[_InputType]):
17 -
    """Rational quadratic kernel.
16 +
    Covariance function defined by
18 17
19 -
    Covariance function defined by :math:`k(x_0, x_1) = \\big(1 + \\frac{\\lVert x_0 -
20 -
    x_1 \\rVert^2}{2\\alpha l^2}\\big)^{-\\alpha}`, where :math:`\\alpha > 0`. For
21 -
    :math:`\\alpha \\rightarrow \\infty` the rational quadratic kernel converges to the
22 -
    :class:`~probnum.kernels.ExpQuad` kernel.
18 +
    .. math::
19 +
        :nowrap:
20 +
21 +
        \begin{equation}
22 +
            k(x_0, x_1)
23 +
            = \left(
24 +
                1 + \frac{\lVert x_0 - x_1 \rVert_2^2}{2 \alpha l^2}
25 +
            \right)^{-\alpha},
26 +
        \end{equation}
27 +
28 +
    where :math:`\alpha > 0`. For :math:`\alpha \rightarrow \infty` the rational
29 +
    quadratic kernel converges to the :class:`~probnum.kernels.ExpQuad` kernel.
23 30
24 31
    Parameters
25 32
    ----------
26 33
    input_dim :
27 34
        Input dimension of the kernel.
28 35
    lengthscale :
29 -
        Lengthscale of the kernel. Describes the input scale on which the process
30 -
        varies.
36 +
        Lengthscale :math:`l` of the kernel. Describes the input scale on which the
37 +
        process varies.
31 38
    alpha :
32 -
        Scale mixture. Positive constant determining the weighting between different
33 -
        lengthscales.
39 +
        Scale mixture :math:`\alpha`. Positive constant determining the weighting
40 +
        between different lengthscales.
34 41
35 42
    See Also
36 43
    --------
@@ -41,7 +48,8 @@
Loading
41 48
    >>> import numpy as np
42 49
    >>> from probnum.kernels import RatQuad
43 50
    >>> K = RatQuad(input_dim=1, lengthscale=0.1, alpha=3)
44 -
    >>> K(np.linspace(0, 1, 3)[:, None])
51 +
    >>> xs = np.linspace(0, 1, 3)[:, None]
52 +
    >>> K(xs[:, None, :], xs[None, :, :])
45 53
    array([[1.00000000e+00, 7.25051190e-03, 1.81357765e-04],
46 54
           [7.25051190e-03, 1.00000000e+00, 7.25051190e-03],
47 55
           [1.81357765e-04, 7.25051190e-03, 1.00000000e+00]])
@@ -57,23 +65,16 @@
Loading
57 65
        self.alpha = _utils.as_numpy_scalar(alpha)
58 66
        if not self.alpha > 0:
59 67
            raise ValueError(f"Scale mixture alpha={self.alpha} must be positive.")
60 -
        super().__init__(input_dim=input_dim, output_dim=1)
61 -
62 -
    def __call__(self, x0: _InputType, x1: Optional[_InputType] = None) -> np.ndarray:
63 -
64 -
        x0, x1, kernshape = self._check_and_reshape_inputs(x0, x1)
68 +
        super().__init__(input_dim=input_dim)
65 69
66 -
        # Compute pairwise euclidean distances ||x0 - x1|| / l
70 +
    def _evaluate(self, x0: np.ndarray, x1: Optional[np.ndarray] = None) -> np.ndarray:
67 71
        if x1 is None:
68 -
            pdists = scipy.spatial.distance.squareform(
69 -
                scipy.spatial.distance.pdist(x0 / self.lengthscale, metric="euclidean")
70 -
            )
71 -
        else:
72 -
            pdists = scipy.spatial.distance.cdist(
73 -
                x0 / self.lengthscale, x1 / self.lengthscale, metric="euclidean"
74 -
            )
72 +
            return np.ones_like(x0[..., 0])
75 73
76 -
        # Kernel matrix
77 -
        kernmat = (1.0 + pdists ** 2 / (2.0 * self.alpha)) ** (-self.alpha)
78 -
79 -
        return Kernel._reshape_kernelmatrix(kernmat, newshape=kernshape)
74 +
        return (
75 +
            1.0
76 +
            + (
77 +
                self._squared_euclidean_distances(x0, x1)
78 +
                / (2.0 * self.alpha * self.lengthscale ** 2)
79 +
            )
80 +
        ) ** -self.alpha

@@ -6,7 +6,7 @@
Loading
6 6
"""
7 7
8 8
from ._exponentiated_quadratic import ExpQuad
9 -
from ._kernel import Kernel
9 +
from ._kernel import IsotropicMixin, Kernel
10 10
from ._linear import Linear
11 11
from ._matern import Matern
12 12
from ._polynomial import Polynomial
@@ -16,6 +16,7 @@
Loading
16 16
# Public classes and functions. Order is reflected in documentation.
17 17
__all__ = [
18 18
    "Kernel",
19 +
    "IsotropicMixin",
19 20
    "WhiteNoise",
20 21
    "Linear",
21 22
    "Polynomial",
@@ -26,6 +27,8 @@
Loading
26 27
27 28
# Set correct module paths. Corrects links and module paths in documentation.
28 29
Kernel.__module__ = "probnum.kernels"
30 +
IsotropicMixin.__module__ = "probnum.kernels"
31 +
29 32
WhiteNoise.__module__ = "probnum.kernels"
30 33
Linear.__module__ = "probnum.kernels"
31 34
Polynomial.__module__ = "probnum.kernels"

@@ -9,13 +9,14 @@
Loading
9 9
10 10
from ._kernel import Kernel
11 11
12 -
_InputType = np.ndarray
13 12
13 +
class Polynomial(Kernel):
14 +
    r"""Polynomial kernel.
14 15
15 -
class Polynomial(Kernel[_InputType]):
16 -
    """Polynomial kernel.
16 +
    Covariance function defined by
17 17
18 -
    Covariance function defined by :math:`k(x_0, x_1) = (x_0^\\top x_1 + c)^q`.
18 +
    .. math ::
19 +
        k(x_0, x_1) = (x_0^\top x_1 + c)^q.
19 20
20 21
    Parameters
21 22
    ----------
@@ -35,7 +36,8 @@
Loading
35 36
    >>> import numpy as np
36 37
    >>> from probnum.kernels import Polynomial
37 38
    >>> K = Polynomial(input_dim=2, constant=1.0, exponent=3)
38 -
    >>> K(np.array([[1, -1], [-1, 0]]))
39 +
    >>> xs = np.array([[1, -1], [-1, 0]])
40 +
    >>> K.matrix(xs)
39 41
    array([[27.,  0.],
40 42
           [ 0.,  8.]])
41 43
    """
@@ -48,15 +50,7 @@
Loading
48 50
    ):
49 51
        self.constant = _utils.as_numpy_scalar(constant)
50 52
        self.exponent = _utils.as_numpy_scalar(exponent)
51 -
        super().__init__(input_dim=input_dim, output_dim=1)
53 +
        super().__init__(input_dim=input_dim)
52 54
53 -
    def __call__(self, x0: _InputType, x1: Optional[_InputType] = None) -> np.ndarray:
54 -
55 -
        x0, x1, kernshape = self._check_and_reshape_inputs(x0, x1)
56 -
57 -
        # Compute kernel matrix
58 -
        if x1 is None:
59 -
            x1 = x0
60 -
        kernmat = (x0 @ x1.T + self.constant) ** self.exponent
61 -
62 -
        return Kernel._reshape_kernelmatrix(kernmat, newshape=kernshape)
55 +
    def _evaluate(self, x0: np.ndarray, x1: Optional[np.ndarray] = None) -> np.ndarray:
56 +
        return (self._euclidean_inner_products(x0, x1) + self.constant) ** self.exponent

@@ -53,7 +53,7 @@
Loading
53 53
           [-0.52972512],
54 54
           [ 0.0674298 ],
55 55
           [ 0.72066223]])
56 -
    >>> gp.cov(x)
56 +
    >>> gp.covmatrix(x)
57 57
    array([[1.        , 0.8824969 , 0.60653066, 0.32465247, 0.13533528],
58 58
           [0.8824969 , 1.        , 0.8824969 , 0.60653066, 0.32465247],
59 59
           [0.60653066, 0.8824969 , 1.        , 0.8824969 , 0.60653066],
@@ -71,16 +71,26 @@
Loading
71 71
                "The covariance functions must be implemented as a " "`Kernel`."
72 72
            )
73 73
74 +
        if cov.shape != () and not (
75 +
            len(cov.shape) == 2 and cov.shape[0] == cov.shape[1]
76 +
        ):
77 +
            raise ValueError(
78 +
                "Only kernels with shape `()` or `(D, D)` are allowed as covariance "
79 +
                "functions of random processes."
80 +
            )
81 +
74 82
        self._meanfun = mean
75 83
        self._covfun = cov
76 84
        super().__init__(
77 85
            input_dim=cov.input_dim,
78 -
            output_dim=cov.output_dim,
86 +
            output_dim=None if cov.shape == () else cov.shape[0],
79 87
            dtype=np.dtype(np.float_),
80 88
        )
81 89
82 90
    def __call__(self, args: _InputType) -> randvars.Normal:
83 -
        return randvars.Normal(mean=self.mean(args), cov=self.cov(args))
91 +
        return randvars.Normal(
92 +
            mean=np.array(self.mean(args), copy=False), cov=self.covmatrix(args)
93 +
        )
84 94
85 95
    def mean(self, args: _InputType) -> _OutputType:
86 96
        return self._meanfun(args)
@@ -88,6 +98,11 @@
Loading
88 98
    def cov(self, args0: _InputType, args1: Optional[_InputType] = None) -> _OutputType:
89 99
        return self._covfun(args0, args1)
90 100
101 +
    def covmatrix(
102 +
        self, args0: _InputType, args1: Optional[_InputType] = None
103 +
    ) -> _OutputType:
104 +
        return self._covfun.matrix(args0, args1)
105 +
91 106
    def _sample_at_input(
92 107
        self,
93 108
        rng: np.random.Generator,

@@ -0,0 +1,15 @@
Loading
1 +
"""Fixtures and configuration for doctests."""
2 +
3 +
import numpy as np
4 +
import pytest
5 +
6 +
import probnum as pn
7 +
8 +
9 +
@pytest.fixture(autouse=True)
10 +
def autoimport_packages(doctest_namespace):
11 +
    """This fixture 'imports' standard packages automatically in order to avoid
12 +
    boilerplate code in doctests"""
13 +
14 +
    doctest_namespace["pn"] = pn
15 +
    doctest_namespace["np"] = np

@@ -4,42 +4,35 @@
Loading
4 4
5 5
import numpy as np
6 6
7 -
import probnum.utils as _utils
7 +
from probnum import utils as _utils
8 8
from probnum.typing import IntArgType, ScalarArgType
9 9
10 10
from ._kernel import Kernel
11 11
12 -
_InputType = np.ndarray
13 12
13 +
class WhiteNoise(Kernel):
14 +
    r"""White noise kernel.
14 15
15 -
class WhiteNoise(Kernel[_InputType]):
16 -
    """White noise kernel.
16 +
    Kernel representing independent and identically distributed white noise
17 17
18 -
    Kernel representing independent and identically distributed white noise :math:`k(
19 -
    x_0, x_1) = \\sigma^2 \\delta(x_0, x_1)`.
18 +
    .. math ::
19 +
        k(x_0, x_1) = \sigma^2 \delta(x_0, x_1).
20 20
21 21
    Parameters
22 22
    ----------
23 23
    input_dim :
24 24
        Input dimension of the kernel.
25 25
    sigma :
26 -
        Noise level.
26 +
        Noise level :math:`\sigma`.
27 27
    """
28 28
29 29
    def __init__(self, input_dim: IntArgType, sigma: ScalarArgType = 1.0):
30 30
        self.sigma = _utils.as_numpy_scalar(sigma)
31 -
        super().__init__(input_dim=input_dim, output_dim=1)
31 +
        self._sigma_sq = self.sigma ** 2
32 +
        super().__init__(input_dim=input_dim)
32 33
33 -
    def __call__(self, x0: _InputType, x1: Optional[_InputType] = None) -> np.ndarray:
34 -
35 -
        x0, x1, kernshape = self._check_and_reshape_inputs(x0, x1)
36 -
37 -
        # Compute kernel matrix
34 +
    def _evaluate(self, x0: np.ndarray, x1: Optional[np.ndarray]) -> np.ndarray:
38 35
        if x1 is None:
39 -
            kernmat = self.sigma ** 2 * np.eye(x0.shape[0])
40 -
        else:
41 -
            kernmat = self.sigma ** 2 * np.equal(x0, x1[:, np.newaxis, :]).all(
42 -
                axis=2
43 -
            ).T.astype(float)
36 +
            return np.full_like(x0[..., 0], self._sigma_sq)
44 37
45 -
        return Kernel._reshape_kernelmatrix(kernmat, newshape=kernshape)
38 +
        return self._sigma_sq * np.all(x0 == x1, axis=-1)

@@ -91,7 +91,7 @@
Loading
91 91
        # compute integral mean and variance
92 92
        # Define kernel embedding
93 93
        kernel_embedding = KernelEmbedding(self.kernel, measure)
94 -
        gram = self.kernel(nodes, nodes)
94 +
        gram = self.kernel.matrix(nodes)
95 95
        kernel_mean = kernel_embedding.kernel_mean(nodes)
96 96
        initial_error = kernel_embedding.kernel_variance()
97 97

@@ -9,13 +9,14 @@
Loading
9 9
10 10
from ._kernel import Kernel
11 11
12 -
_InputType = np.ndarray
13 12
13 +
class Linear(Kernel):
14 +
    r"""Linear kernel.
14 15
15 -
class Linear(Kernel[_InputType]):
16 -
    """Linear kernel.
16 +
    Linear covariance function defined by
17 17
18 -
    Linear covariance function defined by :math:`k(x_0, x_1) = x_0^\\top x_1 + c`.
18 +
    .. math ::
19 +
        k(x_0, x_1) = x_0^\top x_1 + c.
19 20
20 21
    Parameters
21 22
    ----------
@@ -33,22 +34,15 @@
Loading
33 34
    >>> import numpy as np
34 35
    >>> from probnum.kernels import Linear
35 36
    >>> K = Linear(input_dim=2)
36 -
    >>> K(np.array([[1, 2], [2, 3]]))
37 +
    >>> xs = np.array([[1, 2], [2, 3]])
38 +
    >>> K.matrix(xs)
37 39
    array([[ 5.,  8.],
38 40
           [ 8., 13.]])
39 41
    """
40 42
41 43
    def __init__(self, input_dim: IntArgType, constant: ScalarArgType = 0.0):
42 44
        self.constant = _utils.as_numpy_scalar(constant)
43 -
        super().__init__(input_dim=input_dim, output_dim=1)
45 +
        super().__init__(input_dim=input_dim)
44 46
45 -
    def __call__(self, x0: _InputType, x1: Optional[_InputType] = None) -> np.ndarray:
46 -
47 -
        x0, x1, kernshape = self._check_and_reshape_inputs(x0, x1)
48 -
49 -
        # Compute kernel matrix
50 -
        if x1 is None:
51 -
            x1 = x0
52 -
        kernmat = x0 @ x1.T + self.constant
53 -
54 -
        return Kernel._reshape_kernelmatrix(kernmat, newshape=kernshape)
47 +
    def _evaluate(self, x0: np.ndarray, x1: Optional[np.ndarray]) -> np.ndarray:
48 +
        return self._euclidean_inner_products(x0, x1) + self.constant
Files Coverage
src/probnum 87.64%
Project Totals (161 files) 87.64%
1
coverage:
2
  precision: 2
3
  status:
4
    project:
5
      default:
6
        target: auto
7
        threshold: 1%
8
    patch:
9
      default:
10
        target: 90%
11
        threshold: 1%
12

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