Skip to content

Kernels

UML Overview

AbstractKernel

Bases: object

Source code in qmcpy/kernel/abstract_kernel.py 
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
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)

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

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

Source code in qmcpy/kernel/abstract_kernel.py 
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
def __call__(self, 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(\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.
    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)
    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)
        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.
            k_base = self.compiled_parsed___call__(x0g,x1g,batch_params)
            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 
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
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 
184
185
186
187
188
189
190
191
192
193
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 
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
def __init__(self,
        d, 
        scale = 1.,
        lengthscales = 1.,
        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.tf_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.batch_param_names.append("scale")
    self.raw_lengthscales,self.tf_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.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.])
>>> kernel_torch.single_integral_01d(xtorch)
tensor([10.], grad_fn=<SelectBackward0>)

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)
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
(4, 3, 6, 6, 5, 5)

Derivatives 

>>> 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)
>>> y
tensor([ 3003.7945, -1517.1556,   701.6692,  1470.4241], 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([ 3003.7946, -1517.1557,   701.6692,  1470.4241], dtype=torch.float64)
>>> 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())
>>> ynp
array([ 3003.79938927, -1517.15808005,   701.6700331 ,  1470.42580601])
>>> 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)\).

 1.0
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 
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
def __init__(self,
        d, 
        scale = 1., 
        lengthscales = 1.,
        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_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,
    )
    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

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) &= \sum_{k \in \mathbb{N}} \frac{\mathrm{wal}_k(x)}{2^{2 \lfloor \log_2(x) \rfloor}} = 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.])
>>> kernel_torch.single_integral_01d(xtorch)
tensor([10.], grad_fn=<SelectBackward0>)

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)
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
(4, 3, 6, 6, 5, 5)

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

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

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

  6. 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)\).

 1.0
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 
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
def __init__(self,
        d,
        t=None,
        scale = 1., 
        lengthscales = 1.,
        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_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.set_t(t)
    assert all(1<=int(alphaj)<=4 for alphaj in self.alpha)

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 
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
def __init__(self,
        d, 
        scale = 1.,
        lengthscales = 1.,
        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.tf_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.batch_param_names.append("scale")
    self.raw_lengthscales,self.tf_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.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 
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
def __init__(self,
        d, 
        scale = 1.,
        lengthscales = 1.,
        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.tf_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.batch_param_names.append("scale")
    self.raw_lengthscales,self.tf_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.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}\left\lVert\frac{\boldsymbol{x}-\boldsymbol{z}}{\sqrt{2}\boldsymbol{\gamma}}\right\rVert_2^2\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 
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
def __init__(self,
        d, 
        scale = 1., 
        lengthscales = 1.,
        alpha = 1.,
        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.tf_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.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 
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
def __init__(self,
        d, 
        scale = 1.,
        lengthscales = 1.,
        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.tf_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.batch_param_names.append("scale")
    self.raw_lengthscales,self.tf_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.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 
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
def __init__(self,
        d, 
        scale = 1.,
        lengthscales = 1.,
        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.tf_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.batch_param_names.append("scale")
    self.raw_lengthscales,self.tf_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.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 
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
def __init__(self,
        d, 
        scale = 1.,
        lengthscales = 1.,
        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.tf_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.batch_param_names.append("scale")
    self.raw_lengthscales,self.tf_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.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)

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

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
Source code in qmcpy/kernel/multitask_kernel.py 
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
def __init__(self,
        base_kernel,
        num_tasks, 
        factor = 1.,
        diag =  1.,
        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
        ):
    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`.
    """
    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.raw_factor,self.tf_factor = self.parse_assign_param(
        pname = "factor",
        param = factor,
        shape_param = [self.num_tasks,rank_factor] if shape_factor is None else shape_factor,
        requires_grad_param = requires_grad_factor,
        tfs_param = tfs_factor,
        endsize_ops = list(range(self.num_tasks+1)),
        constraints = [])
    assert self.raw_factor.shape[-2]==self.num_tasks
    self.raw_diag,self.tf_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 = [1,self.num_tasks],
        constraints = ["NON-NEGATIVE"])
    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 
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
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 
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
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 
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
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