Kernels

UML Overview

AbstractKernel

Bases: object

Source code in qmcpy/kernel/abstract_kernel.py

def __init__(self, d, torchify, device, compile_call, comiple_call_kwargs):
    super().__init__()
    # dimension
    assert d % 1 == 0 and d > 0, "dimension d must be a positive int"
    self.d = d
    # torchify
    self.torchify = torchify
    if self.torchify:
        import torch

        self.npt = torch
        self.nptarray = torch.tensor
        self.nptarraytype = torch.Tensor
        self.device = torch.device(device)
        self.nptkwargs = {"device": device}
    else:
        self.npt = np
        self.nptarray = np.array
        self.nptarraytype = np.ndarray
        self.device = None
        self.nptkwargs = {}
    self.batch_param_names = []
    if compile_call:
        assert self.torchify, "compile_call requires torchify is True"
        import torch

        self.compiled_parsed___call__ = torch.compile(
            self.parsed___call__, **comiple_call_kwargs
        )
    else:
        self.compiled_parsed___call__ = self.parsed___call__

__call__

__call__(x0, x1, beta0=None, beta1=None, c=None, **kwargs)

Evaluate the kernel with (optional) partial derivatives

\[\sum_{\ell=1}^p c_{\ell} \partial_{\boldsymbol{x}_0}^{\boldsymbol{\beta}_{\ell 0}} \partial_{\boldsymbol{x}_1}^{\boldsymbol{\beta}_{\ell 1}} K(\boldsymbol{x}_0,\boldsymbol{x}_1).\]

Parameters:

Name	Type	Description	Default
x0	Union[ndarray, Tensor]	Shape x0.shape=(...,d) first input to kernel with	required
x1	Union[ndarray, Tensor]	Shape x1.shape=(...,d) second input to kernel with	required
beta0	Union[ndarray, Tensor]	Shape beta0.shape=(p,d) derivative orders with respect to first inputs, \(\boldsymbol{\beta}_0\).	None
beta1	Union[ndarray, Tensor]	Shape beta1.shape=(p,d) derivative orders with respect to first inputs, \(\boldsymbol{\beta}_1\).	None
c	Union[ndarray, Tensor]	Shape c.shape=(p,) coefficients of derivatives.	None
kwargs	dict	keyword arguments to parsed call	{}

Returns: k (Union[np.ndarray, torch.Tensor]): Shape y.shape=(x0+x1).shape[:-1] kernel evaluations.

Source code in qmcpy/kernel/abstract_kernel.py

def __call__(self, x0, x1, beta0=None, beta1=None, c=None, **kwargs):
    r"""
    Evaluate the kernel with (optional) partial derivatives

    $$\sum_{\ell=1}^p c_{\ell} \partial_{\boldsymbol{x}_0}^{\boldsymbol{\beta}_{\ell 0}} \partial_{\boldsymbol{x}_1}^{\boldsymbol{\beta}_{\ell 1}} K(\boldsymbol{x}_0,\boldsymbol{x}_1).$$

    Args:
        x0 (Union[np.ndarray, torch.Tensor]): Shape `x0.shape=(...,d)` first input to kernel with
        x1 (Union[np.ndarray, torch.Tensor]): Shape `x1.shape=(...,d)` second input to kernel with
        beta0 (Union[np.ndarray, torch.Tensor]): Shape `beta0.shape=(p,d)` derivative orders with respect to first inputs, $\boldsymbol{\beta}_0$.
        beta1 (Union[np.ndarray, torch.Tensor]): Shape `beta1.shape=(p,d)` derivative orders with respect to first inputs, $\boldsymbol{\beta}_1$.
        c (Union[np.ndarray, torch.Tensor]): Shape `c.shape=(p,)` coefficients of derivatives.
        kwargs (dict): keyword arguments to parsed call
    Returns:
        k (Union[np.ndarray, torch.Tensor]): Shape `y.shape=(x0+x1).shape[:-1]` kernel evaluations.
    """
    assert isinstance(x0, self.nptarraytype)
    assert isinstance(x0, self.nptarraytype)
    assert (
        x0.shape[-1] == self.d
    ), "the size of the last dimension of x0 must equal d=%d, got x0.shape=%s" % (
        self.d,
        str(tuple(x0.shape)),
    )
    assert (
        x1.shape[-1] == self.d
    ), "the size of the last dimension of x1 must equal d=%d, got x1.shape=%s" % (
        self.d,
        str(tuple(x1.shape)),
    )
    if beta0 is None:
        beta0 = self.npt.zeros((1, self.d), dtype=int, **self.nptkwargs)
    if beta1 is None:
        beta1 = self.npt.zeros((1, self.d), dtype=int, **self.nptkwargs)
    if not isinstance(beta0, self.nptarraytype):
        beta0 = self.nptarray(beta0)
    if not isinstance(beta1, self.nptarraytype):
        beta1 = self.nptarray(beta1)
    beta0 = self.npt.atleast_2d(beta0)
    beta1 = self.npt.atleast_2d(beta1)
    assert (
        beta0.ndim == 2 and beta1.ndim == 2
    ), "beta0 and beta1 must both be 2 dimensional"
    p = beta0.shape[0]
    assert beta0.shape == (
        p,
        self.d,
    ), "expected beta0.shape=(%d,%d) but got beta0.shape=%s" % (
        p,
        self.d,
        str(tuple(beta0.shape)),
    )
    assert beta1.shape == (
        p,
        self.d,
    ), "expected beta1.shape=(%d,%d) but got beta1.shape=%s" % (
        p,
        self.d,
        str(tuple(beta1.shape)),
    )
    assert (beta0 % 1 == 0).all() and (beta0 >= 0).all(), "require int beta0 >= 0"
    assert (beta1 % 1 == 0).all() and (beta1 >= 0).all(), "require int beta1 >= 0"
    if c is None:
        c = self.npt.ones(p, **self.nptkwargs)
    if not isinstance(c, self.nptarraytype):
        c = self.nptarray(c)
    c = self.npt.atleast_1d(c)
    assert c.shape == (p,), "expected c.shape=(%d,) but got c.shape=%s" % (
        p,
        str(tuple(c.shape)),
    )
    if not self.AUTOGRADKERNEL:
        batch_params = self.get_batch_params(max(x0.ndim - 1, x1.ndim - 1))
        k = self.compiled_parsed___call__(
            x0, x1, beta0, beta1, c, batch_params, **kwargs
        )
    else:
        if (beta0 == 0).all() and (beta1 == 0).all():
            batch_params = self.get_batch_params(max(x0.ndim - 1, x1.ndim - 1))
            k = c.sum() * self.compiled_parsed___call__(
                x0, x1, batch_params, **kwargs
            )
        else:  # requires autograd, so self.npt=torch
            assert self.torchify, "autograd requires torchify=True"
            import torch

            incoming_grad_enabled = torch.is_grad_enabled()
            torch.set_grad_enabled(True)
            incoming_grad_enabled_params = {
                pname: param.requires_grad
                for pname, param in self.named_parameters()
            }
            if not incoming_grad_enabled:
                for pname, param in self.named_parameters():
                    param.requires_grad_(False)
            if (beta0 > 0).any():
                tileshapex0 = tuple(
                    self.npt.ceil(
                        self.npt.tensor(x1.shape[:-1])
                        / self.npt.tensor(x0.shape[:-1])
                    ).to(int)
                )
                x0gs = [
                    self.npt.tile(
                        x0[..., j].clone().requires_grad_(True), tileshapex0
                    )
                    for j in range(self.d)
                ]
                [x0gj.requires_grad_(True) for x0gj in x0gs]
                x0g = self.npt.stack(x0gs, dim=-1)
            else:
                x0g = x0
            if (beta1 > 0).any():
                tileshapex1 = tuple(
                    self.npt.ceil(
                        self.npt.tensor(x0.shape[:-1])
                        / self.npt.tensor(x1.shape[:-1])
                    ).to(int)
                )
                x1gs = [
                    self.npt.tile(
                        x1[..., j].clone().requires_grad_(True), tileshapex1
                    )
                    for j in range(self.d)
                ]
                [x1gj.requires_grad_(True) for x1gj in x1gs]
                x1g = self.npt.stack(x1gs, dim=-1)
            else:
                x1g = x1
            batch_params = self.get_batch_params(max(x0.ndim - 1, x1.ndim - 1))
            k = 0.0
            k_base = self.compiled_parsed___call__(x0g, x1g, batch_params, **kwargs)
            for l in range(p):
                if (beta0[l] > 0).any() or (beta1[l] > 0).any():
                    k_part = k_base.clone()
                    for j0 in range(self.d):
                        for _ in range(beta0[l, j0]):
                            k_part = torch.autograd.grad(
                                k_part,
                                x0gs[j0],
                                grad_outputs=torch.ones_like(
                                    k_part, requires_grad=True
                                ),
                                create_graph=True,
                            )[0]
                    for j1 in range(self.d):
                        for _ in range(beta1[l, j1]):
                            k_part = torch.autograd.grad(
                                k_part,
                                x1gs[j1],
                                grad_outputs=torch.ones_like(
                                    k_part, requires_grad=True
                                ),
                                create_graph=True,
                            )[0]
                else:
                    k_part = k_base
                k += c[l] * k_part
            if not incoming_grad_enabled:
                for pname, param in self.named_parameters():
                    param.requires_grad_(incoming_grad_enabled_params[pname])
            if (not incoming_grad_enabled) or (
                (not any(incoming_grad_enabled_params.values()))
                and (not x0.requires_grad)
                and (not x1.requires_grad)
            ):
                k = k.detach()
            torch.set_grad_enabled(incoming_grad_enabled)
    return k

single_integral_01d

single_integral_01d(x)

Evaluate the integral of the kernel over the unit cube

\[\tilde{K}(\boldsymbol{x}) = \int_{[0,1]^d} K(\boldsymbol{x},\boldsymbol{z}) \; \mathrm{d} \boldsymbol{z}.\]

Parameters:

Name	Type	Description	Default
x	Union[ndarray, Tensor]	Shape x0.shape=(...,d) first input to kernel with	required

Returns:

Name	Type	Description
tildek	Union[ndarray, Tensor]	Shape y.shape=x.shape[:-1] integral kernel evaluations.

Source code in qmcpy/kernel/abstract_kernel.py

def single_integral_01d(self, x):
    r"""
    Evaluate the integral of the kernel over the unit cube

    $$\tilde{K}(\boldsymbol{x}) = \int_{[0,1]^d} K(\boldsymbol{x},\boldsymbol{z}) \; \mathrm{d} \boldsymbol{z}.$$

    Args:
        x (Union[np.ndarray, torch.Tensor]): Shape `x0.shape=(...,d)` first input to kernel with

    Returns:
        tildek (Union[np.ndarray, torch.Tensor]): Shape `y.shape=x.shape[:-1]` integral kernel evaluations.
    """
    if self.npt == np:
        assert isinstance(x, np.ndarray)
    else:  # self.npt==torch
        assert isinstance(x, self.npt.Tensor)
    assert (
        x.shape[-1] == self.d
    ), "the size of the last dimension of x must equal d=%d, got x.shape=%s" % (
        self.d,
        str(tuple(x.shape)),
    )
    batch_params = self.get_batch_params(x.ndim - 1)
    return self.parsed_single_integral_01d(x, batch_params)

double_integral_01d

double_integral_01d()

Evaluate the integral of the kernel over the unit cube

\[\tilde{K} = \int_{[0,1]^d} \int_{[0,1]^d} K(\boldsymbol{x},\boldsymbol{z}) \; \mathrm{d} \boldsymbol{x} \; \mathrm{d} \boldsymbol{z}.\]

Returns:

Name	Type	Description
tildek	Union[ndarray, Tensor]	Double integral kernel evaluations.

Source code in qmcpy/kernel/abstract_kernel.py

def double_integral_01d(self):
    r"""
    Evaluate the integral of the kernel over the unit cube

    $$\tilde{K} = \int_{[0,1]^d} \int_{[0,1]^d} K(\boldsymbol{x},\boldsymbol{z}) \; \mathrm{d} \boldsymbol{x} \; \mathrm{d} \boldsymbol{z}.$$

    Returns:
        tildek (Union[np.ndarray, torch.Tensor]): Double integral kernel evaluations.
    """
    raise MethodImplementationError(self, "double_integral_01d")

AbstractKernelScaleLengthscales

Bases: AbstractKernel

Parameters:

Name	Type	Description	Default
d	int	Dimension.	required
scale	Union[ndarray, Tensor]	Scaling factor \(S\).	1.0
lengthscales	Union[ndarray, Tensor]	Lengthscales \(\boldsymbol{\gamma}\).	1.0
shape_scale	list	Shape of scale when np.isscalar(scale).	[1]
shape_lengthscales	list	Shape of lengthscales when np.isscalar(lengthscales)	None
tfs_scale	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_exp_eps_inv, tf_exp_eps)
tfs_lengthscales	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_exp_eps_inv, tf_exp_eps)
torchify	bool	If True, use the torch backend. Set to True if computing gradients with respect to inputs and/or hyperparameters.	False
requires_grad_scale	bool	If True and torchify, set requires_grad=True for scale.	True
requires_grad_lengthscales	bool	If True and torchify, set requires_grad=True for lengthscales.	True
device	device	If torchify, put things onto this device.	'cpu'
compile_call	bool	If True, torch.compile the parsed___call__ method.	False
comiple_call_kwargs	dict	When compile_call is True, pass these keyword arguments to torch.compile.	{}

Source code in qmcpy/kernel/abstract_kernel.py

def __init__(
    self,
    d,
    scale=1.0,
    lengthscales=1.0,
    shape_scale=[1],
    shape_lengthscales=None,
    tfs_scale=(tf_exp_eps_inv, tf_exp_eps),
    tfs_lengthscales=(tf_exp_eps_inv, tf_exp_eps),
    torchify=False,
    requires_grad_scale=True,
    requires_grad_lengthscales=True,
    device="cpu",
    compile_call=False,
    comiple_call_kwargs={},
):
    r"""
    Args:
        d (int): Dimension.
        scale (Union[np.ndarray, torch.Tensor]): Scaling factor $S$.
        lengthscales (Union[np.ndarray, torch.Tensor]): Lengthscales $\boldsymbol{\gamma}$.
        shape_scale (list): Shape of `scale` when `np.isscalar(scale)`.
        shape_lengthscales (list): Shape of `lengthscales` when `np.isscalar(lengthscales)`
        tfs_scale (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_lengthscales (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        torchify (bool): If `True`, use the `torch` backend. Set to `True` if computing gradients with respect to inputs and/or hyperparameters.
        requires_grad_scale (bool): If `True` and `torchify`, set `requires_grad=True` for `scale`.
        requires_grad_lengthscales (bool): If `True` and `torchify`, set `requires_grad=True` for `lengthscales`.
        device (torch.device): If `torchify`, put things onto this device.
        compile_call (bool): If `True`, `torch.compile` the `parsed___call__` method.
        comiple_call_kwargs (dict): When `compile_call` is `True`, pass these keyword arguments to `torch.compile`.
    """
    super().__init__(
        d=d,
        torchify=torchify,
        device=device,
        compile_call=compile_call,
        comiple_call_kwargs=comiple_call_kwargs,
    )
    self.raw_scale = self.parse_assign_param(
        pname="scale",
        param=scale,
        shape_param=shape_scale,
        requires_grad_param=requires_grad_scale,
        tfs_param=tfs_scale,
        endsize_ops=[1],
        constraints=["POSITIVE"],
    )
    self.tfs_scale = tfs_scale
    self.batch_param_names.append("scale")
    self.raw_lengthscales = self.parse_assign_param(
        pname="lengthscales",
        param=lengthscales,
        shape_param=[self.d] if shape_lengthscales is None else shape_lengthscales,
        requires_grad_param=requires_grad_lengthscales,
        tfs_param=tfs_lengthscales,
        endsize_ops=[1, self.d],
        constraints=["POSITIVE"],
    )
    self.tfs_lengthscales = tfs_lengthscales
    self.batch_param_names.append("lengthscales")

KernelShiftInvar

Bases: AbstractSIDSIKernel

Shift invariant kernel with smoothness \(\boldsymbol{\alpha}\), product weights (lengthscales) \(\boldsymbol{\gamma}\), and scale \(S\):

\[\begin{aligned} K(\boldsymbol{x},\boldsymbol{z}) &= S \prod_{j=1}^d \left(1+ \gamma_j \tilde{K}_{\alpha_j}((x_j - z_j) \mod 1))\right), \\ \tilde{K}_\alpha(x) &= (-1)^{\alpha+1}\frac{(2 \pi)^{2 \alpha}}{(2\alpha)!} B_{2\alpha}(x) \end{aligned}\]

where \(B_n\) is the \(n^\text{th}\) Bernoulli polynomial.

Examples:

>>> from qmcpy import Lattice, fftbr, ifftbr
>>> n = 8
>>> d = 4
>>> lat = Lattice(d,seed=11)
>>> x = lat(n)
>>> x.shape
(8, 4)
>>> x.dtype
dtype('float64')
>>> kernel = KernelShiftInvar(
...     d = d, 
...     alpha = list(range(1,d+1)),
...     scale = 10,
...     lengthscales = [1/j**2 for j in range(1,d+1)])
>>> k00 = kernel(x[0],x[0])
>>> k00.item()
91.23444453396341
>>> k0 = kernel(x,x[0])
>>> with np.printoptions(precision=2):
...     print(k0)
[91.23 -2.32  5.69  5.69 12.7  -4.78 -4.78 12.7 ]
>>> assert k0[0]==k00
>>> kmat = kernel(x[:,None,:],x[None,:,:])
>>> with np.printoptions(precision=2):
...     print(kmat)
[[91.23 -2.32  5.69  5.69 12.7  -4.78 -4.78 12.7 ]
 [-2.32 91.23  5.69  5.69 -4.78 12.7  12.7  -4.78]
 [ 5.69  5.69 91.23 -2.32 12.7  -4.78 12.7  -4.78]
 [ 5.69  5.69 -2.32 91.23 -4.78 12.7  -4.78 12.7 ]
 [12.7  -4.78 12.7  -4.78 91.23 -2.32  5.69  5.69]
 [-4.78 12.7  -4.78 12.7  -2.32 91.23  5.69  5.69]
 [-4.78 12.7  12.7  -4.78  5.69  5.69 91.23 -2.32]
 [12.7  -4.78 -4.78 12.7   5.69  5.69 -2.32 91.23]]
>>> assert (kmat[:,0]==k0).all()
>>> lam = np.sqrt(n)*fftbr(k0)
>>> y = np.random.Generator(np.random.PCG64(7)).uniform(low=0,high=1,size=(n))
>>> np.allclose(ifftbr(fftbr(y)*lam),kmat@y)
True
>>> np.allclose(ifftbr(fftbr(y)/lam),np.linalg.solve(kmat,y))
True
>>> import torch 
>>> xtorch = torch.from_numpy(x)
>>> kernel_torch = KernelShiftInvar(
...     d = d, 
...     alpha = list(range(1,d+1)),
...     scale = 10,
...     lengthscales = [1/j**2 for j in range(1,d+1)],
...     torchify = True)
>>> kmat_torch = kernel_torch(xtorch[:,None,:],xtorch[None,:,:])
>>> np.allclose(kmat_torch.detach().numpy(),kmat)
True
>>> kernel.single_integral_01d(x)
array([10., 10., 10., 10., 10., 10., 10., 10.])
>>> kernel_torch.single_integral_01d(xtorch)
tensor([10., 10., 10., 10., 10., 10., 10., 10.], dtype=torch.float64,
       grad_fn=<AddBackward0>)

Batch Params

>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kernel = KernelShiftInvar(
...     d = 2, 
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,2])
>>> x = rng.uniform(low=0,high=1,size=(6,5,2))
>>> kernel(x,x).shape 
(4, 3, 6, 5)
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
(4, 3, 6, 5, 5)
>>> kfast = kernel(x[:,None,:,None,:],x[None,:,None,:,:])
>>> kfast.shape
(4, 3, 6, 6, 5, 5)
>>> kstable = kernel(x[:,None,:,None,:],x[None,:,None,:,:],stable=True)
>>> np.abs(kfast-kstable).max()
np.float64(4.440892098500626e-16)

Derivatives

>>> rng = np.random.Generator(np.random.PCG64(7))
>>> scale = rng.uniform(low=0,high=1,size=(1,))
>>> lengthscales = rng.uniform(low=0,high=1,size=(3,))
>>> kernel = KernelShiftInvar(
...     d = 3,
...     alpha = 3,
...     torchify = True,
...     scale = torch.from_numpy(scale),
...     lengthscales = torch.from_numpy(lengthscales))
>>> x0 = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,))).requires_grad_(True)
>>> x1 = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,))).requires_grad_(True)
>>> x2 = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,))).requires_grad_(True)
>>> x = torch.stack([x0,x1,x2],axis=-1)
>>> z0 = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,))).requires_grad_(True)
>>> z1 = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,))).requires_grad_(True)
>>> z2 = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,))).requires_grad_(True)
>>> z = torch.stack([z0,z1,z2],axis=-1)
>>> c = torch.from_numpy(rng.uniform(low=0,high=1,size=(2,)))
>>> beta0 = torch.tensor([
...     [1,0,0],
...     [0,2,0]])
>>> beta1 = torch.tensor([
...     [0,0,2],
...     [2,1,0]])
>>> with torch.no_grad():
...     y = kernel(x,z,beta0,beta1,c)
>>> with np.printoptions(formatter={"float": lambda x: "%.2f"%x}):
...     y.numpy()
array([1455.14, 9475.57, 7807.08, 2785.47])
>>> y_no_deriv = kernel(x,z)
>>> y_first = y_no_deriv.clone()
>>> y_first = torch.autograd.grad(y_first,x0,grad_outputs=torch.ones_like(y_first,requires_grad=True),create_graph=True)[0]
>>> y_first = torch.autograd.grad(y_first,z2,grad_outputs=torch.ones_like(y_first,requires_grad=True),create_graph=True)[0]
>>> y_first = torch.autograd.grad(y_first,z2,grad_outputs=torch.ones_like(y_first,requires_grad=True),create_graph=True)[0]
>>> y_second = y_no_deriv.clone()
>>> y_second = torch.autograd.grad(y_second,x1,grad_outputs=torch.ones_like(y_second,requires_grad=True),create_graph=True)[0]
>>> y_second = torch.autograd.grad(y_second,x1,grad_outputs=torch.ones_like(y_second,requires_grad=True),create_graph=True)[0]
>>> y_second = torch.autograd.grad(y_second,z0,grad_outputs=torch.ones_like(y_second,requires_grad=True),create_graph=True)[0]
>>> y_second = torch.autograd.grad(y_second,z0,grad_outputs=torch.ones_like(y_second,requires_grad=True),create_graph=True)[0]
>>> y_second = torch.autograd.grad(y_second,z1,grad_outputs=torch.ones_like(y_second,requires_grad=True),create_graph=True)[0]
>>> yhat = (y_first*c[0]+y_second*c[1]).detach()
>>> with np.printoptions(formatter={"float": lambda x: "%.3f"%x}):
...     yhat.numpy()
array([1455.140, 9475.570, 7807.076, 2785.473])
>>> torch.allclose(y,yhat)
True
>>> kernel = KernelShiftInvar(
...     d = 3,
...     alpha = 3,
...     scale = scale,
...     lengthscales = lengthscales)
>>> ynp = kernel(x.detach().numpy(),z.detach().numpy(),beta0.numpy(),beta1.numpy(),c.numpy())
>>> with np.printoptions(formatter={"float": lambda x: "%.2f"%x}):
...     ynp
array([1455.14, 9475.59, 7807.09, 2785.48])
>>> np.allclose(ynp,y.numpy())
True

References:

1. Kaarnioja, Vesa, Frances Y. Kuo, and Ian H. Sloan.
   "Lattice-based kernel approximation and serendipitous weights for parametric PDEs in very high dimensions."
   International Conference on Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing. Cham: Springer International Publishing, 2022.

Parameters:

Name	Type	Description	Default
d	int	Dimension.	required
scale	Union[ndarray, Tensor]	Scaling factor \(S\).	1.0
lengthscales	Union[ndarray, Tensor]	Product weights \((\gamma_1,\dots,\gamma_d)\).	None
alpha	Union[ndarray, Tensor]	Smoothness parameters \((\alpha_1,\dots,\alpha_d)\) where \(\alpha_j \geq 1\) for \(j=1,\dots,d\).	2
shape_scale	list	Shape of scale when np.isscalar(scale).	[1]
shape_lengthscales	list	Shape of lengthscales when np.isscalar(lengthscales)	None
tfs_scale	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_exp_eps_inv, tf_exp_eps)
tfs_lengthscales	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_exp_eps_inv, tf_exp_eps)
torchify	bool	If True, use the torch backend. Set to True if computing gradients with respect to inputs and/or hyperparameters.	False
requires_grad_scale	bool	If True and torchify, set requires_grad=True for scale.	True
requires_grad_lengthscales	bool	If True and torchify, set requires_grad=True for lengthscales.	True
device	device	If torchify, put things onto this device.	'cpu'
compile_call	bool	If True, torch.compile the parsed___call__ method.	False
comiple_call_kwargs	dict	When compile_call is True, pass these keyword arguments to torch.compile.	{}

Source code in qmcpy/kernel/si_dsi_kernels.py

def __init__(
    self,
    d,
    scale=1.0,
    lengthscales=None,
    alpha=2,
    shape_scale=[1],
    shape_lengthscales=None,
    tfs_scale=(tf_exp_eps_inv, tf_exp_eps),
    tfs_lengthscales=(tf_exp_eps_inv, tf_exp_eps),
    torchify=False,
    requires_grad_scale=True,
    requires_grad_lengthscales=True,
    device="cpu",
    compile_call=False,
    comiple_call_kwargs={},
):
    r"""
    Args:
        d (int): Dimension.
        scale (Union[np.ndarray, torch.Tensor]): Scaling factor $S$.
        lengthscales (Union[np.ndarray, torch.Tensor]): Product weights $(\gamma_1,\dots,\gamma_d)$.
        alpha (Union[np.ndarray, torch.Tensor]): Smoothness parameters $(\alpha_1,\dots,\alpha_d)$ where $\alpha_j \geq 1$ for $j=1,\dots,d$.
        shape_scale (list): Shape of `scale` when `np.isscalar(scale)`.
        shape_lengthscales (list): Shape of `lengthscales` when `np.isscalar(lengthscales)`
        tfs_scale (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_lengthscales (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        torchify (bool): If `True`, use the `torch` backend. Set to `True` if computing gradients with respect to inputs and/or hyperparameters.
        requires_grad_scale (bool): If `True` and `torchify`, set `requires_grad=True` for `scale`.
        requires_grad_lengthscales (bool): If `True` and `torchify`, set `requires_grad=True` for `lengthscales`.
        device (torch.device): If `torchify`, put things onto this device.
        compile_call (bool): If `True`, `torch.compile` the `parsed___call__` method.
        comiple_call_kwargs (dict): When `compile_call` is `True`, pass these keyword arguments to `torch.compile`.
    """
    super().__init__(
        d=d,
        scale=scale,
        lengthscales=lengthscales,
        alpha=alpha,
        shape_alpha=[d],
        alpha_endsize_ops=[d],
        shape_scale=shape_scale,
        shape_lengthscales=shape_lengthscales,
        tfs_alpha=(tf_identity, tf_identity),
        tfs_scale=tfs_scale,
        requires_grad_alpha=False,
        tfs_lengthscales=tfs_lengthscales,
        torchify=torchify,
        requires_grad_scale=requires_grad_scale,
        requires_grad_lengthscales=requires_grad_lengthscales,
        device=device,
        compile_call=compile_call,
        comiple_call_kwargs=comiple_call_kwargs,
    )
    assert self.alpha.shape == (self.d,)
    assert all(int(alphaj) in BERNOULLIPOLYSDICT for alphaj in self.alpha)
    if self.torchify:
        import torch

        self.lgamma = torch.lgamma
    else:
        self.lgamma = scipy.special.loggamma

KernelShiftInvarCombined

Bases: AbstractSIDSIKernel

Shift invariant kernel with combination weights \(\boldsymbol{\alpha}_1,\dots,\boldsymbol{\alpha}_d \in \mathbb{R}_{>0}^4\), product weights (lengthscales) \(\boldsymbol{\gamma}\), and scale \(S\):

\[\begin{aligned} K(\boldsymbol{x},\boldsymbol{z}) &= S \prod_{j=1}^d \left(1+ \gamma_j \left(\sum_{p=1}^4 \alpha_{jp} \tilde{K}_p(x_j \mod 1 z_j)\right)\right) \end{aligned}\]

where, \(\tilde{K}_p\) are defined in KernelShiftInvar for \(p \in \{1,2,3,4\}\)

Examples:

>>> from qmcpy import Lattice, fftbr, ifftbr
>>> n = 8
>>> d = 4
>>> lat = Lattice(d,seed=11)
>>> x = lat(n)
>>> x.shape
(8, 4)
>>> x.dtype
dtype('float64')
>>> kernel = KernelShiftInvarCombined(
...     d = d,
...     scale = 10,
...     lengthscales = [1/j**2 for j in range(1,d+1)])
>>> k00 = kernel(x[0],x[0])
>>> k00.item()
117.22427096475315
>>> k0 = kernel(x,x[0])
>>> with np.printoptions(precision=2):
...     print(k0)
[ 117.22  153.84   28.59   36.1  -139.83   77.1    -8.51   35.22]
>>> assert k0[0]==k00
>>> kmat = kernel(x[:,None,:],x[None,:,:])
>>> with np.printoptions(precision=2):
...     print(kmat)
[[ 117.22  153.84   36.1    28.59   35.22   -8.51   77.1  -139.83]
 [ 153.84  117.22   28.59   36.1    -8.51   35.22 -139.83   77.1 ]
 [  28.59   36.1   117.22  153.84 -139.83   77.1    35.22   -8.51]
 [  36.1    28.59  153.84  117.22   77.1  -139.83   -8.51   35.22]
 [-139.83   77.1    35.22   -8.51  117.22  153.84   36.1    28.59]
 [  77.1  -139.83   -8.51   35.22  153.84  117.22   28.59   36.1 ]
 [  -8.51   35.22 -139.83   77.1    28.59   36.1   117.22  153.84]
 [  35.22   -8.51   77.1  -139.83   36.1    28.59  153.84  117.22]]
>>> assert (kmat[:,0]==k0).all()
>>> lam = np.sqrt(n)*fftbr(k0)
>>> y = np.random.Generator(np.random.PCG64(7)).uniform(low=0,high=1,size=(n))
>>> np.allclose(ifftbr(fftbr(y)*lam),kmat@y)
True
>>> np.allclose(ifftbr(fftbr(y)/lam),np.linalg.solve(kmat,y))
True
>>> import torch
>>> xtorch = torch.from_numpy(x)
>>> kernel_torch = KernelShiftInvarCombined(
...     d = d,
...     scale = 10,
...     lengthscales = [1/j**2 for j in range(1,d+1)],
...     torchify = True)
>>> kmat_torch = kernel_torch(xtorch[:,None,:],xtorch[None,:,:])
>>> np.allclose(kmat_torch.detach().numpy(),kmat)
True
>>> kernel.single_integral_01d(x)
array([10., 10., 10., 10., 10., 10., 10., 10.])
>>> kernel_torch.single_integral_01d(xtorch)
tensor([10., 10., 10., 10., 10., 10., 10., 10.], dtype=torch.float64,
       grad_fn=<AddBackward0>)

Batch Params

>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kernel = KernelShiftInvarCombined(
...     d = 2,
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,2])
>>> x = rng.uniform(low=0,high=1,size=(6,5,2))
>>> kernel(x,x).shape
(4, 3, 6, 5)
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
(4, 3, 6, 5, 5)
>>> kfast = kernel(x[:,None,:,None,:],x[None,:,None,:,:])
>>> kfast.shape
(4, 3, 6, 6, 5, 5)
>>> kstable = kernel(x[:,None,:,None,:],x[None,:,None,:,:],stable=True)
>>> np.abs(kfast-kstable).max()
np.float64(3.552713678800501e-15)

References:

1. Kaarnioja, Vesa, Frances Y. Kuo, and Ian H. Sloan. "Lattice-based kernel approximation and serendipitous weights for parametric PDEs in very high dimensions." International Conference on Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing. Cham: Springer International Publishing, 2022.

Parameters:

Name	Type	Description	Default
d	int	Dimension.	required
scale	Union[ndarray, Tensor]	Scaling factor \(S\).	1.0
lengthscales	Union[ndarray, Tensor]	Product weights \((\gamma_1,\dots,\gamma_d)\).	None
alpha	Union[ndarray, Tensor]	Weights \(\boldsymbol{\alpha}_1,\dots,\boldsymbol{\alpha}_d \in \mathbb{R}_{>0}^4\).	1
shape_scale	list	Shape of scale when np.isscalar(scale).	[1]
shape_lengthscales	list	Shape of lengthscales when np.isscalar(lengthscales)	None
shape_alpha	list	Shape of alpha when np.isscalar(alpha)	None
tfs_scale	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_exp_eps_inv, tf_exp_eps)
tfs_lengthscales	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_exp_eps_inv, tf_exp_eps)
tfs_alpha	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_exp_eps_inv, tf_exp_eps)
torchify	bool	If True, use the torch backend. Set to True if computing gradients with respect to inputs and/or hyperparameters.	False
requires_grad_scale	bool	If True and torchify, set requires_grad=True for scale.	True
requires_grad_lengthscales	bool	If True and torchify, set requires_grad=True for lengthscales.	True
requires_grad_alpha	bool	If True and torchify, set requires_grad=True for alpha.	True
device	device	If torchify, put things onto this device.	'cpu'
compile_call	bool	If True, torch.compile the parsed___call__ method.	False
comiple_call_kwargs	dict	When compile_call is True, pass these keyword arguments to torch.compile.	{}

Source code in qmcpy/kernel/si_dsi_kernels.py

def __init__(
    self,
    d,
    scale=1.0,
    lengthscales=None,
    alpha=1,
    shape_scale=[1],
    shape_lengthscales=None,
    shape_alpha=None,
    tfs_scale=(tf_exp_eps_inv, tf_exp_eps),
    tfs_lengthscales=(tf_exp_eps_inv, tf_exp_eps),
    tfs_alpha=(tf_exp_eps_inv, tf_exp_eps),
    torchify=False,
    requires_grad_scale=True,
    requires_grad_lengthscales=True,
    requires_grad_alpha=True,
    device="cpu",
    compile_call=False,
    comiple_call_kwargs={},
):
    r"""
    Args:
        d (int): Dimension.
        scale (Union[np.ndarray, torch.Tensor]): Scaling factor $S$.
        lengthscales (Union[np.ndarray, torch.Tensor]): Product weights $(\gamma_1,\dots,\gamma_d)$.
        alpha (Union[np.ndarray, torch.Tensor]): Weights $\boldsymbol{\alpha}_1,\dots,\boldsymbol{\alpha}_d \in \mathbb{R}_{>0}^4$.
        shape_scale (list): Shape of `scale` when `np.isscalar(scale)`.
        shape_lengthscales (list): Shape of `lengthscales` when `np.isscalar(lengthscales)`
        shape_alpha (list): Shape of `alpha` when `np.isscalar(alpha)`
        tfs_scale (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_lengthscales (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_alpha (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        torchify (bool): If `True`, use the `torch` backend. Set to `True` if computing gradients with respect to inputs and/or hyperparameters.
        requires_grad_scale (bool): If `True` and `torchify`, set `requires_grad=True` for `scale`.
        requires_grad_lengthscales (bool): If `True` and `torchify`, set `requires_grad=True` for `lengthscales`.
        requires_grad_alpha (bool): If `True` and `torchify`, set `requires_grad=True` for `alpha`.
        device (torch.device): If `torchify`, put things onto this device.
        compile_call (bool): If `True`, `torch.compile` the `parsed___call__` method.
        comiple_call_kwargs (dict): When `compile_call` is `True`, pass these keyword arguments to `torch.compile`.
    """
    super().__init__(
        d=d,
        scale=scale,
        lengthscales=lengthscales,
        alpha=alpha,
        alpha_endsize_ops=[d],
        shape_scale=shape_scale,
        shape_lengthscales=shape_lengthscales,
        shape_alpha=[4, d] if shape_alpha is None else shape_alpha,
        tfs_scale=tfs_scale,
        tfs_alpha=tfs_alpha,
        tfs_lengthscales=tfs_lengthscales,
        torchify=torchify,
        requires_grad_scale=requires_grad_scale,
        requires_grad_lengthscales=requires_grad_lengthscales,
        requires_grad_alpha=requires_grad_alpha,
        device=device,
        compile_call=compile_call,
        comiple_call_kwargs=comiple_call_kwargs,
    )
    assert self.alpha.shape[-2:] == (4, d)
    if self.torchify:
        import torch

        self.lgamma = torch.lgamma
    else:
        self.lgamma = scipy.special.loggamma
    a = self.npt.arange(1, 5, **self.nptkwargs)
    self.coeffs = (-1) ** (a + 1) * self.npt.exp(
        2 * a * np.log(2 * np.pi) - self.lgamma(2 * a + 1)
    )

KernelDigShiftInvar

Bases: AbstractSIDSIKernel

Digitally shift invariant kernel in base \(b=2\) with smoothness \(\boldsymbol{\alpha}\), product weights \(\boldsymbol{\gamma}\), and scale \(S\):

\[\begin{aligned} K(\boldsymbol{x},\boldsymbol{z}) &= S \prod_{j=1}^d \left(1+ \gamma_j \tilde{K}_{\alpha_j}(x_j \oplus z_j)\right), \qquad\mathrm{where} \\ \tilde{K}_1(x) &= 6 \left(\frac{1}{6} - 2^{\lfloor \log_2(x) \rfloor -1}\right), \\ \tilde{K}_2(x) &= \sum_{k \in \mathbb{N}} \frac{\mathrm{wal}_k(x)}{2^{\mu_2(k)}} = -\beta(x) x + \frac{5}{2}\left[1-t_1(x)\right]-1, \\ \tilde{K}_3(x) &= \sum_{k \in \mathbb{N}} \frac{\mathrm{wal}_k(x)}{2^{\mu_3(k)}} = \beta(x)x^2-5\left[1-t_1(x)\right]x+\frac{43}{18}\left[1-t_2(x)\right]-1, \\ \tilde{K}_4(x) &= \sum_{k \in \mathbb{N}} \frac{\mathrm{wal}_k(x)}{2^{\mu_4(k)}} = - \frac{2}{3}\beta(x)x^3+5\left[1-t_1(x)\right]x^2 - \frac{43}{9}\left[1-t_2(x)\right]x +\frac{701}{294}\left[1-t_3(x)\right]+\beta(x)\left[\frac{1}{48}\sum_{a=0}^\infty \frac{\mathrm{wal}_{2^a}(x)}{2^{3a}} - \frac{1}{42}\right] - 1. \end{aligned}\]

where

• \(x \oplus z\) is XOR between bits,
• \(\mathrm{wal}_k\) is the \(k^\text{th}\) Walsh function,
• \(\beta(x) = - \lfloor \log_2(x) \rfloor\) and \(t_\nu(x) = 2^{-\nu \beta(x)}\) where \(\beta(0)=t_\nu(0) = 0\), and
• and \(\mu_\alpha\) is the Dick weight function which sums the first \(\alpha\) largest indices of \(1\) bits in the binary expansion of \(k\) e.g. \(k=13=1101_2\) has 1-bit indexes \((4,3,1)\) so

\[\mu_1(k) = 4, \mu_2(k) = 4+3, \mu_3(k) = 4+3+1 = \mu_4(k) = \mu_5(k) = \dots.\]

Examples:

>>> from qmcpy import DigitalNetB2, fwht
>>> n = 8
>>> d = 4
>>> dnb2 = DigitalNetB2(d,seed=11)
>>> x = dnb2(n,return_binary=True)
>>> x.shape
(8, 4)
>>> x.dtype
dtype('uint64')
>>> kernel = KernelDigShiftInvar(
...     d = d, 
...     t = dnb2.t,
...     alpha = list(range(1,d+1)),
...     scale = 10,
...     lengthscales = [1/j**2 for j in range(1,d+1)])
>>> k00 = kernel(x[0],x[0])
>>> k00.item()
34.490370029184525
>>> k0 = kernel(x,x[0])
>>> with np.printoptions(precision=2):
...     print(k0)
[34.49  4.15  9.59  4.98 15.42  5.45 11.99  4.51]
>>> assert k0[0]==k00
>>> kmat = kernel(x[:,None,:],x[None,:,:])
>>> with np.printoptions(precision=2):
...     print(kmat)
[[34.49  4.15  9.59  4.98 15.42  5.45 11.99  4.51]
 [ 4.15 34.49  4.98  9.59  5.45 15.42  4.51 11.99]
 [ 9.59  4.98 34.49  4.15 11.99  4.51 15.42  5.45]
 [ 4.98  9.59  4.15 34.49  4.51 11.99  5.45 15.42]
 [15.42  5.45 11.99  4.51 34.49  4.15  9.59  4.98]
 [ 5.45 15.42  4.51 11.99  4.15 34.49  4.98  9.59]
 [11.99  4.51 15.42  5.45  9.59  4.98 34.49  4.15]
 [ 4.51 11.99  5.45 15.42  4.98  9.59  4.15 34.49]]
>>> assert (kmat[:,0]==k0).all()
>>> lam = np.sqrt(n)*fwht(k0)
>>> y = np.random.Generator(np.random.PCG64(7)).uniform(low=0,high=1,size=(n))
>>> np.allclose(fwht(fwht(y)*lam),kmat@y)
True
>>> np.allclose(fwht(fwht(y)/lam),np.linalg.solve(kmat,y))
True
>>> import torch 
>>> xtorch = bin_from_numpy_to_torch(x)
>>> kernel_torch = KernelDigShiftInvar(
...     d = d, 
...     t = dnb2.t,
...     alpha = list(range(1,d+1)),
...     scale = 10,
...     lengthscales = [1/j**2 for j in range(1,d+1)],
...     torchify = True)
>>> kmat_torch = kernel_torch(xtorch[:,None,:],xtorch[None,:,:])
>>> np.allclose(kmat_torch.detach().numpy(),kmat)
True
>>> xf = to_float(x,dnb2.t)
>>> kmat_from_floats = kernel(xf[:,None,:],xf[None,:,:])
>>> np.allclose(kmat,kmat_from_floats)
True
>>> xftorch = to_float(xtorch,dnb2.t)
>>> xftorch.dtype
torch.float32
>>> kmat_torch_from_floats = kernel_torch(xftorch[:,None,:],xftorch[None,:,:])
>>> torch.allclose(kmat_torch_from_floats,kmat_torch)
True
>>> kernel.single_integral_01d(x)
array([10., 10., 10., 10., 10., 10., 10., 10.])
>>> kernel_torch.single_integral_01d(xtorch)
tensor([10., 10., 10., 10., 10., 10., 10., 10.], grad_fn=<AddBackward0>)

Batch Params

>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kernel = KernelDigShiftInvar(
...     d = 2, 
...     t = 10,
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,2])
>>> x = rng.uniform(low=0,high=1,size=(6,5,2))
>>> kernel(x,x).shape 
(4, 3, 6, 5)
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
(4, 3, 6, 5, 5)
>>> kfast = kernel(x[:,None,:,None,:],x[None,:,None,:,:])
>>> kfast.shape
(4, 3, 6, 6, 5, 5)
>>> kstable = kernel(x[:,None,:,None,:],x[None,:,None,:,:],stable=True)
>>> np.abs(kfast-kstable).max()
np.float64(4.440892098500626e-16)

References:

1. Dick, Josef.
   "Walsh spaces containing smooth functions and quasi-Monte Carlo rules of arbitrary high order."
   SIAM Journal on Numerical Analysis 46.3 (2008): 1519-1553.

2. Dick, Josef.
   "The decay of the Walsh coefficients of smooth functions."
   Bulletin of the Australian Mathematical Society 80.3 (2009): 430-453.

3. Jagadeeswaran, Rathinavel, and Fred J. Hickernell.
   "Fast automatic Bayesian cubature using Sobol' sampling."
   Advances in Modeling and Simulation: Festschrift for Pierre L'Ecuyer. Cham: Springer International Publishing, 2022. 301-318.

4. Rathinavel, Jagadeeswaran.
   Fast automatic Bayesian cubature using matching kernels and designs.
   Illinois Institute of Technology, 2019.

5. Sorokin, Aleksei.
   "A Unified Implementation of Quasi-Monte Carlo Generators, Randomization Routines, and Fast Kernel Methods."
   arXiv preprint arXiv:2502.14256 (2025).

Parameters:

Name	Type	Description	Default
d	int	Dimension.	required
t	int	number of bits in binary represtnations. Typically dnb2.t where isinstance(dnb2,DigitalNetB2).	None
scale	Union[ndarray, Tensor]	Scaling factor \(S\).	1.0
lengthscales	Union[ndarray, Tensor]	Product weights \((\gamma_1,\dots,\gamma_d)\).	None
alpha	Union[ndarray, Tensor]	Smoothness parameters \((\alpha_1,\dots,\alpha_d)\) where \(\alpha_j \geq 1\) for \(j=1,\dots,d\).	2
shape_scale	list	Shape of scale when np.isscalar(scale).	[1]
shape_lengthscales	list	Shape of lengthscales when np.isscalar(lengthscales)	None
tfs_scale	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_exp_eps_inv, tf_exp_eps)
tfs_lengthscales	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_exp_eps_inv, tf_exp_eps)
torchify	bool	If True, use the torch backend. Set to True if computing gradients with respect to inputs and/or hyperparameters.	False
requires_grad_scale	bool	If True and torchify, set requires_grad=True for scale.	True
requires_grad_lengthscales	bool	If True and torchify, set requires_grad=True for lengthscales.	True
device	device	If torchify, put things onto this device.	'cpu'
compile_call	bool	If True, torch.compile the parsed___call__ method.	False
comiple_call_kwargs	dict	When compile_call is True, pass these keyword arguments to torch.compile.	{}

Source code in qmcpy/kernel/si_dsi_kernels.py

def __init__(
    self,
    d,
    t=None,
    scale=1.0,
    lengthscales=None,
    alpha=2,
    shape_scale=[1],
    shape_lengthscales=None,
    tfs_scale=(tf_exp_eps_inv, tf_exp_eps),
    tfs_lengthscales=(tf_exp_eps_inv, tf_exp_eps),
    torchify=False,
    requires_grad_scale=True,
    requires_grad_lengthscales=True,
    device="cpu",
    compile_call=False,
    comiple_call_kwargs={},
):
    r"""
    Args:
        d (int): Dimension.
        t (int): number of bits in binary represtnations. Typically `dnb2.t` where `isinstance(dnb2,DigitalNetB2)`.
        scale (Union[np.ndarray, torch.Tensor]): Scaling factor $S$.
        lengthscales (Union[np.ndarray, torch.Tensor]): Product weights $(\gamma_1,\dots,\gamma_d)$.
        alpha (Union[np.ndarray, torch.Tensor]): Smoothness parameters $(\alpha_1,\dots,\alpha_d)$ where $\alpha_j \geq 1$ for $j=1,\dots,d$.
        shape_scale (list): Shape of `scale` when `np.isscalar(scale)`.
        shape_lengthscales (list): Shape of `lengthscales` when `np.isscalar(lengthscales)`
        tfs_scale (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_lengthscales (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        torchify (bool): If `True`, use the `torch` backend. Set to `True` if computing gradients with respect to inputs and/or hyperparameters.
        requires_grad_scale (bool): If `True` and `torchify`, set `requires_grad=True` for `scale`.
        requires_grad_lengthscales (bool): If `True` and `torchify`, set `requires_grad=True` for `lengthscales`.
        device (torch.device): If `torchify`, put things onto this device.
        compile_call (bool): If `True`, `torch.compile` the `parsed___call__` method.
        comiple_call_kwargs (dict): When `compile_call` is `True`, pass these keyword arguments to `torch.compile`.
    """
    super().__init__(
        d=d,
        scale=scale,
        lengthscales=lengthscales,
        alpha=alpha,
        shape_alpha=[d],
        alpha_endsize_ops=[d],
        shape_scale=shape_scale,
        shape_lengthscales=shape_lengthscales,
        tfs_alpha=(tf_identity, tf_identity),
        tfs_scale=tfs_scale,
        tfs_lengthscales=tfs_lengthscales,
        torchify=torchify,
        requires_grad_alpha=False,
        requires_grad_scale=requires_grad_scale,
        requires_grad_lengthscales=requires_grad_lengthscales,
        device=device,
        compile_call=compile_call,
        comiple_call_kwargs=comiple_call_kwargs,
    )
    assert self.alpha.shape == (self.d,)
    self.set_t(t)
    assert all(1 <= int(alphaj) <= 4 for alphaj in self.alpha)

KernelDigShiftInvarAdaptiveAlpha

Bases: AbstractSIDSIKernel

Digitally shift invariant kernel in base \(b=2\) with smoothness \(\boldsymbol{\alpha} \geq \boldsymbol{0}\), product weights \(\boldsymbol{\gamma}\), and scale \(S\):

\[\begin{aligned} K(\boldsymbol{x},\boldsymbol{z}) &= S \prod_{j=1}^d \left(1+ \gamma_j \tilde{K}_{\alpha_j}(x_j \oplus z_j)\right), \qquad\mathrm{where} \\ \tilde{K}_\alpha(x) &= \sum_{k \in \mathbb{N}} \frac{\mathrm{wal}_k(x)}{2^{{\alpha+1} (\mu_1(k)-1)}} = \frac{2^{\alpha+1}}{2^{\alpha+1}-2} - \left(\frac{2^{\alpha+1}}{2^{\alpha+1}-2}+1\right) 2^{\alpha(\lfloor \log_2(x) \rfloor+1)}, \\ \end{aligned}\]

where

• \(x \oplus z\) is XOR between bits,
• \(\mathrm{wal}_k\) is the \(k^\text{th}\) Walsh function,
• \(\beta(x) = - \lfloor \log_2(x) \rfloor\) and \(t_\nu(x) = 2^{-\nu \beta(x)}\) where \(\beta(0)=t_\nu(0) = 0\), and
• and \(\mu_\alpha\) is the Dick weight function which sums the first \(\alpha\) largest indices of \(1\) bits in the binary expansion of \(k\) e.g. \(k=13=1101_2\) has 1-bit indexes \((4,3,1)\) so

\[\mu_1(k) = 4, \mu_2(k) = 4+3, \mu_3(k) = 4+3+1 = \mu_4(k) = \mu_5(k) = \dots.\]

Examples:

>>> from qmcpy import DigitalNetB2, fwht
>>> n = 8
>>> d = 4
>>> dnb2 = DigitalNetB2(d,seed=11)
>>> x = dnb2(n,return_binary=True)
>>> x.shape
(8, 4)
>>> x.dtype
dtype('uint64')
>>> kernel = KernelDigShiftInvarAdaptiveAlpha(
...     d = d, 
...     t = dnb2.t,
...     alpha = list(range(1,d+1)),
...     scale = 10,
...     lengthscales = [1/j**2 for j in range(1,d+1)])
>>> k00 = kernel(x[0],x[0])
>>> k00.item()
48.084656084656096
>>> k0 = kernel(x,x[0])
>>> with np.printoptions(precision=2):
...     print(k0)
[48.08  0.    9.38  0.   19.74  0.   14.84  0.  ]
>>> assert k0[0]==k00
>>> kmat = kernel(x[:,None,:],x[None,:,:])
>>> with np.printoptions(precision=2):
...     print(kmat)
[[48.08  0.    9.38  0.   19.74  0.   14.84  0.  ]
 [ 0.   48.08  0.    9.38  0.   19.74  0.   14.84]
 [ 9.38  0.   48.08  0.   14.84  0.   19.74  0.  ]
 [ 0.    9.38  0.   48.08  0.   14.84  0.   19.74]
 [19.74  0.   14.84  0.   48.08  0.    9.38  0.  ]
 [ 0.   19.74  0.   14.84  0.   48.08  0.    9.38]
 [14.84  0.   19.74  0.    9.38  0.   48.08  0.  ]
 [ 0.   14.84  0.   19.74  0.    9.38  0.   48.08]]
>>> assert (kmat[:,0]==k0).all()
>>> lam = np.sqrt(n)*fwht(k0)
>>> y = np.random.Generator(np.random.PCG64(7)).uniform(low=0,high=1,size=(n))
>>> np.allclose(fwht(fwht(y)*lam),kmat@y)
True
>>> np.allclose(fwht(fwht(y)/lam),np.linalg.solve(kmat,y))
True
>>> import torch 
>>> xtorch = bin_from_numpy_to_torch(x)
>>> kernel_torch = KernelDigShiftInvarAdaptiveAlpha(
...     d = d, 
...     t = dnb2.t,
...     alpha = list(range(1,d+1)),
...     scale = 10,
...     lengthscales = [1/j**2 for j in range(1,d+1)],
...     torchify = True)
>>> kmat_torch = kernel_torch(xtorch[:,None,:],xtorch[None,:,:])
>>> np.allclose(kmat_torch.detach().numpy(),kmat,atol=1e-5)
True
>>> xf = to_float(x,dnb2.t)
>>> kmat_from_floats = kernel(xf[:,None,:],xf[None,:,:])
>>> np.allclose(kmat,kmat_from_floats)
True
>>> xftorch = to_float(xtorch,dnb2.t)
>>> xftorch.dtype
torch.float32
>>> kmat_torch_from_floats = kernel_torch(xftorch[:,None,:],xftorch[None,:,:])
>>> torch.allclose(kmat_torch_from_floats,kmat_torch)
True
>>> kernel.single_integral_01d(x)
array([10., 10., 10., 10., 10., 10., 10., 10.])
>>> kernel_torch.single_integral_01d(xtorch)
tensor([10., 10., 10., 10., 10., 10., 10., 10.], grad_fn=<AddBackward0>)

Batch Params

>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kernel = KernelDigShiftInvarAdaptiveAlpha(
...     d = 2, 
...     t = 10,
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,2])
>>> x = rng.uniform(low=0,high=1,size=(6,5,2))
>>> kernel(x,x).shape 
(4, 3, 6, 5)
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
(4, 3, 6, 5, 5)
>>> kfast = kernel(x[:,None,:,None,:],x[None,:,None,:,:])
>>> kfast.shape
(4, 3, 6, 6, 5, 5)
>>> kstable = kernel(x[:,None,:,None,:],x[None,:,None,:,:],stable=True)
>>> np.abs(kfast-kstable).max()
np.float64(4.440892098500626e-16)

References:

1. Dick, Josef, and Friedrich Pillichshammer.
   "Multivariate integration in weighted Hilbert spaces based on Walsh functions and weighted Sobolev spaces."
   Journal of Complexity 21.2 (2005): 149-195.

Parameters:

Name	Type	Description	Default
d	int	Dimension.	required
t	int	number of bits in binary represtnations. Typically dnb2.t where isinstance(dnb2,DigitalNetB2).	None
scale	Union[ndarray, Tensor]	Scaling factor \(S\).	1.0
lengthscales	Union[ndarray, Tensor]	Product weights \((\gamma_1,\dots,\gamma_d)\).	None
alpha	Union[ndarray, Tensor]	Smoothness parameters \((\alpha_1,\dots,\alpha_d)\) where \(\alpha_j \geq 1\) for \(j=1,\dots,d\).	1
shape_alpha	list	Shape of alpha when np.isscalar(alpha)	None
shape_scale	list	Shape of scale when np.isscalar(scale).	[1]
shape_lengthscales	list	Shape of lengthscales when np.isscalar(lengthscales)	None
tfs_scale	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_exp_eps_inv, tf_exp_eps)
tfs_lengthscales	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_exp_eps_inv, tf_exp_eps)
tfs_alpha	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_exp_eps_inv, tf_exp_eps)
torchify	bool	If True, use the torch backend. Set to True if computing gradients with respect to inputs and/or hyperparameters.	False
requires_grad_scale	bool	If True and torchify, set requires_grad=True for scale.	True
requires_grad_lengthscales	bool	If True and torchify, set requires_grad=True for lengthscales.	True
requires_grad_alpha	bool	If True and torchify, set requires_grad=True for alpha.	True
device	device	If torchify, put things onto this device.	'cpu'
compile_call	bool	If True, torch.compile the parsed___call__ method.	False
comiple_call_kwargs	dict	When compile_call is True, pass these keyword arguments to torch.compile.	{}

Source code in qmcpy/kernel/si_dsi_kernels.py

def __init__(
    self,
    d,
    t=None,
    scale=1.0,
    lengthscales=None,
    alpha=1,
    shape_scale=[1],
    shape_lengthscales=None,
    shape_alpha=None,
    tfs_scale=(tf_exp_eps_inv, tf_exp_eps),
    tfs_lengthscales=(tf_exp_eps_inv, tf_exp_eps),
    tfs_alpha=(tf_exp_eps_inv, tf_exp_eps),
    torchify=False,
    requires_grad_scale=True,
    requires_grad_lengthscales=True,
    requires_grad_alpha=True,
    device="cpu",
    compile_call=False,
    comiple_call_kwargs={},
):
    r"""
    Args:
        d (int): Dimension.
        t (int): number of bits in binary represtnations. Typically `dnb2.t` where `isinstance(dnb2,DigitalNetB2)`.
        scale (Union[np.ndarray, torch.Tensor]): Scaling factor $S$.
        lengthscales (Union[np.ndarray, torch.Tensor]): Product weights $(\gamma_1,\dots,\gamma_d)$.
        alpha (Union[np.ndarray, torch.Tensor]): Smoothness parameters $(\alpha_1,\dots,\alpha_d)$ where $\alpha_j \geq 1$ for $j=1,\dots,d$.
        shape_alpha (list): Shape of `alpha` when `np.isscalar(alpha)`
        shape_scale (list): Shape of `scale` when `np.isscalar(scale)`.
        shape_lengthscales (list): Shape of `lengthscales` when `np.isscalar(lengthscales)`
        tfs_scale (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_lengthscales (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_alpha (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        torchify (bool): If `True`, use the `torch` backend. Set to `True` if computing gradients with respect to inputs and/or hyperparameters.
        requires_grad_scale (bool): If `True` and `torchify`, set `requires_grad=True` for `scale`.
        requires_grad_lengthscales (bool): If `True` and `torchify`, set `requires_grad=True` for `lengthscales`.
        requires_grad_alpha (bool): If `True` and `torchify`, set `requires_grad=True` for `alpha`.
        device (torch.device): If `torchify`, put things onto this device.
        compile_call (bool): If `True`, `torch.compile` the `parsed___call__` method.
        comiple_call_kwargs (dict): When `compile_call` is `True`, pass these keyword arguments to `torch.compile`.
    """
    super().__init__(
        d=d,
        scale=scale,
        lengthscales=lengthscales,
        alpha=alpha,
        shape_alpha=[d] if shape_alpha is None else shape_alpha,
        alpha_endsize_ops=[d],
        shape_scale=shape_scale,
        shape_lengthscales=shape_lengthscales,
        tfs_alpha=tfs_alpha,
        tfs_scale=tfs_scale,
        tfs_lengthscales=tfs_lengthscales,
        torchify=torchify,
        requires_grad_alpha=requires_grad_alpha,
        requires_grad_scale=requires_grad_scale,
        requires_grad_lengthscales=requires_grad_lengthscales,
        device=device,
        compile_call=compile_call,
        comiple_call_kwargs=comiple_call_kwargs,
    )
    self.set_t(t)
    self.batch_param_names.append("alpha")

KernelDigShiftInvarCombined

Bases: AbstractSIDSIKernel

Digitally shift invariant kernel in base \(b=2\) with combination weights \(\boldsymbol{\alpha}_1,\dots,\boldsymbol{\alpha}_d \in \mathbb{R}_{>0}^4\), smoothness \(\boldsymbol{\alpha}\), product weights \(\boldsymbol{\gamma}\), and scale \(S\):

\[\begin{aligned} K(\boldsymbol{x},\boldsymbol{z}) &= S \prod_{j=1}^d \left(1+ \gamma_j \left(\sum_{p=1}^4 \alpha_{jp} \tilde{K}_p(x_j \oplus z_j)\right)\right) \end{aligned}\]

where, \(\oplus\) is defined in the docs for KernelDigShiftInvar and so are \(\tilde{K}_p\) for \(p \in \{1,2,3,4\}\)

Examples:

>>> from qmcpy import DigitalNetB2, fwht
>>> n = 8
>>> d = 4
>>> dnb2 = DigitalNetB2(d,seed=11)
>>> x = dnb2(n,return_binary=True)
>>> x.shape
(8, 4)
>>> x.dtype
dtype('uint64')
>>> kernel = KernelDigShiftInvarCombined(
...     d = d,
...     t = dnb2.t,
...     scale = 10,
...     lengthscales = [1/j**2 for j in range(1,d+1)])
>>> k00 = kernel(x[0],x[0])
>>> k00.item()
306.66030731930863
>>> k0 = kernel(x,x[0])
>>> with np.printoptions(precision=2):
...     print(k0)
[306.66  -1.65   6.88 -13.5   22.57  -9.34   7.99  -7.54]
>>> assert k0[0]==k00
>>> kmat = kernel(x[:,None,:],x[None,:,:])
>>> with np.printoptions(precision=2):
...     print(kmat)
[[306.66  -1.65   6.88 -13.5   22.57  -9.34   7.99  -7.54]
 [ -1.65 306.66 -13.5    6.88  -9.34  22.57  -7.54   7.99]
 [  6.88 -13.5  306.66  -1.65   7.99  -7.54  22.57  -9.34]
 [-13.5    6.88  -1.65 306.66  -7.54   7.99  -9.34  22.57]
 [ 22.57  -9.34   7.99  -7.54 306.66  -1.65   6.88 -13.5 ]
 [ -9.34  22.57  -7.54   7.99  -1.65 306.66 -13.5    6.88]
 [  7.99  -7.54  22.57  -9.34   6.88 -13.5  306.66  -1.65]
 [ -7.54   7.99  -9.34  22.57 -13.5    6.88  -1.65 306.66]]
>>> assert (kmat[:,0]==k0).all()
>>> lam = np.sqrt(n)*fwht(k0)
>>> y = np.random.Generator(np.random.PCG64(7)).uniform(low=0,high=1,size=(n))
>>> np.allclose(fwht(fwht(y)*lam),kmat@y)
True
>>> np.allclose(fwht(fwht(y)/lam),np.linalg.solve(kmat,y))
True
>>> import torch
>>> xtorch = bin_from_numpy_to_torch(x)
>>> kernel_torch = KernelDigShiftInvarCombined(
...     d = d,
...     t = dnb2.t,
...     scale = 10,
...     lengthscales = [1/j**2 for j in range(1,d+1)],
...     torchify = True)
>>> kmat_torch = kernel_torch(xtorch[:,None,:],xtorch[None,:,:])
>>> np.allclose(kmat_torch.detach().numpy(),kmat)
True
>>> xf = to_float(x,dnb2.t)
>>> kmat_from_floats = kernel(xf[:,None,:],xf[None,:,:])
>>> np.allclose(kmat,kmat_from_floats)
True
>>> xftorch = to_float(xtorch,dnb2.t)
>>> xftorch.dtype
torch.float32
>>> kmat_torch_from_floats = kernel_torch(xftorch[:,None,:],xftorch[None,:,:])
>>> torch.allclose(kmat_torch_from_floats,kmat_torch)
True
>>> kernel.single_integral_01d(x)
array([10., 10., 10., 10., 10., 10., 10., 10.])
>>> kernel_torch.single_integral_01d(xtorch)
tensor([10., 10., 10., 10., 10., 10., 10., 10.], grad_fn=<AddBackward0>)

Batch Params

>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kernel = KernelDigShiftInvarCombined(
...     d = 2,
...     t = 10,
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,2])
>>> x = rng.uniform(low=0,high=1,size=(6,5,2))
>>> kernel(x,x).shape
(4, 3, 6, 5)
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
(4, 3, 6, 5, 5)
>>> kfast = kernel(x[:,None,:,None,:],x[None,:,None,:,:])
>>> kfast.shape
(4, 3, 6, 6, 5, 5)
>>> kstable = kernel(x[:,None,:,None,:],x[None,:,None,:,:],stable=True)
>>> np.abs(kfast-kstable).max()
np.float64(8.881784197001252e-16)

Parameters:

Name	Type	Description	Default
d	int	Dimension.	required
t	int	number of bits in binary represtnations. Typically dnb2.t where isinstance(dnb2,DigitalNetB2).	None
scale	Union[ndarray, Tensor]	Scaling factor \(S\).	1.0
lengthscales	Union[ndarray, Tensor]	Product weights \((\gamma_1,\dots,\gamma_d)\).	None
alpha	Union[ndarray, Tensor]	Weights \(\boldsymbol{\alpha}_1,\dots,\boldsymbol{\alpha}_d \in \mathbb{R}_{>0}^4\).	1.0
shape_scale	list	Shape of scale when np.isscalar(scale).	[1]
shape_lengthscales	list	Shape of lengthscales when np.isscalar(lengthscales)	None
shape_alpha	list	Shape of alpha when np.isscalar(alpha)	None
tfs_scale	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_exp_eps_inv, tf_exp_eps)
tfs_lengthscales	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_exp_eps_inv, tf_exp_eps)
torchify	bool	If True, use the torch backend. Set to True if computing gradients with respect to inputs and/or hyperparameters.	False
requires_grad_scale	bool	If True and torchify, set requires_grad=True for scale.	True
requires_grad_lengthscales	bool	If True and torchify, set requires_grad=True for lengthscales.	True
requires_grad_alpha	bool	If True and torchify, set requires_grad=True for alpha.	True
device	device	If torchify, put things onto this device.	'cpu'
compile_call	bool	If True, torch.compile the parsed___call__ method.	False
comiple_call_kwargs	dict	When compile_call is True, pass these keyword arguments to torch.compile.	{}

Source code in qmcpy/kernel/si_dsi_kernels.py

def __init__(
    self,
    d,
    t=None,
    scale=1.0,
    lengthscales=None,
    alpha=1.0,
    shape_scale=[1],
    shape_lengthscales=None,
    shape_alpha=None,
    tfs_scale=(tf_exp_eps_inv, tf_exp_eps),
    tfs_lengthscales=(tf_exp_eps_inv, tf_exp_eps),
    tfs_alpha=(tf_exp_eps_inv, tf_exp_eps),
    torchify=False,
    requires_grad_scale=True,
    requires_grad_lengthscales=True,
    requires_grad_alpha=True,
    device="cpu",
    compile_call=False,
    comiple_call_kwargs={},
):
    r"""
    Args:
        d (int): Dimension.
        t (int): number of bits in binary represtnations. Typically `dnb2.t` where `isinstance(dnb2,DigitalNetB2)`.
        scale (Union[np.ndarray, torch.Tensor]): Scaling factor $S$.
        lengthscales (Union[np.ndarray, torch.Tensor]): Product weights $(\gamma_1,\dots,\gamma_d)$.
        alpha (Union[np.ndarray, torch.Tensor]): Weights $\boldsymbol{\alpha}_1,\dots,\boldsymbol{\alpha}_d \in \mathbb{R}_{>0}^4$.
        shape_scale (list): Shape of `scale` when `np.isscalar(scale)`.
        shape_lengthscales (list): Shape of `lengthscales` when `np.isscalar(lengthscales)`
        shape_alpha (list): Shape of `alpha` when `np.isscalar(alpha)`
        tfs_scale (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_lengthscales (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        torchify (bool): If `True`, use the `torch` backend. Set to `True` if computing gradients with respect to inputs and/or hyperparameters.
        requires_grad_scale (bool): If `True` and `torchify`, set `requires_grad=True` for `scale`.
        requires_grad_lengthscales (bool): If `True` and `torchify`, set `requires_grad=True` for `lengthscales`.
        requires_grad_alpha (bool): If `True` and `torchify`, set `requires_grad=True` for `alpha`.
        device (torch.device): If `torchify`, put things onto this device.
        compile_call (bool): If `True`, `torch.compile` the `parsed___call__` method.
        comiple_call_kwargs (dict): When `compile_call` is `True`, pass these keyword arguments to `torch.compile`.
    """
    super().__init__(
        d=d,
        scale=scale,
        lengthscales=lengthscales,
        alpha=alpha,
        shape_alpha=[4, d] if shape_alpha is None else shape_alpha,
        alpha_endsize_ops=[d],
        shape_scale=shape_scale,
        shape_lengthscales=shape_lengthscales,
        tfs_alpha=tfs_alpha,
        tfs_scale=tfs_scale,
        tfs_lengthscales=tfs_lengthscales,
        torchify=torchify,
        requires_grad_alpha=requires_grad_alpha,
        requires_grad_scale=requires_grad_scale,
        requires_grad_lengthscales=requires_grad_lengthscales,
        device=device,
        compile_call=compile_call,
        comiple_call_kwargs=comiple_call_kwargs,
    )
    self.set_t(t)
    assert self.alpha.shape[-2:] == (4, d)

KernelGaussian

Bases: AbstractKernelGaussianSE

Gaussian / Squared Exponential kernel implemented using the product of exponentials.

\[K(\boldsymbol{x},\boldsymbol{z}) = S \prod_{j=1}^d \exp\left(-\left(\frac{x_j-z_j}{\sqrt{2} \gamma_j}\right)^2\right)\]

Examples:

>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kernel = KernelGaussian(d=2)
>>> x = rng.uniform(low=0,high=1,size=(4,2))
>>> kernel(x[0],x[0]).item()
1.0
>>> kernel(x,x)
array([1., 1., 1., 1.])
>>> kernel(x[:,None,:],x[None,:,:])
array([[1.        , 0.78888466, 0.9483142 , 0.8228498 ],
       [0.78888466, 1.        , 0.72380317, 0.62226176],
       [0.9483142 , 0.72380317, 1.        , 0.95613874],
       [0.8228498 , 0.62226176, 0.95613874, 1.        ]])

Multiple randomizations

>>> x = rng.uniform(low=0,high=1,size=(6,5,2))
>>> kernel(x,x).shape
(6, 5)
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
(6, 5, 5)
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
(6, 6, 5, 5)

Batch hyperparameters

>>> kernel = KernelGaussian(
...     d = 2,
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,1])
>>> kernel(x,x).shape
(4, 3, 6, 5)
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
(4, 3, 6, 5, 5)
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
(4, 3, 6, 6, 5, 5)

Integrals

>>> kernel = KernelGaussian(
...     d = 2,
...     scale = rng.uniform(low=0,high=1,size=(3,1,)),
...     lengthscales = rng.uniform(low=0,high=1,size=(1,2)))
>>> kintint = kernel.double_integral_01d()
>>> kintint
array([0.50079567, 0.11125229, 0.34760005])
>>> x_qmc_4d = DigitalNetB2(4,seed=7)(2**16)
>>> kintint_qmc = kernel(x_qmc_4d[:,:2],x_qmc_4d[:,2:]).mean(1)
>>> kintint_qmc
array([0.5007959 , 0.11125234, 0.34760021])
>>> with np.printoptions(formatter={"float":lambda x: "%.1e"%x}):
...     np.abs(kintint-kintint_qmc)
array([2.3e-07, 5.1e-08, 1.6e-07])
>>> x = rng.uniform(low=0,high=1,size=(4,2))
>>> kint = kernel.single_integral_01d(x)
>>> kint
array([[0.54610372, 0.4272801 , 0.43001936, 0.44688778],
       [0.12131753, 0.09492073, 0.09552926, 0.0992766 ],
       [0.37904817, 0.29657323, 0.29847453, 0.31018283]])
>>> x_qmc_2d = DigitalNetB2(2,seed=7)(2**16)
>>> kint_qmc = kernel(x[:,None,:],x_qmc_2d).mean(-1)
>>> kint_qmc
array([[0.54610372, 0.4272801 , 0.43001936, 0.44688778],
       [0.12131753, 0.09492073, 0.09552926, 0.0992766 ],
       [0.37904817, 0.29657323, 0.29847453, 0.31018283]])
>>> with np.printoptions(formatter={"float":lambda x: "%.1e"%x}):
...     np.abs(kint-kint_qmc)
array([[2.2e-13, 4.0e-13, 4.2e-12, 3.9e-12],
       [4.8e-14, 8.9e-14, 9.4e-13, 8.8e-13],
       [1.5e-13, 2.8e-13, 2.9e-12, 2.7e-12]])
>>> k_1l = KernelGaussian(d=2,lengthscales=[.5])
>>> k_2l = KernelGaussian(d=2,lengthscales=[.5,.5])
>>> print("%.5f"%k_2l.double_integral_01d())
0.58363
>>> print("%.5f"%k_1l.double_integral_01d())
0.58363
>>> k_1l.single_integral_01d(x)
array([0.58119655, 0.46559577, 0.57603494, 0.53494393])
>>> k_2l.single_integral_01d(x)
array([0.58119655, 0.46559577, 0.57603494, 0.53494393])

PyTorch

>>> import torch
>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kernel = KernelGaussian(d=2,torchify=True)
>>> x = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,2)))
>>> kernel(x[0],x[0]).item()
1.0
>>> kernel(x,x)
tensor([1., 1., 1., 1.], dtype=torch.float64, grad_fn=<MulBackward0>)
>>> kernel(x[:,None,:],x[None,:,:])
tensor([[1.0000, 0.7889, 0.9483, 0.8228],
        [0.7889, 1.0000, 0.7238, 0.6223],
        [0.9483, 0.7238, 1.0000, 0.9561],
        [0.8228, 0.6223, 0.9561, 1.0000]], dtype=torch.float64,
       grad_fn=<MulBackward0>)

Multiple randomizations

>>> x = torch.from_numpy(rng.uniform(low=0,high=1,size=(6,5,2)))
>>> kernel(x,x).shape
torch.Size([6, 5])
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
torch.Size([6, 5, 5])
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
torch.Size([6, 6, 5, 5])

Batch hyperparameters

>>> kernel = KernelGaussian(
...     d = 2,
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,1],
...     torchify = True)
>>> kernel(x,x).shape
torch.Size([4, 3, 6, 5])
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
torch.Size([4, 3, 6, 5, 5])
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
torch.Size([4, 3, 6, 6, 5, 5])

Integrals

>>> kernel = KernelGaussian(
...     d = 2,
...     torchify = True,
...     scale = torch.from_numpy(rng.uniform(low=0,high=1,size=(3,1,))),
...     lengthscales = torch.from_numpy(rng.uniform(low=0,high=1,size=(1,2))),
...     requires_grad_scale = False,
...     requires_grad_lengthscales = False)
>>> kintint = kernel.double_integral_01d()
>>> kintint
tensor([0.5008, 0.1113, 0.3476], dtype=torch.float64)
>>> x_qmc_4d = torch.from_numpy(DigitalNetB2(4,seed=7)(2**16))
>>> kintint_qmc = kernel(x_qmc_4d[:,:2],x_qmc_4d[:,2:]).mean(1)
>>> kintint_qmc
tensor([0.5008, 0.1113, 0.3476], dtype=torch.float64)
>>> with np.printoptions(formatter={"float":lambda x: "%.1e"%x}):
...     torch.abs(kintint-kintint_qmc).numpy()
array([2.3e-07, 5.1e-08, 1.6e-07])
>>> x = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,2)))
>>> kint = kernel.single_integral_01d(x)
>>> kint
tensor([[0.5461, 0.4273, 0.4300, 0.4469],
        [0.1213, 0.0949, 0.0955, 0.0993],
        [0.3790, 0.2966, 0.2985, 0.3102]], dtype=torch.float64)
>>> x_qmc_2d = torch.from_numpy(DigitalNetB2(2,seed=7)(2**16))
>>> kint_qmc = kernel(x[:,None,:],x_qmc_2d).mean(-1)
>>> kint_qmc
tensor([[0.5461, 0.4273, 0.4300, 0.4469],
        [0.1213, 0.0949, 0.0955, 0.0993],
        [0.3790, 0.2966, 0.2985, 0.3102]], dtype=torch.float64)
>>> with np.printoptions(formatter={"float":lambda x: "%.1e"%x}):
...     torch.abs(kint-kint_qmc).numpy()
array([[2.2e-13, 4.0e-13, 4.2e-12, 3.9e-12],
       [4.8e-14, 8.9e-14, 9.4e-13, 8.8e-13],
       [1.5e-13, 2.8e-13, 2.9e-12, 2.7e-12]])
>>> k_1l = KernelGaussian(d=2,lengthscales=[.5],torchify=True)
>>> k_2l = KernelGaussian(d=2,lengthscales=[.5,.5],torchify=True)
>>> k_2l.double_integral_01d()
tensor(0.5836, grad_fn=<MulBackward0>)
>>> k_1l.double_integral_01d()
tensor(0.5836, grad_fn=<MulBackward0>)
>>> k_1l.single_integral_01d(x)
tensor([0.5812, 0.4656, 0.5760, 0.5349], dtype=torch.float64,
       grad_fn=<MulBackward0>)
>>> k_2l.single_integral_01d(x)
tensor([0.5812, 0.4656, 0.5760, 0.5349], dtype=torch.float64,
       grad_fn=<MulBackward0>)

Derivatives

>>> kernel = KernelGaussian(
...     d = 3,
...     torchify = True,
...     scale = torch.from_numpy(rng.uniform(low=0,high=1,size=(1,))),
...     lengthscales = torch.from_numpy(rng.uniform(low=0,high=1,size=(3,))))
>>> x0 = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,))).requires_grad_(True)
>>> x1 = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,))).requires_grad_(True)
>>> x2 = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,))).requires_grad_(True)
>>> x = torch.stack([x0,x1,x2],axis=-1)
>>> z0 = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,))).requires_grad_(True)
>>> z1 = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,))).requires_grad_(True)
>>> z2 = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,))).requires_grad_(True)
>>> z = torch.stack([z0,z1,z2],axis=-1)
>>> c = torch.from_numpy(rng.uniform(low=0,high=1,size=(2,)))
>>> beta0 = torch.tensor([
...     [1,0,0],
...     [0,2,0]])
>>> beta1 = torch.tensor([
...     [0,0,2],
...     [2,1,0]])
>>> with torch.no_grad():
...     y = kernel(x,z,beta0,beta1,c)
>>> y
tensor([ 97.6657,   3.8621, -65.9329,  -1.1932], dtype=torch.float64)
>>> y_no_deriv = kernel(x,z)
>>> y_first = y_no_deriv.clone()
>>> y_first = torch.autograd.grad(y_first,x0,grad_outputs=torch.ones_like(y_first,requires_grad=True),create_graph=True)[0]
>>> y_first = torch.autograd.grad(y_first,z2,grad_outputs=torch.ones_like(y_first,requires_grad=True),create_graph=True)[0]
>>> y_first = torch.autograd.grad(y_first,z2,grad_outputs=torch.ones_like(y_first,requires_grad=True),create_graph=True)[0]
>>> y_second = y_no_deriv.clone()
>>> y_second = torch.autograd.grad(y_second,x1,grad_outputs=torch.ones_like(y_second,requires_grad=True),create_graph=True)[0]
>>> y_second = torch.autograd.grad(y_second,x1,grad_outputs=torch.ones_like(y_second,requires_grad=True),create_graph=True)[0]
>>> y_second = torch.autograd.grad(y_second,z0,grad_outputs=torch.ones_like(y_second,requires_grad=True),create_graph=True)[0]
>>> y_second = torch.autograd.grad(y_second,z0,grad_outputs=torch.ones_like(y_second,requires_grad=True),create_graph=True)[0]
>>> y_second = torch.autograd.grad(y_second,z1,grad_outputs=torch.ones_like(y_second,requires_grad=True),create_graph=True)[0]
>>> yhat = (y_first*c[0]+y_second*c[1]).detach()
>>> yhat
tensor([ 97.6657,   3.8621, -65.9329,  -1.1932], dtype=torch.float64)
>>> torch.allclose(y,yhat)
True

Source code in qmcpy/kernel/abstract_kernel.py

def __init__(
    self,
    d,
    scale=1.0,
    lengthscales=1.0,
    shape_scale=[1],
    shape_lengthscales=None,
    tfs_scale=(tf_exp_eps_inv, tf_exp_eps),
    tfs_lengthscales=(tf_exp_eps_inv, tf_exp_eps),
    torchify=False,
    requires_grad_scale=True,
    requires_grad_lengthscales=True,
    device="cpu",
    compile_call=False,
    comiple_call_kwargs={},
):
    r"""
    Args:
        d (int): Dimension.
        scale (Union[np.ndarray, torch.Tensor]): Scaling factor $S$.
        lengthscales (Union[np.ndarray, torch.Tensor]): Lengthscales $\boldsymbol{\gamma}$.
        shape_scale (list): Shape of `scale` when `np.isscalar(scale)`.
        shape_lengthscales (list): Shape of `lengthscales` when `np.isscalar(lengthscales)`
        tfs_scale (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_lengthscales (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        torchify (bool): If `True`, use the `torch` backend. Set to `True` if computing gradients with respect to inputs and/or hyperparameters.
        requires_grad_scale (bool): If `True` and `torchify`, set `requires_grad=True` for `scale`.
        requires_grad_lengthscales (bool): If `True` and `torchify`, set `requires_grad=True` for `lengthscales`.
        device (torch.device): If `torchify`, put things onto this device.
        compile_call (bool): If `True`, `torch.compile` the `parsed___call__` method.
        comiple_call_kwargs (dict): When `compile_call` is `True`, pass these keyword arguments to `torch.compile`.
    """
    super().__init__(
        d=d,
        torchify=torchify,
        device=device,
        compile_call=compile_call,
        comiple_call_kwargs=comiple_call_kwargs,
    )
    self.raw_scale = self.parse_assign_param(
        pname="scale",
        param=scale,
        shape_param=shape_scale,
        requires_grad_param=requires_grad_scale,
        tfs_param=tfs_scale,
        endsize_ops=[1],
        constraints=["POSITIVE"],
    )
    self.tfs_scale = tfs_scale
    self.batch_param_names.append("scale")
    self.raw_lengthscales = self.parse_assign_param(
        pname="lengthscales",
        param=lengthscales,
        shape_param=[self.d] if shape_lengthscales is None else shape_lengthscales,
        requires_grad_param=requires_grad_lengthscales,
        tfs_param=tfs_lengthscales,
        endsize_ops=[1, self.d],
        constraints=["POSITIVE"],
    )
    self.tfs_lengthscales = tfs_lengthscales
    self.batch_param_names.append("lengthscales")

KernelSquaredExponential

Bases: AbstractKernelGaussianSE

Gaussian / Squared Exponential kernel implemented using the pairwise distance function. Please use KernelGaussian when using derivative information.

\[K(\boldsymbol{x},\boldsymbol{z}) = S \exp\left(-d_{\boldsymbol{\gamma}}^2(\boldsymbol{x},\boldsymbol{z})\right), \qquad d_{\boldsymbol{\gamma}}(\boldsymbol{x},\boldsymbol{z}) = \left\lVert\frac{\boldsymbol{x}-\boldsymbol{z}}{\sqrt{2}\boldsymbol{\gamma}}\right\rVert_2.\]

Examples:

>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kernel = KernelSquaredExponential(d=2)
>>> x = rng.uniform(low=0,high=1,size=(4,2))
>>> kernel(x[0],x[0]).item()
1.0
>>> kernel(x,x)
array([1., 1., 1., 1.])
>>> kernel(x[:,None,:],x[None,:,:])
array([[1.        , 0.78888466, 0.9483142 , 0.8228498 ],
       [0.78888466, 1.        , 0.72380317, 0.62226176],
       [0.9483142 , 0.72380317, 1.        , 0.95613874],
       [0.8228498 , 0.62226176, 0.95613874, 1.        ]])

Multiple randomizations

>>> x = rng.uniform(low=0,high=1,size=(6,5,2))
>>> kernel(x,x).shape
(6, 5)
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
(6, 5, 5)
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
(6, 6, 5, 5)

Batch hyperparameters

>>> kernel = KernelSquaredExponential(
...     d = 2,
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,1])
>>> kernel(x,x).shape
(4, 3, 6, 5)
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
(4, 3, 6, 5, 5)
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
(4, 3, 6, 6, 5, 5)

PyTorch

>>> import torch
>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kernel = KernelSquaredExponential(d=2,torchify=True)
>>> x = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,2)))
>>> kernel(x[0],x[0]).item()
1.0
>>> kernel(x,x)
tensor([1., 1., 1., 1.], dtype=torch.float64, grad_fn=<MulBackward0>)
>>> kernel(x[:,None,:],x[None,:,:])
tensor([[1.0000, 0.7889, 0.9483, 0.8228],
        [0.7889, 1.0000, 0.7238, 0.6223],
        [0.9483, 0.7238, 1.0000, 0.9561],
        [0.8228, 0.6223, 0.9561, 1.0000]], dtype=torch.float64,
       grad_fn=<MulBackward0>)

Multiple randomizations

>>> x = torch.from_numpy(rng.uniform(low=0,high=1,size=(6,5,2)))
>>> kernel(x,x).shape
torch.Size([6, 5])
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
torch.Size([6, 5, 5])
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
torch.Size([6, 6, 5, 5])

Batch hyperparameters

>>> kernel = KernelSquaredExponential(
...     d = 2,
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,1],
...     torchify = True)
>>> kernel(x,x).shape
torch.Size([4, 3, 6, 5])
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
torch.Size([4, 3, 6, 5, 5])
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
torch.Size([4, 3, 6, 6, 5, 5])

Source code in qmcpy/kernel/abstract_kernel.py

def __init__(
    self,
    d,
    scale=1.0,
    lengthscales=1.0,
    shape_scale=[1],
    shape_lengthscales=None,
    tfs_scale=(tf_exp_eps_inv, tf_exp_eps),
    tfs_lengthscales=(tf_exp_eps_inv, tf_exp_eps),
    torchify=False,
    requires_grad_scale=True,
    requires_grad_lengthscales=True,
    device="cpu",
    compile_call=False,
    comiple_call_kwargs={},
):
    r"""
    Args:
        d (int): Dimension.
        scale (Union[np.ndarray, torch.Tensor]): Scaling factor $S$.
        lengthscales (Union[np.ndarray, torch.Tensor]): Lengthscales $\boldsymbol{\gamma}$.
        shape_scale (list): Shape of `scale` when `np.isscalar(scale)`.
        shape_lengthscales (list): Shape of `lengthscales` when `np.isscalar(lengthscales)`
        tfs_scale (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_lengthscales (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        torchify (bool): If `True`, use the `torch` backend. Set to `True` if computing gradients with respect to inputs and/or hyperparameters.
        requires_grad_scale (bool): If `True` and `torchify`, set `requires_grad=True` for `scale`.
        requires_grad_lengthscales (bool): If `True` and `torchify`, set `requires_grad=True` for `lengthscales`.
        device (torch.device): If `torchify`, put things onto this device.
        compile_call (bool): If `True`, `torch.compile` the `parsed___call__` method.
        comiple_call_kwargs (dict): When `compile_call` is `True`, pass these keyword arguments to `torch.compile`.
    """
    super().__init__(
        d=d,
        torchify=torchify,
        device=device,
        compile_call=compile_call,
        comiple_call_kwargs=comiple_call_kwargs,
    )
    self.raw_scale = self.parse_assign_param(
        pname="scale",
        param=scale,
        shape_param=shape_scale,
        requires_grad_param=requires_grad_scale,
        tfs_param=tfs_scale,
        endsize_ops=[1],
        constraints=["POSITIVE"],
    )
    self.tfs_scale = tfs_scale
    self.batch_param_names.append("scale")
    self.raw_lengthscales = self.parse_assign_param(
        pname="lengthscales",
        param=lengthscales,
        shape_param=[self.d] if shape_lengthscales is None else shape_lengthscales,
        requires_grad_param=requires_grad_lengthscales,
        tfs_param=tfs_lengthscales,
        endsize_ops=[1, self.d],
        constraints=["POSITIVE"],
    )
    self.tfs_lengthscales = tfs_lengthscales
    self.batch_param_names.append("lengthscales")

KernelRationalQuadratic

Bases: AbstractKernelScaleLengthscales

Rational Quadratic kernel

\[K(\boldsymbol{x},\boldsymbol{z}) = S \left(1+\frac{d_{\boldsymbol{\gamma}}^2(\boldsymbol{x},\boldsymbol{z})}{\alpha}\right)^{-\alpha}, \qquad d_{\boldsymbol{\gamma}}(\boldsymbol{x},\boldsymbol{z}) = \left\lVert\frac{\boldsymbol{x}-\boldsymbol{z}}{\sqrt{2}\boldsymbol{\gamma}}\right\rVert_2.\]

Examples:

>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kernel = KernelRationalQuadratic(d=2)
>>> x = rng.uniform(low=0,high=1,size=(4,2))
>>> kernel(x[0],x[0]).item()
1.0
>>> kernel(x,x)
array([1., 1., 1., 1.])
>>> kernel(x[:,None,:],x[None,:,:])
array([[1.        , 0.80831912, 0.94960504, 0.83683297],
       [0.80831912, 1.        , 0.75572321, 0.67824456],
       [0.94960504, 0.75572321, 1.        , 0.95707312],
       [0.83683297, 0.67824456, 0.95707312, 1.        ]])

Multiple randomizations

>>> x = rng.uniform(low=0,high=1,size=(6,5,2))
>>> kernel(x,x).shape
(6, 5)
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
(6, 5, 5)
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
(6, 6, 5, 5)

Batch hyperparameters

>>> kernel = KernelRationalQuadratic(
...     d = 2,
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,1])
>>> kernel(x,x).shape
(4, 3, 6, 5)
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
(4, 3, 6, 5, 5)
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
(4, 3, 6, 6, 5, 5)

PyTorch

>>> import torch
>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kernel = KernelRationalQuadratic(d=2,torchify=True)
>>> x = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,2)))
>>> kernel(x[0],x[0]).item()
1.0
>>> kernel(x,x)
tensor([1., 1., 1., 1.], dtype=torch.float64, grad_fn=<MulBackward0>)
>>> kernel(x[:,None,:],x[None,:,:])
tensor([[1.0000, 0.8083, 0.9496, 0.8368],
        [0.8083, 1.0000, 0.7557, 0.6782],
        [0.9496, 0.7557, 1.0000, 0.9571],
        [0.8368, 0.6782, 0.9571, 1.0000]], dtype=torch.float64,
       grad_fn=<MulBackward0>)

Multiple randomizations

>>> x = torch.from_numpy(rng.uniform(low=0,high=1,size=(6,5,2)))
>>> kernel(x,x).shape
torch.Size([6, 5])
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
torch.Size([6, 5, 5])
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
torch.Size([6, 6, 5, 5])

Batch hyperparameters

>>> kernel = KernelRationalQuadratic(
...     d = 2,
...     shape_alpha = [4,3,1],
...     shape_lengthscales = [3,1],
...     torchify = True)
>>> kernel(x,x).shape
torch.Size([4, 3, 6, 5])
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
torch.Size([4, 3, 6, 5, 5])
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
torch.Size([4, 3, 6, 6, 5, 5])

Parameters:

Name	Type	Description	Default
d	int	Dimension.	required
scale	Union[ndarray, Tensor]	Scaling factor \(S\).	1.0
lengthscales	Union[ndarray, Tensor]	Lengthscales \(\boldsymbol{\gamma}\).	1.0
alpha	Union[ndarray, Tensor]	Scale mixture parameter \(\alpha\).	1.0
shape_scale	list	Shape of scale when np.isscalar(scale).	[1]
shape_lengthscales	list	Shape of lengthscales when np.isscalar(lengthscales)	None
shape_alpha	list	Shape of alpha when np.isscalar(alpha)	[1]
tfs_scale	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_exp_eps_inv, tf_exp_eps)
tfs_lengthscales	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_exp_eps_inv, tf_exp_eps)
tfs_alpha	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_exp_eps_inv, tf_exp_eps)
torchify	bool	If True, use the torch backend. Set to True if computing gradients with respect to inputs and/or hyperparameters.	False
requires_grad_scale	bool	If True and torchify, set requires_grad=True for scale.	True
requires_grad_lengthscales	bool	If True and torchify, set requires_grad=True for lengthscales.	True
requires_grad_alpha	bool	If True and torchify, set requires_grad=True for alpha.	True
device	device	If torchify, put things onto this device.	'cpu'
compile_call	bool	If True, torch.compile the parsed___call__ method.	False
comiple_call_kwargs	dict	When compile_call is True, pass these keyword arguments to torch.compile.	{}

Source code in qmcpy/kernel/common_kernels.py

def __init__(
    self,
    d,
    scale=1.0,
    lengthscales=1.0,
    alpha=1.0,
    shape_scale=[1],
    shape_lengthscales=None,
    shape_alpha=[1],
    tfs_scale=(tf_exp_eps_inv, tf_exp_eps),
    tfs_lengthscales=(tf_exp_eps_inv, tf_exp_eps),
    tfs_alpha=(tf_exp_eps_inv, tf_exp_eps),
    torchify=False,
    requires_grad_scale=True,
    requires_grad_lengthscales=True,
    requires_grad_alpha=True,
    device="cpu",
    compile_call=False,
    comiple_call_kwargs={},
):
    r"""
    Args:
        d (int): Dimension.
        scale (Union[np.ndarray, torch.Tensor]): Scaling factor $S$.
        lengthscales (Union[np.ndarray, torch.Tensor]): Lengthscales $\boldsymbol{\gamma}$.
        alpha (Union[np.ndarray, torch.Tensor]): Scale mixture parameter $\alpha$.
        shape_scale (list): Shape of `scale` when `np.isscalar(scale)`.
        shape_lengthscales (list): Shape of `lengthscales` when `np.isscalar(lengthscales)`
        shape_alpha (list): Shape of `alpha` when `np.isscalar(alpha)`
        tfs_scale (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_lengthscales (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_alpha (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        torchify (bool): If `True`, use the `torch` backend. Set to `True` if computing gradients with respect to inputs and/or hyperparameters.
        requires_grad_scale (bool): If `True` and `torchify`, set `requires_grad=True` for `scale`.
        requires_grad_lengthscales (bool): If `True` and `torchify`, set `requires_grad=True` for `lengthscales`.
        requires_grad_alpha (bool): If `True` and `torchify`, set `requires_grad=True` for `alpha`.
        device (torch.device): If `torchify`, put things onto this device.
        compile_call (bool): If `True`, `torch.compile` the `parsed___call__` method.
        comiple_call_kwargs (dict): When `compile_call` is `True`, pass these keyword arguments to `torch.compile`.
    """
    super().__init__(
        d=d,
        scale=scale,
        lengthscales=lengthscales,
        shape_scale=shape_scale,
        shape_lengthscales=shape_lengthscales,
        tfs_scale=tfs_scale,
        tfs_lengthscales=tfs_lengthscales,
        torchify=torchify,
        requires_grad_scale=requires_grad_scale,
        requires_grad_lengthscales=requires_grad_lengthscales,
        device=device,
        compile_call=compile_call,
        comiple_call_kwargs=comiple_call_kwargs,
    )
    self.raw_alpha = self.parse_assign_param(
        pname="alpha",
        param=alpha,
        shape_param=shape_alpha,
        requires_grad_param=requires_grad_alpha,
        tfs_param=tfs_alpha,
        endsize_ops=[1],
        constraints=["POSITIVE"],
    )
    self.tfs_alpha = tfs_alpha
    self.batch_param_names.append("alpha")

KernelMatern12

Bases: AbstractKernelScaleLengthscales

Matern kernel with \(\alpha=1/2\).

\[K(\boldsymbol{x},\boldsymbol{z}) = S \exp\left(-d_{\boldsymbol{\gamma}}(\boldsymbol{x},\boldsymbol{z})\right), \qquad d_{\boldsymbol{\gamma}}(\boldsymbol{x},\boldsymbol{z}) = \left\lVert\frac{\boldsymbol{x}-\boldsymbol{z}}{\sqrt{2}\boldsymbol{\gamma}}\right\rVert_2.\]

Examples:

>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kernel = KernelMatern12(d=2)
>>> x = rng.uniform(low=0,high=1,size=(4,2))
>>> kernel(x[0],x[0]).item()
1.0
>>> kernel(x,x)
array([1., 1., 1., 1.])
>>> kernel(x[:,None,:],x[None,:,:])
array([[1.        , 0.61448839, 0.79424131, 0.64302787],
       [0.61448839, 1.        , 0.56635268, 0.50219691],
       [0.79424131, 0.56635268, 1.        , 0.80913986],
       [0.64302787, 0.50219691, 0.80913986, 1.        ]])

Multiple randomizations

>>> x = rng.uniform(low=0,high=1,size=(6,5,2))
>>> kernel(x,x).shape
(6, 5)
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
(6, 5, 5)
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
(6, 6, 5, 5)

Batch hyperparameters

>>> kernel = KernelMatern12(
...     d = 2,
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,1])
>>> kernel(x,x).shape
(4, 3, 6, 5)
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
(4, 3, 6, 5, 5)
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
(4, 3, 6, 6, 5, 5)

PyTorch

>>> import torch
>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kernel = KernelMatern12(d=2,torchify=True)
>>> x = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,2)))
>>> kernel(x[0],x[0]).item()
1.0
>>> kernel(x,x)
tensor([1., 1., 1., 1.], dtype=torch.float64, grad_fn=<MulBackward0>)
>>> kernel(x[:,None,:],x[None,:,:])
tensor([[1.0000, 0.6145, 0.7942, 0.6430],
        [0.6145, 1.0000, 0.5664, 0.5022],
        [0.7942, 0.5664, 1.0000, 0.8091],
        [0.6430, 0.5022, 0.8091, 1.0000]], dtype=torch.float64,
       grad_fn=<MulBackward0>)

Multiple randomizations

>>> x = torch.from_numpy(rng.uniform(low=0,high=1,size=(6,5,2)))
>>> kernel(x,x).shape
torch.Size([6, 5])
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
torch.Size([6, 5, 5])
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
torch.Size([6, 6, 5, 5])

Batch hyperparameters

>>> kernel = KernelMatern12(
...     d = 2,
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,1],
...     torchify = True)
>>> kernel(x,x).shape
torch.Size([4, 3, 6, 5])
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
torch.Size([4, 3, 6, 5, 5])
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
torch.Size([4, 3, 6, 6, 5, 5])

Source code in qmcpy/kernel/abstract_kernel.py

def __init__(
    self,
    d,
    scale=1.0,
    lengthscales=1.0,
    shape_scale=[1],
    shape_lengthscales=None,
    tfs_scale=(tf_exp_eps_inv, tf_exp_eps),
    tfs_lengthscales=(tf_exp_eps_inv, tf_exp_eps),
    torchify=False,
    requires_grad_scale=True,
    requires_grad_lengthscales=True,
    device="cpu",
    compile_call=False,
    comiple_call_kwargs={},
):
    r"""
    Args:
        d (int): Dimension.
        scale (Union[np.ndarray, torch.Tensor]): Scaling factor $S$.
        lengthscales (Union[np.ndarray, torch.Tensor]): Lengthscales $\boldsymbol{\gamma}$.
        shape_scale (list): Shape of `scale` when `np.isscalar(scale)`.
        shape_lengthscales (list): Shape of `lengthscales` when `np.isscalar(lengthscales)`
        tfs_scale (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_lengthscales (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        torchify (bool): If `True`, use the `torch` backend. Set to `True` if computing gradients with respect to inputs and/or hyperparameters.
        requires_grad_scale (bool): If `True` and `torchify`, set `requires_grad=True` for `scale`.
        requires_grad_lengthscales (bool): If `True` and `torchify`, set `requires_grad=True` for `lengthscales`.
        device (torch.device): If `torchify`, put things onto this device.
        compile_call (bool): If `True`, `torch.compile` the `parsed___call__` method.
        comiple_call_kwargs (dict): When `compile_call` is `True`, pass these keyword arguments to `torch.compile`.
    """
    super().__init__(
        d=d,
        torchify=torchify,
        device=device,
        compile_call=compile_call,
        comiple_call_kwargs=comiple_call_kwargs,
    )
    self.raw_scale = self.parse_assign_param(
        pname="scale",
        param=scale,
        shape_param=shape_scale,
        requires_grad_param=requires_grad_scale,
        tfs_param=tfs_scale,
        endsize_ops=[1],
        constraints=["POSITIVE"],
    )
    self.tfs_scale = tfs_scale
    self.batch_param_names.append("scale")
    self.raw_lengthscales = self.parse_assign_param(
        pname="lengthscales",
        param=lengthscales,
        shape_param=[self.d] if shape_lengthscales is None else shape_lengthscales,
        requires_grad_param=requires_grad_lengthscales,
        tfs_param=tfs_lengthscales,
        endsize_ops=[1, self.d],
        constraints=["POSITIVE"],
    )
    self.tfs_lengthscales = tfs_lengthscales
    self.batch_param_names.append("lengthscales")

KernelMatern32

Bases: AbstractKernelScaleLengthscales

Matern kernel with \(\alpha=3/2\).

\[K(\boldsymbol{x},\boldsymbol{z}) = S \left(1+\sqrt{3} d_{\boldsymbol{\gamma}}(\boldsymbol{x},\boldsymbol{z})\right)\exp\left(-\sqrt{3}d_{\boldsymbol{\gamma}}(\boldsymbol{x},\boldsymbol{z})\right), \qquad d_{\boldsymbol{\gamma}}(\boldsymbol{x},\boldsymbol{z}) = \left\lVert\frac{\boldsymbol{x}-\boldsymbol{z}}{\sqrt{2}\boldsymbol{\gamma}}\right\rVert_2.\]

Examples:

>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kernel = KernelMatern32(d=2)
>>> x = rng.uniform(low=0,high=1,size=(4,2))
>>> kernel(x[0],x[0]).item()
1.0
>>> kernel(x,x)
array([1., 1., 1., 1.])
>>> kernel(x[:,None,:],x[None,:,:])
array([[1.        , 0.79309639, 0.93871358, 0.82137958],
       [0.79309639, 1.        , 0.74137353, 0.66516872],
       [0.93871358, 0.74137353, 1.        , 0.9471166 ],
       [0.82137958, 0.66516872, 0.9471166 , 1.        ]])

Multiple randomizations

>>> x = rng.uniform(low=0,high=1,size=(6,5,2))
>>> kernel(x,x).shape
(6, 5)
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
(6, 5, 5)
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
(6, 6, 5, 5)

Batch hyperparameters

>>> kernel = KernelMatern32(
...     d = 2,
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,1])
>>> kernel(x,x).shape
(4, 3, 6, 5)
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
(4, 3, 6, 5, 5)
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
(4, 3, 6, 6, 5, 5)

PyTorch

>>> import torch
>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kernel = KernelMatern32(d=2,torchify=True)
>>> x = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,2)))
>>> kernel(x[0],x[0]).item()
1.0
>>> kernel(x,x)
tensor([1., 1., 1., 1.], dtype=torch.float64, grad_fn=<MulBackward0>)
>>> kernel(x[:,None,:],x[None,:,:])
tensor([[1.0000, 0.7931, 0.9387, 0.8214],
        [0.7931, 1.0000, 0.7414, 0.6652],
        [0.9387, 0.7414, 1.0000, 0.9471],
        [0.8214, 0.6652, 0.9471, 1.0000]], dtype=torch.float64,
       grad_fn=<MulBackward0>)

Multiple randomizations

>>> x = torch.from_numpy(rng.uniform(low=0,high=1,size=(6,5,2)))
>>> kernel(x,x).shape
torch.Size([6, 5])
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
torch.Size([6, 5, 5])
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
torch.Size([6, 6, 5, 5])

Batch hyperparameters

>>> kernel = KernelMatern32(
...     d = 2,
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,1],
...     torchify = True)
>>> kernel(x,x).shape
torch.Size([4, 3, 6, 5])
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
torch.Size([4, 3, 6, 5, 5])
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
torch.Size([4, 3, 6, 6, 5, 5])

Source code in qmcpy/kernel/abstract_kernel.py

def __init__(
    self,
    d,
    scale=1.0,
    lengthscales=1.0,
    shape_scale=[1],
    shape_lengthscales=None,
    tfs_scale=(tf_exp_eps_inv, tf_exp_eps),
    tfs_lengthscales=(tf_exp_eps_inv, tf_exp_eps),
    torchify=False,
    requires_grad_scale=True,
    requires_grad_lengthscales=True,
    device="cpu",
    compile_call=False,
    comiple_call_kwargs={},
):
    r"""
    Args:
        d (int): Dimension.
        scale (Union[np.ndarray, torch.Tensor]): Scaling factor $S$.
        lengthscales (Union[np.ndarray, torch.Tensor]): Lengthscales $\boldsymbol{\gamma}$.
        shape_scale (list): Shape of `scale` when `np.isscalar(scale)`.
        shape_lengthscales (list): Shape of `lengthscales` when `np.isscalar(lengthscales)`
        tfs_scale (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_lengthscales (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        torchify (bool): If `True`, use the `torch` backend. Set to `True` if computing gradients with respect to inputs and/or hyperparameters.
        requires_grad_scale (bool): If `True` and `torchify`, set `requires_grad=True` for `scale`.
        requires_grad_lengthscales (bool): If `True` and `torchify`, set `requires_grad=True` for `lengthscales`.
        device (torch.device): If `torchify`, put things onto this device.
        compile_call (bool): If `True`, `torch.compile` the `parsed___call__` method.
        comiple_call_kwargs (dict): When `compile_call` is `True`, pass these keyword arguments to `torch.compile`.
    """
    super().__init__(
        d=d,
        torchify=torchify,
        device=device,
        compile_call=compile_call,
        comiple_call_kwargs=comiple_call_kwargs,
    )
    self.raw_scale = self.parse_assign_param(
        pname="scale",
        param=scale,
        shape_param=shape_scale,
        requires_grad_param=requires_grad_scale,
        tfs_param=tfs_scale,
        endsize_ops=[1],
        constraints=["POSITIVE"],
    )
    self.tfs_scale = tfs_scale
    self.batch_param_names.append("scale")
    self.raw_lengthscales = self.parse_assign_param(
        pname="lengthscales",
        param=lengthscales,
        shape_param=[self.d] if shape_lengthscales is None else shape_lengthscales,
        requires_grad_param=requires_grad_lengthscales,
        tfs_param=tfs_lengthscales,
        endsize_ops=[1, self.d],
        constraints=["POSITIVE"],
    )
    self.tfs_lengthscales = tfs_lengthscales
    self.batch_param_names.append("lengthscales")

KernelMatern52

Bases: AbstractKernelScaleLengthscales

Matern kernel with \(\alpha=5/2\).

\[K(\boldsymbol{x},\boldsymbol{z}) = S \left(1+\sqrt{5} d_{\boldsymbol{\gamma}}(\boldsymbol{x},\boldsymbol{z}) + \frac{5}{3} d_{\boldsymbol{\gamma}}^2(\boldsymbol{x},\boldsymbol{z})\right)\exp\left(-\sqrt{5}d_{\boldsymbol{\gamma}}(\boldsymbol{x},\boldsymbol{z})\right), \qquad d_{\boldsymbol{\gamma}}(\boldsymbol{x},\boldsymbol{z}) = \left\lVert\frac{\boldsymbol{x}-\boldsymbol{z}}{\sqrt{2}\boldsymbol{\gamma}}\right\rVert_2.\]

Examples:

>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kernel = KernelMatern52(d=2)
>>> x = rng.uniform(low=0,high=1,size=(4,2))
>>> kernel(x[0],x[0]).item()
1.0
>>> kernel(x,x)
array([1., 1., 1., 1.])
>>> kernel(x[:,None,:],x[None,:,:])
array([[1.        , 0.83612941, 0.95801903, 0.861472  ],
       [0.83612941, 1.        , 0.78812397, 0.71396963],
       [0.95801903, 0.78812397, 1.        , 0.96425994],
       [0.861472  , 0.71396963, 0.96425994, 1.        ]])

Multiple randomizations

>>> x = rng.uniform(low=0,high=1,size=(6,5,2))
>>> kernel(x,x).shape
(6, 5)
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
(6, 5, 5)
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
(6, 6, 5, 5)

Batch hyperparameters

>>> kernel = KernelMatern52(
...     d = 2,
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,1])
>>> kernel(x,x).shape
(4, 3, 6, 5)
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
(4, 3, 6, 5, 5)
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
(4, 3, 6, 6, 5, 5)

PyTorch

>>> import torch
>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kernel = KernelMatern52(d=2,torchify=True)
>>> x = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,2)))
>>> kernel(x[0],x[0]).item()
1.0
>>> kernel(x,x)
tensor([1., 1., 1., 1.], dtype=torch.float64, grad_fn=<MulBackward0>)
>>> kernel(x[:,None,:],x[None,:,:])
tensor([[1.0000, 0.8361, 0.9580, 0.8615],
        [0.8361, 1.0000, 0.7881, 0.7140],
        [0.9580, 0.7881, 1.0000, 0.9643],
        [0.8615, 0.7140, 0.9643, 1.0000]], dtype=torch.float64,
       grad_fn=<MulBackward0>)

Multiple randomizations

>>> x = torch.from_numpy(rng.uniform(low=0,high=1,size=(6,5,2)))
>>> kernel(x,x).shape
torch.Size([6, 5])
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
torch.Size([6, 5, 5])
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
torch.Size([6, 6, 5, 5])

Batch hyperparameters

>>> kernel = KernelMatern52(
...     d = 2,
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,1],
...     torchify = True)
>>> kernel(x,x).shape
torch.Size([4, 3, 6, 5])
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
torch.Size([4, 3, 6, 5, 5])
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
torch.Size([4, 3, 6, 6, 5, 5])

Source code in qmcpy/kernel/abstract_kernel.py

def __init__(
    self,
    d,
    scale=1.0,
    lengthscales=1.0,
    shape_scale=[1],
    shape_lengthscales=None,
    tfs_scale=(tf_exp_eps_inv, tf_exp_eps),
    tfs_lengthscales=(tf_exp_eps_inv, tf_exp_eps),
    torchify=False,
    requires_grad_scale=True,
    requires_grad_lengthscales=True,
    device="cpu",
    compile_call=False,
    comiple_call_kwargs={},
):
    r"""
    Args:
        d (int): Dimension.
        scale (Union[np.ndarray, torch.Tensor]): Scaling factor $S$.
        lengthscales (Union[np.ndarray, torch.Tensor]): Lengthscales $\boldsymbol{\gamma}$.
        shape_scale (list): Shape of `scale` when `np.isscalar(scale)`.
        shape_lengthscales (list): Shape of `lengthscales` when `np.isscalar(lengthscales)`
        tfs_scale (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_lengthscales (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        torchify (bool): If `True`, use the `torch` backend. Set to `True` if computing gradients with respect to inputs and/or hyperparameters.
        requires_grad_scale (bool): If `True` and `torchify`, set `requires_grad=True` for `scale`.
        requires_grad_lengthscales (bool): If `True` and `torchify`, set `requires_grad=True` for `lengthscales`.
        device (torch.device): If `torchify`, put things onto this device.
        compile_call (bool): If `True`, `torch.compile` the `parsed___call__` method.
        comiple_call_kwargs (dict): When `compile_call` is `True`, pass these keyword arguments to `torch.compile`.
    """
    super().__init__(
        d=d,
        torchify=torchify,
        device=device,
        compile_call=compile_call,
        comiple_call_kwargs=comiple_call_kwargs,
    )
    self.raw_scale = self.parse_assign_param(
        pname="scale",
        param=scale,
        shape_param=shape_scale,
        requires_grad_param=requires_grad_scale,
        tfs_param=tfs_scale,
        endsize_ops=[1],
        constraints=["POSITIVE"],
    )
    self.tfs_scale = tfs_scale
    self.batch_param_names.append("scale")
    self.raw_lengthscales = self.parse_assign_param(
        pname="lengthscales",
        param=lengthscales,
        shape_param=[self.d] if shape_lengthscales is None else shape_lengthscales,
        requires_grad_param=requires_grad_lengthscales,
        tfs_param=tfs_lengthscales,
        endsize_ops=[1, self.d],
        constraints=["POSITIVE"],
    )
    self.tfs_lengthscales = tfs_lengthscales
    self.batch_param_names.append("lengthscales")

KernelMultiTask

Bases: AbstractKernel

Multi-task kernel

\[K((i,\boldsymbol{x}),(j,\boldsymbol{z})) = K_{\mathrm{task}}(i,j) K_{\mathrm{base}}(\boldsymbol{x},\boldsymbol{z})\]

parameterized for \(T\) tasks by a factor \(\mathsf{F} \in \mathbb{R}^{T \times r}\) and a diagonal \(\boldsymbol{v} \in \mathbb{R}^T\) so that

\[\left[K_{\mathrm{task}}(i,j)\right]_{i,j=1}^T = \mathsf{F} \mathsf{F}^T + \mathrm{diag}(\boldsymbol{v}).\]

Examples:

>>> kmt = KernelMultiTask(KernelGaussian(d=2),num_tasks=3,diag=[1,2,3])
>>> x = np.random.rand(4,2)
>>> task = np.arange(3)
>>> kmt(task[1],task[1],x[0],x[0]).item()
3.0
>>> kmt(task,task,x[0],x[0])
array([2., 3., 4.])
>>> kmt(task[:,None],task[None,:],x[0],x[0]).shape
(3, 3)
>>> kmt(task[1],task[1],x,x,)
array([3., 3., 3., 3.])
>>> kmt(task,task,x[:3],x[:3])
array([2., 3., 4.])
>>> kmt(task[:,None],task[:,None],x,x).shape
(3, 4)
>>> v = kmt(task,task,x[:3,None,:],x[None,:3,:])
>>> v.shape
(3, 3)
>>> np.allclose(v,kmt.base_kernel(x[:3,None,:],x[None,:3,:])*kmt.taskmat[task,task])
True
>>> kmt(task[:,None,None],task[:,None,None],x[:,None,:],x[None,:,:]).shape
(3, 4, 4)
>>> kmt(task[:,None,None,None],task[None,:,None,None],x[:,None,:],x[None,:,:]).shape
(3, 3, 4, 4)

Batched inference

>>> kernel_base = KernelGaussian(d=10,shape_lengthscales=(5,1),shape_scale=(3,5,1))
>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kmt = KernelMultiTask(kernel_base,
...     num_tasks=4,
...     shape_factor=(6,3,5,4,2),
...     diag = rng.uniform(low=0,high=1,size=(3,5,4)))
>>> x = np.random.rand(8,10)
>>> task = np.arange(4)
>>> kmt(task[0],task[0],x[0],x[0]).shape
(6, 3, 5)
>>> kmt(task,task,x[0],x[0]).shape
(6, 3, 5, 4)
>>> kmt(task[0],task[0],x,x).shape
(6, 3, 5, 8)
>>> v = kmt(task,task,x[:4,None,:],x[None,:4,:])
>>> v.shape
(6, 3, 5, 4, 4)
>>> kmat_x = kmt.base_kernel(x[:4,None,:],x[None,:4,:])
>>> kmat_x.shape
(3, 5, 4, 4)
>>> kmat_tasks = kmt.taskmat[...,task,task]
>>> kmat_tasks.shape
(6, 3, 5, 4)
>>> np.allclose(v,kmat_tasks[...,None,:]*kmat_x)
True
>>> np.allclose(v,kmat_tasks[...,:,None]*kmat_x)
False
>>> kmt(task[:,None,None,None],task[None,:,None,None],x[:,None,:],x[None,:,:]).shape
(6, 3, 5, 4, 4, 8, 8)
>>> kmt(task[:,None],task[None,:],x[:,None,None,None,:],x[None,:,None,None,:]).shape
(6, 3, 5, 8, 8, 4, 4)

Integrals

>>> kernel = KernelGaussian(
...     d = 2,
...     scale = rng.uniform(low=0,high=1,size=(3,1)),
...     lengthscales = rng.uniform(low=0,high=1,size=(1,2)))
>>> kmt = KernelMultiTask(
...     base_kernel=kernel,
...     num_tasks=5,
...     diag = rng.uniform(low=0,high=1,size=(6,3,5)),
...     factor = rng.uniform(low=0,high=1,size=(3,5,2)))
>>> task = np.arange(5)
>>> kmt.double_integral_01d(task0=task[0],task1=task[1]).shape
(6, 3)
>>> kmt.double_integral_01d(task0=task,task1=task).shape
(6, 3, 5)
>>> kmt.double_integral_01d(task0=task[:,None],task1=task[None,:]).shape
(6, 3, 5, 5)
>>> x = rng.uniform(low=0,high=1,size=(7,4,2))
>>> kmt.single_integral_01d(task0=task[0],task1=task[1],x=x).shape
(6, 3, 7, 4)
>>> kmt.single_integral_01d(task0=task[:,None,None],task1=task[:,None,None],x=x).shape
(6, 3, 5, 7, 4)
>>> kmt.single_integral_01d(task0=task[:,None,None,None],task1=task[None,:,None,None],x=x).shape
(6, 3, 5, 5, 7, 4)

Cholesky Construction

>>> kmt = KernelMultiTask(
...     KernelGaussian(5),
...     num_tasks = 3,
...     method = "CHOLESKY",
...     factor = 2,
...     diag = 1.5,
...     )
>>> kmt.taskmat
array([[ 2.25,  3.  ,  3.  ],
       [ 3.  ,  6.25,  7.  ],
       [ 3.  ,  7.  , 10.25]])
>>> kmt = KernelMultiTask(
...     KernelGaussian(5),
...     num_tasks = 3,
...     method = "CHOLESKY",
...     factor = 2,
...     diag = 1.5,
...     shape_factor = [2,3*(1+3)//2-3],
...     shape_diag = [3],
...     )
>>> kmt.taskmat
array([[[ 2.25,  3.  ,  3.  ],
        [ 3.  ,  6.25,  7.  ],
        [ 3.  ,  7.  , 10.25]],

       [[ 2.25,  3.  ,  3.  ],
        [ 3.  ,  6.25,  7.  ],
        [ 3.  ,  7.  , 10.25]]])
>>> kmt = KernelMultiTask(
...     KernelGaussian(5),
...     num_tasks = 3,
...     method = "CHOLESKY",
...     factor = 2,
...     diag = 1.5,
...     shape_factor = [3*(1+3)//2-3],
...     shape_diag = [2,3],
...     )
>>> kmt.taskmat
array([[[ 2.25,  3.  ,  3.  ],
        [ 3.  ,  6.25,  7.  ],
        [ 3.  ,  7.  , 10.25]],

       [[ 2.25,  3.  ,  3.  ],
        [ 3.  ,  6.25,  7.  ],
        [ 3.  ,  7.  , 10.25]]])
>>> kmt = KernelMultiTask(
...     KernelGaussian(5),
...     num_tasks = 3,
...     method = "CHOLESKY",
...     factor = 2,
...     diag = 1.5,
...     shape_factor = [2,3*(1+3)//2-3],
...     shape_diag = [2,3],
...     )
>>> kmt.taskmat
array([[[ 2.25,  3.  ,  3.  ],
        [ 3.  ,  6.25,  7.  ],
        [ 3.  ,  7.  , 10.25]],

       [[ 2.25,  3.  ,  3.  ],
        [ 3.  ,  6.25,  7.  ],
        [ 3.  ,  7.  , 10.25]]])

PyTorch

>>> import torch
>>> kmt = KernelMultiTask(KernelGaussian(d=2,torchify=True),num_tasks=3,diag=[1,2,3])
>>> x = torch.from_numpy(np.random.rand(4,2))
>>> task = np.arange(3)
>>> kmt(task[1],task[1],x[0],x[0]).item()
3.0
>>> kmt(task,task,x[0],x[0])
tensor([2., 3., 4.], dtype=torch.float64, grad_fn=<SelectBackward0>)
>>> kmt(task[:,None],task[None,:],x[0],x[0]).shape
torch.Size([3, 3])
>>> kmt(task[1],task[1],x,x)
tensor([3., 3., 3., 3.], dtype=torch.float64, grad_fn=<SelectBackward0>)
>>> kmt(task,task,x[:3],x[:3])
tensor([2., 3., 4.], dtype=torch.float64, grad_fn=<SelectBackward0>)
>>> kmt(task[:,None],task[:,None],x,x).shape
torch.Size([3, 4])
>>> v = kmt(task,task,x[:3,None,:],x[None,:3,:])
>>> v.shape
torch.Size([3, 3])
>>> torch.allclose(v,kmt.base_kernel(x[:3,None,:],x[None,:3,:])*kmt.taskmat[task,task])
True
>>> kmt(task[:,None,None],task[:,None,None],x[:,None,:],x[None,:,:]).shape
torch.Size([3, 4, 4])
>>> kmt(task[:,None,None,None],task[None,:,None,None],x[:,None,:],x[None,:,:]).shape
torch.Size([3, 3, 4, 4])

Batched inference

>>> kernel_base = KernelGaussian(d=10,shape_lengthscales=(5,1),shape_scale=(3,5,1),torchify=True)
>>> rng = np.random.Generator(np.random.PCG64(7))
>>> kmt = KernelMultiTask(kernel_base,
...     num_tasks=4,
...     shape_factor=(6,3,5,4,2),
...     diag = rng.uniform(low=0,high=1,size=(3,5,4)))
>>> x = torch.from_numpy(np.random.rand(8,10))
>>> task = torch.arange(4)
>>> kmt(task[0],task[0],x[0],x[0]).shape
torch.Size([6, 3, 5])
>>> kmt(task,task,x[0],x[0]).shape
torch.Size([6, 3, 5, 4])
>>> kmt(task[0],task[0],x,x).shape
torch.Size([6, 3, 5, 8])
>>> v = kmt(task,task,x[:4,None,:],x[None,:4,:])
>>> v.shape
torch.Size([6, 3, 5, 4, 4])
>>> kmat_x = kmt.base_kernel(x[:4,None,:],x[None,:4,:])
>>> kmat_x.shape
torch.Size([3, 5, 4, 4])
>>> kmat_tasks = kmt.taskmat[...,task,task]
>>> kmat_tasks.shape
torch.Size([6, 3, 5, 4])
>>> torch.allclose(v,kmat_tasks[...,None,:]*kmat_x)
True
>>> torch.allclose(v,kmat_tasks[...,:,None]*kmat_x)
False
>>> kmt(task[:,None,None,None],task[None,:,None,None],x[:,None,:],x[None,:,:]).shape
torch.Size([6, 3, 5, 4, 4, 8, 8])
>>> kmt(task[:,None],task[None,:],x[:,None,None,None,:],x[None,:,None,None,:]).shape
torch.Size([6, 3, 5, 8, 8, 4, 4])
>>> kmt.factor.dtype
torch.float32
>>> kmt.diag.dtype
torch.float64

Integrals

>>> kernel = KernelGaussian(
...     d = 2,
...     torchify = True,
...     scale = rng.uniform(low=0,high=1,size=(3,1)),
...     lengthscales = rng.uniform(low=0,high=1,size=(1,2)))
>>> kmt = KernelMultiTask(
...     base_kernel=kernel,
...     num_tasks=5,
...     diag = rng.uniform(low=0,high=1,size=(6,3,5)),
...     factor = rng.uniform(low=0,high=1,size=(3,5,2)))
>>> task = torch.arange(5)
>>> kmt.double_integral_01d(task0=task[0],task1=task[1]).shape
torch.Size([6, 3])
>>> kmt.double_integral_01d(task0=task,task1=task).shape
torch.Size([6, 3, 5])
>>> kmt.double_integral_01d(task0=task[:,None],task1=task[None,:]).shape
torch.Size([6, 3, 5, 5])
>>> x = torch.from_numpy(rng.uniform(low=0,high=1,size=(7,4,2)))
>>> kmt.single_integral_01d(task0=task[0],task1=task[1],x=x).shape
torch.Size([6, 3, 7, 4])
>>> kmt.single_integral_01d(task0=task[:,None,None],task1=task[:,None,None],x=x).shape
torch.Size([6, 3, 5, 7, 4])
>>> kmt.single_integral_01d(task0=task[:,None,None,None],task1=task[None,:,None,None],x=x).shape
torch.Size([6, 3, 5, 5, 7, 4])

Cholesky Construction

>>> kmt = KernelMultiTask(
...     KernelGaussian(5,torchify=True),
...     num_tasks = 3,
...     method = "CHOLESKY",
...     factor = 2,
...     diag = 1.5,
...     )
>>> kmt.taskmat
tensor([[ 2.2500,  3.0000,  3.0000],
        [ 3.0000,  6.2500,  7.0000],
        [ 3.0000,  7.0000, 10.2500]], grad_fn=<ViewBackward0>)
>>> kmt = KernelMultiTask(
...     KernelGaussian(5,torchify=True),
...     num_tasks = 3,
...     method = "CHOLESKY",
...     factor = 2,
...     diag = 1.5,
...     shape_factor = [2,3*(1+3)//2-3],
...     shape_diag = [3],
...     )
>>> kmt.taskmat
tensor([[[ 2.2500,  3.0000,  3.0000],
         [ 3.0000,  6.2500,  7.0000],
         [ 3.0000,  7.0000, 10.2500]],

        [[ 2.2500,  3.0000,  3.0000],
         [ 3.0000,  6.2500,  7.0000],
         [ 3.0000,  7.0000, 10.2500]]], grad_fn=<ViewBackward0>)
>>> kmt = KernelMultiTask(
...     KernelGaussian(5,torchify=True),
...     num_tasks = 3,
...     method = "CHOLESKY",
...     factor = 2,
...     diag = 1.5,
...     shape_factor = [3*(1+3)//2-3],
...     shape_diag = [2,3],
...     )
>>> kmt.taskmat
tensor([[[ 2.2500,  3.0000,  3.0000],
         [ 3.0000,  6.2500,  7.0000],
         [ 3.0000,  7.0000, 10.2500]],

        [[ 2.2500,  3.0000,  3.0000],
         [ 3.0000,  6.2500,  7.0000],
         [ 3.0000,  7.0000, 10.2500]]], grad_fn=<ViewBackward0>)
>>> kmt = KernelMultiTask(
...     KernelGaussian(5,torchify=True),
...     num_tasks = 3,
...     method = "CHOLESKY",
...     factor = 2,
...     diag = 1.5,
...     shape_factor = [2,3*(1+3)//2-3],
...     shape_diag = [2,3],
...     )
>>> kmt.taskmat
tensor([[[ 2.2500,  3.0000,  3.0000],
         [ 3.0000,  6.2500,  7.0000],
         [ 3.0000,  7.0000, 10.2500]],

        [[ 2.2500,  3.0000,  3.0000],
         [ 3.0000,  6.2500,  7.0000],
         [ 3.0000,  7.0000, 10.2500]]], grad_fn=<ViewBackward0>)

Parameters:

Name	Type	Description	Default
base_kernel	AbstractKernel	\(K_{\mathrm{base}}\).	required
num_tasks	int	Number of tasks \(T>1\).	required
factor	Union[ndarray, Tensor]	Factor \(\mathsf{F}\).	1.0
diag	Union[ndarray, Tensor]	Diagonal parameter \(\boldsymbol{v}\).	1.0
shape_factor	list	Shape of factor when np.isscalar(factor).	None
shape_diag	list	Shape of diag when np.isscalar(diag).	None
tfs_factor	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_identity, tf_identity)
tfs_diag	Tuple[callable, callable]	The first argument transforms to the raw value to be optimized; the second applies the inverse transform.	(tf_exp_eps_inv, tf_exp_eps)
requires_grad_factor	bool	If True and torchify, set requires_grad=True for factor.	True
requires_grad_diag	bool	If True and torchify, set requires_grad=True for diag.	True
method	str	"LOW RANK" or "CHOLESKY"	'LOW RANK'

Source code in qmcpy/kernel/multitask_kernel.py

def __init__(
    self,
    base_kernel,
    num_tasks,
    factor=1.0,
    diag=1.0,
    shape_factor=None,
    shape_diag=None,
    tfs_factor=(tf_identity, tf_identity),
    tfs_diag=(tf_exp_eps_inv, tf_exp_eps),
    requires_grad_factor=True,
    requires_grad_diag=True,
    rank_factor=1,
    method="LOW RANK",
):
    r"""
    Args:
        base_kernel (AbstractKernel): $K_{\mathrm{base}}$.
        num_tasks (int): Number of tasks $T>1$.
        factor (Union[np.ndarray, torch.Tensor]): Factor $\mathsf{F}$.
        diag (Union[np.ndarray, torch.Tensor]): Diagonal parameter $\boldsymbol{v}$.
        shape_factor (list): Shape of `factor` when `np.isscalar(factor)`.
        shape_diag (list): Shape of `diag` when `np.isscalar(diag)`.
        tfs_factor (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_diag (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        requires_grad_factor (bool): If `True` and `torchify`, set `requires_grad=True` for `factor`.
        requires_grad_diag (bool): If `True` and `torchify`, set `requires_grad=True` for `diag`.
        method (str): `"LOW RANK"` or "CHOLESKY"
    """
    assert isinstance(base_kernel, AbstractKernel)
    super().__init__(
        d=base_kernel.d,
        torchify=base_kernel.torchify,
        device=base_kernel.device,
        compile_call=False,
        comiple_call_kwargs={},
    )
    self.base_kernel = base_kernel
    self.AUTOGRADKERNEL = base_kernel.AUTOGRADKERNEL
    assert np.isscalar(num_tasks) and num_tasks % 1 == 0
    self.num_tasks = num_tasks
    assert (
        np.isscalar(rank_factor)
        and rank_factor % 1 == 0
        and 0 <= rank_factor <= self.num_tasks
    )
    self.method = str(method).upper().replace("_", " ").strip()
    if self.method == "LOW RANK":
        if shape_factor is None:
            shape_factor = [self.num_tasks, rank_factor]
        factor_endsize_ops = list(range(self.num_tasks + 1))
        diag_endsize_ops = [1, self.num_tasks]
    elif self.method == "CHOLESKY":
        if self.torchify:
            self.lti0, self.lti1 = self.npt.tril_indices(num_tasks, num_tasks, -1)
        else:
            self.lti0, self.lti1 = self.npt.tril_indices(num_tasks, -1)
        if shape_factor is None:
            shape_factor = [len(self.lti0)]
        factor_endsize_ops = [len(self.lti0)]
        diag_endsize_ops = [self.num_tasks]
    else:
        raise Exception(
            "invalid method = %s, must be in ['LOW RANK','CHOLESKY']" % self.method
        )
    self.raw_factor = self.parse_assign_param(
        pname="factor",
        param=factor,
        shape_param=shape_factor,
        requires_grad_param=requires_grad_factor,
        tfs_param=tfs_factor,
        endsize_ops=factor_endsize_ops,
        constraints=[],
    )
    self.tfs_factor = tfs_factor
    if self.method == "LOW RANK":
        assert self.raw_factor.shape[-2] == self.num_tasks
    self.raw_diag = self.parse_assign_param(
        pname="diag",
        param=diag,
        shape_param=[self.num_tasks] if shape_diag is None else shape_diag,
        requires_grad_param=requires_grad_diag,
        tfs_param=tfs_diag,
        endsize_ops=diag_endsize_ops,
        constraints=["NON-NEGATIVE"],
    )
    self.tfs_diag = tfs_diag
    self.eye_num_tasks = self.npt.eye(self.num_tasks, **self.nptkwargs)
    self.batch_param_names.append("taskmat")
    self._nbdim_base = None

__call__

__call__(task0, task1, x0, x1, beta0=None, beta1=None, c=None)

Evaluate the kernel with (optional) partial derivatives

\[\sum_{\ell=1}^p c_\ell \partial_{\boldsymbol{x}_0}^{\boldsymbol{\beta}_{\ell,0}} \partial_{\boldsymbol{x}_1}^{\boldsymbol{\beta}_{\ell,1}} K((i_0,\boldsymbol{x}_0),(i_1,\boldsymbol{x}_1)).\]

Parameters:

Name	Type	Description	Default
task0	Union[int, ndarray, Tensor]	First task indices \(i_0\).	required
task1	Union[int, ndarray, Tensor]	Second task indices \(i_1\).	required
x0	Union[ndarray, Tensor]	Shape x0.shape=(...,d) first input to kernel.	required
x1	Union[ndarray, Tensor]	Shape x1.shape=(...,d) second input to kernel.	required
beta0	Union[ndarray, Tensor]	Shape beta0.shape=(p,d) derivative orders with respect to first inputs, \(\boldsymbol{\beta}_0\).	None
beta1	Union[ndarray, Tensor]	Shape beta1.shape=(p,d) derivative orders with respect to first inputs, \(\boldsymbol{\beta}_1\).	None
c	Union[ndarray, Tensor]	Shape c.shape=(p,) coefficients of derivatives.	None

Returns:

Name	Type	Description
k	Union[ndarray, Tensor]	Kernel evaluations with batched shape, see the doctests for examples.

Source code in qmcpy/kernel/multitask_kernel.py

def __call__(self, task0, task1, x0, x1, beta0=None, beta1=None, c=None):
    r"""
    Evaluate the kernel with (optional) partial derivatives

    $$\sum_{\ell=1}^p c_\ell \partial_{\boldsymbol{x}_0}^{\boldsymbol{\beta}_{\ell,0}} \partial_{\boldsymbol{x}_1}^{\boldsymbol{\beta}_{\ell,1}} K((i_0,\boldsymbol{x}_0),(i_1,\boldsymbol{x}_1)).$$

    Args:
        task0 (Union[int, np.ndarray, torch.Tensor]): First task indices $i_0$.
        task1 (Union[int, np.ndarray, torch.Tensor]): Second task indices $i_1$.
        x0 (Union[np.ndarray, torch.Tensor]): Shape `x0.shape=(...,d)` first input to kernel.
        x1 (Union[np.ndarray, torch.Tensor]): Shape `x1.shape=(...,d)` second input to kernel.
        beta0 (Union[np.ndarray, torch.Tensor]): Shape `beta0.shape=(p,d)` derivative orders with respect to first inputs, $\boldsymbol{\beta}_0$.
        beta1 (Union[np.ndarray, torch.Tensor]): Shape `beta1.shape=(p,d)` derivative orders with respect to first inputs, $\boldsymbol{\beta}_1$.
        c (Union[np.ndarray, torch.Tensor]): Shape `c.shape=(p,)` coefficients of derivatives.

    Returns:
        k (Union[np.ndarray, torch.Tensor]): Kernel evaluations with batched shape, see the doctests for examples.
    """
    kmat_x = self.base_kernel.__call__(x0, x1, beta0, beta1, c)
    return self._parsed__call__(task0, task1, kmat_x)

single_integral_01d

single_integral_01d(task0, task1, x)

Evaluate the integral of the kernel over the unit cube

\[\tilde{K}((i_0,\boldsymbol{x}),i_1) = \int_{[0,1]^d} K((i_0,\boldsymbol{x}),(i_1,\boldsymbol{z}) \; \mathrm{d} \boldsymbol{z}.\]

Parameters:

Name	Type	Description	Default
task0	Union[int, ndarray, Tensor]	First task indices \(i_0\).	required
task1	Union[int, ndarray, Tensor]	Second task indices \(i_1\).	required
x	Union[ndarray, Tensor]	Shape x0.shape=(...,d) first input to kernel with	required

Returns:

Name	Type	Description
tildek	Union[ndarray, Tensor]	Shape y.shape=x.shape[:-1] integral kernel evaluations.

Source code in qmcpy/kernel/multitask_kernel.py

def single_integral_01d(self, task0, task1, x):
    r"""
    Evaluate the integral of the kernel over the unit cube

    $$\tilde{K}((i_0,\boldsymbol{x}),i_1) = \int_{[0,1]^d} K((i_0,\boldsymbol{x}),(i_1,\boldsymbol{z}) \; \mathrm{d} \boldsymbol{z}.$$

    Args:
        task0 (Union[int, np.ndarray, torch.Tensor]): First task indices $i_0$.
        task1 (Union[int, np.ndarray, torch.Tensor]): Second task indices $i_1$.
        x (Union[np.ndarray, torch.Tensor]): Shape `x0.shape=(...,d)` first input to kernel with

    Returns:
        tildek (Union[np.ndarray, torch.Tensor]): Shape `y.shape=x.shape[:-1]` integral kernel evaluations.
    """
    kint_x = self.base_kernel.single_integral_01d(x)
    return self._parsed__call__(task0, task1, kint_x)

double_integral_01d

double_integral_01d(task0, task1)

Evaluate the integral of the kernel over the unit cube

\[\tilde{K}(i_0,i_1) = \int_{[0,1]^d} \int_{[0,1]^d} K((i_0,\boldsymbol{x}),(i_1,\boldsymbol{z})) \; \mathrm{d} \boldsymbol{x} \; \mathrm{d} \boldsymbol{z}.\]

Parameters:

Name	Type	Description	Default
task0	Union[int, ndarray, Tensor]	First task indices \(i_0\).	required
task1	Union[int, ndarray, Tensor]	Second task indices \(i_1\).	required

Returns:

Name	Type	Description
tildek	Union[ndarray, Tensor]	Double integral kernel evaluations.

Source code in qmcpy/kernel/multitask_kernel.py

def double_integral_01d(self, task0, task1):
    r"""
    Evaluate the integral of the kernel over the unit cube

    $$\tilde{K}(i_0,i_1) = \int_{[0,1]^d} \int_{[0,1]^d} K((i_0,\boldsymbol{x}),(i_1,\boldsymbol{z})) \; \mathrm{d} \boldsymbol{x} \; \mathrm{d} \boldsymbol{z}.$$

    Args:
        task0 (Union[int, np.ndarray, torch.Tensor]): First task indices $i_0$.
        task1 (Union[int, np.ndarray, torch.Tensor]): Second task indices $i_1$.

    Returns:
        tildek (Union[np.ndarray, torch.Tensor]): Double integral kernel evaluations.
    """
    kint_x = self.base_kernel.double_integral_01d()
    return self._parsed__call__(task0, task1, kint_x)

UML Specific