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, **kwargs)

Evaluate the kernel with (optional) partial derivatives

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

Parameters:

Name Type Description Default
x0 Union[ndarray, Tensor]

Shape x0.shape=(...,d) first input to kernel with

required
x1 Union[ndarray, Tensor]

Shape x1.shape=(...,d) second input to kernel with

required
beta0 Union[ndarray, Tensor]

Shape beta0.shape=(p,d) derivative orders with respect to first inputs, \(\boldsymbol{\beta}_0\).

 None
beta1 Union[ndarray, Tensor]

Shape beta1.shape=(p,d) derivative orders with respect to first inputs, \(\boldsymbol{\beta}_1\).

 None
c Union[ndarray, Tensor]

Shape c.shape=(p,) coefficients of derivatives.

 None
kwargs dict

keyword arguments to parsed call

 {}

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

Source code in qmcpy/kernel/abstract_kernel.py 
 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
157
def __call__(self, x0, x1, beta0=None, beta1=None, c=None, **kwargs):
    r"""
    Evaluate the kernel with (optional) partial derivatives 

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

    Args:
        x0 (Union[np.ndarray,torch.Tensor]): Shape `x0.shape=(...,d)` first input to kernel with 
        x1 (Union[np.ndarray,torch.Tensor]): Shape `x1.shape=(...,d)` second input to kernel with 
        beta0 (Union[np.ndarray,torch.Tensor]): Shape `beta0.shape=(p,d)` derivative orders with respect to first inputs, $\boldsymbol{\beta}_0$.
        beta1 (Union[np.ndarray,torch.Tensor]): Shape `beta1.shape=(p,d)` derivative orders with respect to first inputs, $\boldsymbol{\beta}_1$.
        c (Union[np.ndarray,torch.Tensor]): Shape `c.shape=(p,)` coefficients of derivatives.
        kwargs (dict): keyword arguments to parsed call
    Returns:
        k (Union[np.ndarray,torch.Tensor]): Shape `y.shape=(x0+x1).shape[:-1]` kernel evaluations. 
    """
    assert isinstance(x0,self.nptarraytype) 
    assert isinstance(x0,self.nptarraytype) 
    assert x0.shape[-1]==self.d, "the size of the last dimension of x0 must equal d=%d, got x0.shape=%s"%(self.d,str(tuple(x0.shape)))
    assert x1.shape[-1]==self.d, "the size of the last dimension of x1 must equal d=%d, got x1.shape=%s"%(self.d,str(tuple(x1.shape)))
    if beta0 is None:
        beta0 = self.npt.zeros((1,self.d),dtype=int,**self.nptkwargs)
    if beta1 is None:
        beta1 = self.npt.zeros((1,self.d),dtype=int,**self.nptkwargs)
    if not isinstance(beta0,self.nptarraytype):
        beta0 = self.nptarray(beta0)
    if not isinstance(beta1,self.nptarraytype):
        beta1 = self.nptarray(beta1)
    beta0 = self.npt.atleast_2d(beta0)
    beta1 = self.npt.atleast_2d(beta1)
    assert beta0.ndim==2 and beta1.ndim==2, "beta0 and beta1 must both be 2 dimensional"
    p = beta0.shape[0]
    assert beta0.shape==(p,self.d), "expected beta0.shape=(%d,%d) but got beta0.shape=%s"%(p,self.d,str(tuple(beta0.shape)))
    assert beta1.shape==(p,self.d), "expected beta1.shape=(%d,%d) but got beta1.shape=%s"%(p,self.d,str(tuple(beta1.shape)))
    assert (beta0%1==0).all() and (beta0>=0).all(), "require int beta0 >= 0"
    assert (beta1%1==0).all() and (beta1>=0).all(), "require int beta1 >= 0"
    if c is None:
        c = self.npt.ones(p,**self.nptkwargs)
    if not isinstance(c,self.nptarraytype):
        c = self.nptarray(c) 
    c = self.npt.atleast_1d(c) 
    assert c.shape==(p,), "expected c.shape=(%d,) but got c.shape=%s"%(p,str(tuple(c.shape)))
    if not self.AUTOGRADKERNEL:
        batch_params = self.get_batch_params(max(x0.ndim-1,x1.ndim-1))
        k = self.compiled_parsed___call__(x0,x1,beta0,beta1,c,batch_params,**kwargs)
    else:
        if (beta0==0).all() and (beta1==0).all():
            batch_params = self.get_batch_params(max(x0.ndim-1,x1.ndim-1))
            k = c.sum()*self.compiled_parsed___call__(x0,x1,batch_params,**kwargs)
        else: # requires autograd, so self.npt=torch
            assert self.torchify, "autograd requires torchify=True"
            import torch
            incoming_grad_enabled = torch.is_grad_enabled()
            torch.set_grad_enabled(True)
            incoming_grad_enabled_params = {pname: param.requires_grad for pname,param in self.named_parameters()}
            if not incoming_grad_enabled:
                for pname,param in self.named_parameters():
                    param.requires_grad_(False)
            if (beta0>0).any():
                tileshapex0 = tuple(self.npt.ceil(self.npt.tensor(x1.shape[:-1])/self.npt.tensor(x0.shape[:-1])).to(int))
                x0gs = [self.npt.tile(x0[...,j].clone().requires_grad_(True),tileshapex0) for j in range(self.d)]
                [x0gj.requires_grad_(True) for x0gj in x0gs]
                x0g = self.npt.stack(x0gs,dim=-1)
            else:
                x0g = x0
            if (beta1>0).any():
                tileshapex1 = tuple(self.npt.ceil(self.npt.tensor(x0.shape[:-1])/self.npt.tensor(x1.shape[:-1])).to(int))
                x1gs = [self.npt.tile(x1[...,j].clone().requires_grad_(True),tileshapex1) for j in range(self.d)]
                [x1gj.requires_grad_(True) for x1gj in x1gs]
                x1g = self.npt.stack(x1gs,dim=-1)
            else:
                x1g = x1
            batch_params = self.get_batch_params(max(x0.ndim-1,x1.ndim-1))
            k = 0.
            k_base = self.compiled_parsed___call__(x0g,x1g,batch_params,**kwargs)
            for l in range(p):
                if (beta0[l]>0).any() or (beta1[l]>0).any():
                    k_part = k_base.clone()
                    for j0 in range(self.d):
                        for _ in range(beta0[l,j0]):
                            k_part = torch.autograd.grad(k_part,x0gs[j0],grad_outputs=torch.ones_like(k_part,requires_grad=True),create_graph=True)[0]
                    for j1 in range(self.d):
                        for _ in range(beta1[l,j1]):
                            k_part = torch.autograd.grad(k_part,x1gs[j1],grad_outputs=torch.ones_like(k_part,requires_grad=True),create_graph=True)[0]
                else:
                    k_part = k_base 
                k += c[l]*k_part
            if not incoming_grad_enabled:
                for pname,param in self.named_parameters():
                    param.requires_grad_(incoming_grad_enabled_params[pname])
            if (not incoming_grad_enabled) or ((not any(incoming_grad_enabled_params.values())) and (not x0.requires_grad) and (not x1.requires_grad)):
                k = k.detach()
            torch.set_grad_enabled(incoming_grad_enabled)
    return k

single_integral_01d 

single_integral_01d(x)

Evaluate the integral of the kernel over the unit cube

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

Parameters:

Name Type Description Default
x Union[ndarray, Tensor]

Shape x0.shape=(...,d) first input to kernel with

required

Returns:

Name Type Description
tildek Union[ndarray, Tensor]

Shape y.shape=x.shape[:-1] integral kernel evaluations.
Source code in qmcpy/kernel/abstract_kernel.py 
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
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 
185
186
187
188
189
190
191
192
193
194
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 
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
270
271
272
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.parse_assign_param(
        pname = "scale",
        param = scale, 
        shape_param = shape_scale,
        requires_grad_param = requires_grad_scale,
        tfs_param = tfs_scale,
        endsize_ops = [1],
        constraints = ["POSITIVE"])
    self.tfs_scale = tfs_scale
    self.batch_param_names.append("scale")
    self.raw_lengthscales = self.parse_assign_param(
        pname = "lengthscales",
        param = lengthscales, 
        shape_param = [self.d] if shape_lengthscales is None else shape_lengthscales,
        requires_grad_param = requires_grad_lengthscales,
        tfs_param = tfs_lengthscales,
        endsize_ops = [1,self.d],
        constraints = ["POSITIVE"])
    self.tfs_lengthscales = tfs_lengthscales
    self.batch_param_names.append("lengthscales")

KernelShiftInvar

Bases: AbstractSIDSIKernel

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

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

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

Examples:

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

Batch Params 

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

Derivatives 

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

References:

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

Parameters:

Name Type Description Default
d int

Dimension.

required
scale Union[ndarray, Tensor]

Scaling factor \(S\).

 1.0
lengthscales Union[ndarray, Tensor]

Product weights \((\gamma_1,\dots,\gamma_d)\).

 None
alpha Union[ndarray, Tensor]

Smoothness parameters \((\alpha_1,\dots,\alpha_d)\) where \(\alpha_j \geq 1\) for \(j=1,\dots,d\).

 2
shape_scale list

Shape of scale when np.isscalar(scale).

[1]
shape_lengthscales list

Shape of lengthscales when np.isscalar(lengthscales)

 None
tfs_scale Tuple[callable, callable]

The first argument transforms to the raw value to be optimized; the second applies the inverse transform.

 (tf_exp_eps_inv, tf_exp_eps)
tfs_lengthscales Tuple[callable, callable]

The first argument transforms to the raw value to be optimized; the second applies the inverse transform.

 (tf_exp_eps_inv, tf_exp_eps)
torchify bool

If True, use the torch backend. Set to True if computing gradients with respect to inputs and/or hyperparameters.

 False
requires_grad_scale bool

If True and torchify, set requires_grad=True for scale.

 True
requires_grad_lengthscales bool

If True and torchify, set requires_grad=True for lengthscales.

 True
device device

If torchify, put things onto this device.

 'cpu'
compile_call bool

If True, torch.compile the parsed___call__ method.

False
comiple_call_kwargs dict

When compile_call is True, pass these keyword arguments to torch.compile.

 {}
Source code in qmcpy/kernel/si_dsi_kernels.py 
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
def __init__(self,
        d, 
        scale = 1., 
        lengthscales = None,
        alpha = 2,
        shape_scale = [1],
        shape_lengthscales = None, 
        tfs_scale = (tf_exp_eps_inv,tf_exp_eps),
        tfs_lengthscales = (tf_exp_eps_inv,tf_exp_eps),
        torchify = False, 
        requires_grad_scale = True, 
        requires_grad_lengthscales = True, 
        device = "cpu",
        compile_call = False,
        comiple_call_kwargs = {},
        ):
    r"""
    Args:
        d (int): Dimension. 
        scale (Union[np.ndarray,torch.Tensor]): Scaling factor $S$.
        lengthscales (Union[np.ndarray,torch.Tensor]): Product weights $(\gamma_1,\dots,\gamma_d)$.
        alpha (Union[np.ndarray,torch.Tensor]): Smoothness parameters $(\alpha_1,\dots,\alpha_d)$ where $\alpha_j \geq 1$ for $j=1,\dots,d$.
        shape_scale (list): Shape of `scale` when `np.isscalar(scale)`. 
        shape_lengthscales (list): Shape of `lengthscales` when `np.isscalar(lengthscales)`
        tfs_scale (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_lengthscales (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        torchify (bool): If `True`, use the `torch` backend. Set to `True` if computing gradients with respect to inputs and/or hyperparameters.
        requires_grad_scale (bool): If `True` and `torchify`, set `requires_grad=True` for `scale`.
        requires_grad_lengthscales (bool): If `True` and `torchify`, set `requires_grad=True` for `lengthscales`.
        device (torch.device): If `torchify`, put things onto this device.
        compile_call (bool): If `True`, `torch.compile` the `parsed___call__` method. 
        comiple_call_kwargs (dict): When `compile_call` is `True`, pass these keyword arguments to `torch.compile`.
    """
    super().__init__(
        d = d, 
        scale = scale, 
        lengthscales = lengthscales,
        alpha = alpha, 
        shape_alpha = [d], 
        alpha_endsize_ops = [d],
        shape_scale = shape_scale,
        shape_lengthscales = shape_lengthscales, 
        tfs_alpha = (tf_identity,tf_identity), 
        tfs_scale = tfs_scale,
        requires_grad_alpha = False,
        tfs_lengthscales = tfs_lengthscales, 
        torchify = torchify,
        requires_grad_scale = requires_grad_scale,
        requires_grad_lengthscales = requires_grad_lengthscales,
        device = device,
        compile_call = compile_call,
        comiple_call_kwargs = comiple_call_kwargs,
    )
    assert self.alpha.shape==(self.d,)
    assert all(int(alphaj) in BERNOULLIPOLYSDICT for alphaj in self.alpha)
    if self.torchify:
        import torch 
        self.lgamma = torch.lgamma 
    else:
        self.lgamma = scipy.special.loggamma

KernelShiftInvarCombined

Bases: AbstractSIDSIKernel

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

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

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

Examples:

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

Batch Params 

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

References:

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

Parameters:

Name Type Description Default
d int

Dimension.

required
scale Union[ndarray, Tensor]

Scaling factor \(S\).

 1.0
lengthscales Union[ndarray, Tensor]

Product weights \((\gamma_1,\dots,\gamma_d)\).

 None
alpha Union[ndarray, Tensor]

Weights \(\boldsymbol{\alpha}_1,\dots,\boldsymbol{\alpha}_d \in \mathbb{R}_{>0}^4\).

1
shape_scale list

Shape of scale when np.isscalar(scale).

[1]
shape_lengthscales list

Shape of lengthscales when np.isscalar(lengthscales)

 None
shape_alpha list

Shape of alpha when np.isscalar(alpha)

 None
tfs_scale Tuple[callable, callable]

The first argument transforms to the raw value to be optimized; the second applies the inverse transform.

 (tf_exp_eps_inv, tf_exp_eps)
tfs_lengthscales Tuple[callable, callable]

The first argument transforms to the raw value to be optimized; the second applies the inverse transform.

 (tf_exp_eps_inv, tf_exp_eps)
tfs_alpha Tuple[callable, callable]

The first argument transforms to the raw value to be optimized; the second applies the inverse transform.

 (tf_exp_eps_inv, tf_exp_eps)
torchify bool

If True, use the torch backend. Set to True if computing gradients with respect to inputs and/or hyperparameters.

 False
requires_grad_scale bool

If True and torchify, set requires_grad=True for scale.

 True
requires_grad_lengthscales bool

If True and torchify, set requires_grad=True for lengthscales.

 True
requires_grad_alpha bool

If True and torchify, set requires_grad=True for alpha.

 True
device device

If torchify, put things onto this device.

 'cpu'
compile_call bool

If True, torch.compile the parsed___call__ method.

False
comiple_call_kwargs dict

When compile_call is True, pass these keyword arguments to torch.compile.

 {}
Source code in qmcpy/kernel/si_dsi_kernels.py 
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
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
def __init__(self,
        d, 
        scale = 1., 
        lengthscales = None,
        alpha = 1,
        shape_scale = [1],
        shape_lengthscales = None, 
        shape_alpha = None,
        tfs_scale = (tf_exp_eps_inv,tf_exp_eps),
        tfs_lengthscales = (tf_exp_eps_inv,tf_exp_eps),
        tfs_alpha = (tf_exp_eps_inv,tf_exp_eps),
        torchify = False, 
        requires_grad_scale = True, 
        requires_grad_lengthscales = True, 
        requires_grad_alpha = True,
        device = "cpu",
        compile_call = False,
        comiple_call_kwargs = {},
        ):
    r"""
    Args:
        d (int): Dimension. 
        scale (Union[np.ndarray,torch.Tensor]): Scaling factor $S$.
        lengthscales (Union[np.ndarray,torch.Tensor]): Product weights $(\gamma_1,\dots,\gamma_d)$.
        alpha (Union[np.ndarray,torch.Tensor]): Weights $\boldsymbol{\alpha}_1,\dots,\boldsymbol{\alpha}_d \in \mathbb{R}_{>0}^4$. 
        shape_scale (list): Shape of `scale` when `np.isscalar(scale)`. 
        shape_lengthscales (list): Shape of `lengthscales` when `np.isscalar(lengthscales)`
        shape_alpha (list): Shape of `alpha` when `np.isscalar(alpha)`
        tfs_scale (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_lengthscales (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_alpha (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        torchify (bool): If `True`, use the `torch` backend. Set to `True` if computing gradients with respect to inputs and/or hyperparameters.
        requires_grad_scale (bool): If `True` and `torchify`, set `requires_grad=True` for `scale`.
        requires_grad_lengthscales (bool): If `True` and `torchify`, set `requires_grad=True` for `lengthscales`.
        requires_grad_alpha (bool): If `True` and `torchify`, set `requires_grad=True` for `alpha`.
        device (torch.device): If `torchify`, put things onto this device.
        compile_call (bool): If `True`, `torch.compile` the `parsed___call__` method. 
        comiple_call_kwargs (dict): When `compile_call` is `True`, pass these keyword arguments to `torch.compile`.
    """
    super().__init__(
        d = d, 
        scale = scale, 
        lengthscales = lengthscales,
        alpha = alpha, 
        alpha_endsize_ops = [d],
        shape_scale = shape_scale,
        shape_lengthscales = shape_lengthscales, 
        shape_alpha = [4,d] if shape_alpha is None else shape_alpha, 
        tfs_scale = tfs_scale,
        tfs_alpha = tfs_alpha, 
        tfs_lengthscales = tfs_lengthscales, 
        torchify = torchify,
        requires_grad_scale = requires_grad_scale,
        requires_grad_lengthscales = requires_grad_lengthscales,
        requires_grad_alpha = requires_grad_alpha,
        device = device,
        compile_call = compile_call,
        comiple_call_kwargs = comiple_call_kwargs,
    )
    assert self.alpha.shape[-2:]==(4,d)
    if self.torchify:
        import torch 
        self.lgamma = torch.lgamma 
    else:
        self.lgamma = scipy.special.loggamma
    a = self.npt.arange(1,5,**self.nptkwargs)
    self.coeffs = (-1)**(a+1)*self.npt.exp(2*a*np.log(2*np.pi)-self.lgamma(2*a+1))

KernelDigShiftInvar

Bases: AbstractSIDSIKernel

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

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

where

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

Examples:

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

Batch Params 

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

References:

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

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

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

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

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

Parameters:

Name Type Description Default
d int

Dimension.

required
t int

number of bits in binary represtnations. Typically dnb2.t where isinstance(dnb2,DigitalNetB2).

 None
scale Union[ndarray, Tensor]

Scaling factor \(S\).

 1.0
lengthscales Union[ndarray, Tensor]

Product weights \((\gamma_1,\dots,\gamma_d)\).

 None
alpha Union[ndarray, Tensor]

Smoothness parameters \((\alpha_1,\dots,\alpha_d)\) where \(\alpha_j \geq 1\) for \(j=1,\dots,d\).

 2
shape_scale list

Shape of scale when np.isscalar(scale).

[1]
shape_lengthscales list

Shape of lengthscales when np.isscalar(lengthscales)

 None
tfs_scale Tuple[callable, callable]

The first argument transforms to the raw value to be optimized; the second applies the inverse transform.

 (tf_exp_eps_inv, tf_exp_eps)
tfs_lengthscales Tuple[callable, callable]

The first argument transforms to the raw value to be optimized; the second applies the inverse transform.

 (tf_exp_eps_inv, tf_exp_eps)
torchify bool

If True, use the torch backend. Set to True if computing gradients with respect to inputs and/or hyperparameters.

 False
requires_grad_scale bool

If True and torchify, set requires_grad=True for scale.

 True
requires_grad_lengthscales bool

If True and torchify, set requires_grad=True for lengthscales.

 True
device device

If torchify, put things onto this device.

 'cpu'
compile_call bool

If True, torch.compile the parsed___call__ method.

False
comiple_call_kwargs dict

When compile_call is True, pass these keyword arguments to torch.compile.

 {}
Source code in qmcpy/kernel/si_dsi_kernels.py 
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
def __init__(self,
        d,
        t = None,
        scale = 1., 
        lengthscales = None,
        alpha = 2,
        shape_scale = [1],
        shape_lengthscales = None, 
        tfs_scale = (tf_exp_eps_inv,tf_exp_eps),
        tfs_lengthscales = (tf_exp_eps_inv,tf_exp_eps),
        torchify = False, 
        requires_grad_scale = True, 
        requires_grad_lengthscales = True, 
        device = "cpu",
        compile_call = False,
        comiple_call_kwargs = {},
        ):
    r"""
    Args:
        d (int): Dimension. 
        t (int): number of bits in binary represtnations. Typically `dnb2.t` where `isinstance(dnb2,DigitalNetB2)`.
        scale (Union[np.ndarray,torch.Tensor]): Scaling factor $S$.
        lengthscales (Union[np.ndarray,torch.Tensor]): Product weights $(\gamma_1,\dots,\gamma_d)$.
        alpha (Union[np.ndarray,torch.Tensor]): Smoothness parameters $(\alpha_1,\dots,\alpha_d)$ where $\alpha_j \geq 1$ for $j=1,\dots,d$.
        shape_scale (list): Shape of `scale` when `np.isscalar(scale)`. 
        shape_lengthscales (list): Shape of `lengthscales` when `np.isscalar(lengthscales)`
        tfs_scale (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_lengthscales (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        torchify (bool): If `True`, use the `torch` backend. Set to `True` if computing gradients with respect to inputs and/or hyperparameters.
        requires_grad_scale (bool): If `True` and `torchify`, set `requires_grad=True` for `scale`.
        requires_grad_lengthscales (bool): If `True` and `torchify`, set `requires_grad=True` for `lengthscales`.
        device (torch.device): If `torchify`, put things onto this device.
        compile_call (bool): If `True`, `torch.compile` the `parsed___call__` method. 
        comiple_call_kwargs (dict): When `compile_call` is `True`, pass these keyword arguments to `torch.compile`.
    """
    super().__init__(
        d = d, 
        scale = scale, 
        lengthscales = lengthscales,
        alpha = alpha,
        shape_alpha = [d], 
        alpha_endsize_ops = [d],
        shape_scale = shape_scale,
        shape_lengthscales = shape_lengthscales, 
        tfs_alpha = (tf_identity,tf_identity), 
        tfs_scale = tfs_scale, 
        tfs_lengthscales = tfs_lengthscales, 
        torchify = torchify,
        requires_grad_alpha = False,
        requires_grad_scale = requires_grad_scale,
        requires_grad_lengthscales = requires_grad_lengthscales,
        device = device,
        compile_call = compile_call,
        comiple_call_kwargs = comiple_call_kwargs,
    )
    assert self.alpha.shape==(self.d,)
    self.set_t(t)
    assert all(1<=int(alphaj)<=4 for alphaj in self.alpha)

KernelDigShiftInvarAdaptiveAlpha

Bases: AbstractSIDSIKernel

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

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

where

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

Examples:

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

Batch Params 

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

References:

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

Parameters:

Name Type Description Default
d int

Dimension.

required
t int

number of bits in binary represtnations. Typically dnb2.t where isinstance(dnb2,DigitalNetB2).

 None
scale Union[ndarray, Tensor]

Scaling factor \(S\).

 1.0
lengthscales Union[ndarray, Tensor]

Product weights \((\gamma_1,\dots,\gamma_d)\).

 None
alpha Union[ndarray, Tensor]

Smoothness parameters \((\alpha_1,\dots,\alpha_d)\) where \(\alpha_j \geq 1\) for \(j=1,\dots,d\).

 1
shape_alpha list

Shape of alpha when np.isscalar(alpha)

 None
shape_scale list

Shape of scale when np.isscalar(scale).

[1]
shape_lengthscales list

Shape of lengthscales when np.isscalar(lengthscales)

 None
tfs_scale Tuple[callable, callable]

The first argument transforms to the raw value to be optimized; the second applies the inverse transform.

 (tf_exp_eps_inv, tf_exp_eps)
tfs_lengthscales Tuple[callable, callable]

The first argument transforms to the raw value to be optimized; the second applies the inverse transform.

 (tf_exp_eps_inv, tf_exp_eps)
tfs_alpha Tuple[callable, callable]

The first argument transforms to the raw value to be optimized; the second applies the inverse transform.

 (tf_exp_eps_inv, tf_exp_eps)
torchify bool

If True, use the torch backend. Set to True if computing gradients with respect to inputs and/or hyperparameters.

 False
requires_grad_scale bool

If True and torchify, set requires_grad=True for scale.

 True
requires_grad_lengthscales bool

If True and torchify, set requires_grad=True for lengthscales.

 True
requires_grad_alpha bool

If True and torchify, set requires_grad=True for alpha.

 True
device device

If torchify, put things onto this device.

 'cpu'
compile_call bool

If True, torch.compile the parsed___call__ method.

False
comiple_call_kwargs dict

When compile_call is True, pass these keyword arguments to torch.compile.

 {}
Source code in qmcpy/kernel/si_dsi_kernels.py 
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
def __init__(self,
        d,
        t=None,
        scale = 1., 
        lengthscales = None,
        alpha = 1,
        shape_scale = [1],
        shape_lengthscales = None, 
        shape_alpha = None,
        tfs_scale = (tf_exp_eps_inv,tf_exp_eps),
        tfs_lengthscales = (tf_exp_eps_inv,tf_exp_eps),
        tfs_alpha = (tf_exp_eps_inv,tf_exp_eps),
        torchify = False, 
        requires_grad_scale = True, 
        requires_grad_lengthscales = True, 
        requires_grad_alpha = True, 
        device = "cpu",
        compile_call = False,
        comiple_call_kwargs = {},
        ):
    r"""
    Args:
        d (int): Dimension. 
        t (int): number of bits in binary represtnations. Typically `dnb2.t` where `isinstance(dnb2,DigitalNetB2)`.
        scale (Union[np.ndarray,torch.Tensor]): Scaling factor $S$.
        lengthscales (Union[np.ndarray,torch.Tensor]): Product weights $(\gamma_1,\dots,\gamma_d)$.
        alpha (Union[np.ndarray,torch.Tensor]): Smoothness parameters $(\alpha_1,\dots,\alpha_d)$ where $\alpha_j \geq 1$ for $j=1,\dots,d$.
        shape_alpha (list): Shape of `alpha` when `np.isscalar(alpha)`
        shape_scale (list): Shape of `scale` when `np.isscalar(scale)`. 
        shape_lengthscales (list): Shape of `lengthscales` when `np.isscalar(lengthscales)`
        tfs_scale (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_lengthscales (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_alpha (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        torchify (bool): If `True`, use the `torch` backend. Set to `True` if computing gradients with respect to inputs and/or hyperparameters.
        requires_grad_scale (bool): If `True` and `torchify`, set `requires_grad=True` for `scale`.
        requires_grad_lengthscales (bool): If `True` and `torchify`, set `requires_grad=True` for `lengthscales`.
        requires_grad_alpha (bool): If `True` and `torchify`, set `requires_grad=True` for `alpha`.
        device (torch.device): If `torchify`, put things onto this device.
        compile_call (bool): If `True`, `torch.compile` the `parsed___call__` method. 
        comiple_call_kwargs (dict): When `compile_call` is `True`, pass these keyword arguments to `torch.compile`.
    """
    super().__init__(
        d = d, 
        scale = scale, 
        lengthscales = lengthscales,
        alpha = alpha, 
        shape_alpha = [d] if shape_alpha is None else shape_alpha,
        alpha_endsize_ops = [d],
        shape_scale = shape_scale,
        shape_lengthscales = shape_lengthscales,
        tfs_alpha = tfs_alpha, 
        tfs_scale = tfs_scale, 
        tfs_lengthscales = tfs_lengthscales, 
        torchify = torchify,
        requires_grad_alpha = requires_grad_alpha,
        requires_grad_scale = requires_grad_scale,
        requires_grad_lengthscales = requires_grad_lengthscales,
        device = device,
        compile_call = compile_call,
        comiple_call_kwargs = comiple_call_kwargs,
    )
    self.set_t(t)
    self.batch_param_names.append("alpha")

KernelDigShiftInvarCombined

Bases: AbstractSIDSIKernel

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

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

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

Examples:

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

Batch Params 

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

Parameters:

Name Type Description Default
d int

Dimension.

required
t int

number of bits in binary represtnations. Typically dnb2.t where isinstance(dnb2,DigitalNetB2).

 None
scale Union[ndarray, Tensor]

Scaling factor \(S\).

 1.0
lengthscales Union[ndarray, Tensor]

Product weights \((\gamma_1,\dots,\gamma_d)\).

 None
alpha Union[ndarray, Tensor]

Weights \(\boldsymbol{\alpha}_1,\dots,\boldsymbol{\alpha}_d \in \mathbb{R}_{>0}^4\).

1.0
shape_scale list

Shape of scale when np.isscalar(scale).

[1]
shape_lengthscales list

Shape of lengthscales when np.isscalar(lengthscales)

 None
shape_alpha list

Shape of alpha when np.isscalar(alpha)

 None
tfs_scale Tuple[callable, callable]

The first argument transforms to the raw value to be optimized; the second applies the inverse transform.

 (tf_exp_eps_inv, tf_exp_eps)
tfs_lengthscales Tuple[callable, callable]

The first argument transforms to the raw value to be optimized; the second applies the inverse transform.

 (tf_exp_eps_inv, tf_exp_eps)
torchify bool

If True, use the torch backend. Set to True if computing gradients with respect to inputs and/or hyperparameters.

 False
requires_grad_scale bool

If True and torchify, set requires_grad=True for scale.

 True
requires_grad_lengthscales bool

If True and torchify, set requires_grad=True for lengthscales.

 True
requires_grad_alpha bool

If True and torchify, set requires_grad=True for alpha.

 True
device device

If torchify, put things onto this device.

 'cpu'
compile_call bool

If True, torch.compile the parsed___call__ method.

False
comiple_call_kwargs dict

When compile_call is True, pass these keyword arguments to torch.compile.

 {}
Source code in qmcpy/kernel/si_dsi_kernels.py 
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
def __init__(self,
        d,
        t = None,
        scale = 1., 
        lengthscales = None,
        alpha = 1.,
        shape_scale = [1],
        shape_lengthscales = None,
        shape_alpha = None,
        tfs_scale = (tf_exp_eps_inv,tf_exp_eps),
        tfs_lengthscales = (tf_exp_eps_inv,tf_exp_eps),
        tfs_alpha = (tf_exp_eps_inv,tf_exp_eps),
        torchify = False, 
        requires_grad_scale = True, 
        requires_grad_lengthscales = True,
        requires_grad_alpha = True,
        device = "cpu",
        compile_call = False,
        comiple_call_kwargs = {},
        ):
    r"""
    Args:
        d (int): Dimension. 
        t (int): number of bits in binary represtnations. Typically `dnb2.t` where `isinstance(dnb2,DigitalNetB2)`.
        scale (Union[np.ndarray,torch.Tensor]): Scaling factor $S$.
        lengthscales (Union[np.ndarray,torch.Tensor]): Product weights $(\gamma_1,\dots,\gamma_d)$.
        alpha (Union[np.ndarray,torch.Tensor]): Weights $\boldsymbol{\alpha}_1,\dots,\boldsymbol{\alpha}_d \in \mathbb{R}_{>0}^4$. 
        shape_scale (list): Shape of `scale` when `np.isscalar(scale)`. 
        shape_lengthscales (list): Shape of `lengthscales` when `np.isscalar(lengthscales)`
        shape_alpha (list): Shape of `alpha` when `np.isscalar(alpha)`
        tfs_scale (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_lengthscales (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        torchify (bool): If `True`, use the `torch` backend. Set to `True` if computing gradients with respect to inputs and/or hyperparameters.
        requires_grad_scale (bool): If `True` and `torchify`, set `requires_grad=True` for `scale`.
        requires_grad_lengthscales (bool): If `True` and `torchify`, set `requires_grad=True` for `lengthscales`.
        requires_grad_alpha (bool): If `True` and `torchify`, set `requires_grad=True` for `alpha`.
        device (torch.device): If `torchify`, put things onto this device.
        compile_call (bool): If `True`, `torch.compile` the `parsed___call__` method. 
        comiple_call_kwargs (dict): When `compile_call` is `True`, pass these keyword arguments to `torch.compile`.
    """
    super().__init__(
        d = d, 
        scale = scale, 
        lengthscales = lengthscales,
        alpha = alpha,
        shape_alpha = [4,d] if shape_alpha is None else shape_alpha, 
        alpha_endsize_ops = [d],
        shape_scale = shape_scale,
        shape_lengthscales = shape_lengthscales, 
        tfs_alpha = tfs_alpha, 
        tfs_scale = tfs_scale, 
        tfs_lengthscales = tfs_lengthscales, 
        torchify = torchify,
        requires_grad_alpha = requires_grad_alpha,
        requires_grad_scale = requires_grad_scale,
        requires_grad_lengthscales = requires_grad_lengthscales,
        device = device,
        compile_call = compile_call,
        comiple_call_kwargs = comiple_call_kwargs,
    )
    self.set_t(t)
    assert self.alpha.shape[-2:]==(4,d)

KernelGaussian

Bases: AbstractKernelGaussianSE

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

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

Examples:

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

Multiple randomizations 

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

Batch hyperparameters 

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

Integrals 

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

PyTorch

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

Multiple randomizations 

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

Batch hyperparameters 

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

Integrals 

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

Derivatives 

>>> kernel = KernelGaussian(
...     d = 3,
...     torchify = True,
...     scale = torch.from_numpy(rng.uniform(low=0,high=1,size=(1,))),
...     lengthscales = torch.from_numpy(rng.uniform(low=0,high=1,size=(3,))))
>>> x0 = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,))).requires_grad_(True)
>>> x1 = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,))).requires_grad_(True)
>>> x2 = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,))).requires_grad_(True)
>>> x = torch.stack([x0,x1,x2],axis=-1)
>>> z0 = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,))).requires_grad_(True)
>>> z1 = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,))).requires_grad_(True)
>>> z2 = torch.from_numpy(rng.uniform(low=0,high=1,size=(4,))).requires_grad_(True)
>>> z = torch.stack([z0,z1,z2],axis=-1)
>>> c = torch.from_numpy(rng.uniform(low=0,high=1,size=(2,)))
>>> beta0 = torch.tensor([
...     [1,0,0],
...     [0,2,0]])
>>> beta1 = torch.tensor([
...     [0,0,2],
...     [2,1,0]])
>>> with torch.no_grad():
...     y = kernel(x,z,beta0,beta1,c)
>>> y
tensor([ 97.6657,   3.8621, -65.9329,  -1.1932], dtype=torch.float64)
>>> y_no_deriv = kernel(x,z)
>>> y_first = y_no_deriv.clone()
>>> y_first = torch.autograd.grad(y_first,x0,grad_outputs=torch.ones_like(y_first,requires_grad=True),create_graph=True)[0]
>>> y_first = torch.autograd.grad(y_first,z2,grad_outputs=torch.ones_like(y_first,requires_grad=True),create_graph=True)[0]
>>> y_first = torch.autograd.grad(y_first,z2,grad_outputs=torch.ones_like(y_first,requires_grad=True),create_graph=True)[0]
>>> y_second = y_no_deriv.clone()
>>> y_second = torch.autograd.grad(y_second,x1,grad_outputs=torch.ones_like(y_second,requires_grad=True),create_graph=True)[0]
>>> y_second = torch.autograd.grad(y_second,x1,grad_outputs=torch.ones_like(y_second,requires_grad=True),create_graph=True)[0]
>>> y_second = torch.autograd.grad(y_second,z0,grad_outputs=torch.ones_like(y_second,requires_grad=True),create_graph=True)[0]
>>> y_second = torch.autograd.grad(y_second,z0,grad_outputs=torch.ones_like(y_second,requires_grad=True),create_graph=True)[0]
>>> y_second = torch.autograd.grad(y_second,z1,grad_outputs=torch.ones_like(y_second,requires_grad=True),create_graph=True)[0]
>>> yhat = (y_first*c[0]+y_second*c[1]).detach()
>>> yhat
tensor([ 97.6657,   3.8621, -65.9329,  -1.1932], dtype=torch.float64)
>>> torch.allclose(y,yhat)
True
Source code in qmcpy/kernel/abstract_kernel.py 
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
270
271
272
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.parse_assign_param(
        pname = "scale",
        param = scale, 
        shape_param = shape_scale,
        requires_grad_param = requires_grad_scale,
        tfs_param = tfs_scale,
        endsize_ops = [1],
        constraints = ["POSITIVE"])
    self.tfs_scale = tfs_scale
    self.batch_param_names.append("scale")
    self.raw_lengthscales = self.parse_assign_param(
        pname = "lengthscales",
        param = lengthscales, 
        shape_param = [self.d] if shape_lengthscales is None else shape_lengthscales,
        requires_grad_param = requires_grad_lengthscales,
        tfs_param = tfs_lengthscales,
        endsize_ops = [1,self.d],
        constraints = ["POSITIVE"])
    self.tfs_lengthscales = tfs_lengthscales
    self.batch_param_names.append("lengthscales")

KernelSquaredExponential

Bases: AbstractKernelGaussianSE

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

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

Examples:

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

Multiple randomizations 

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

Batch hyperparameters 

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

PyTorch

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

Multiple randomizations 

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

Batch hyperparameters 

>>> kernel = KernelSquaredExponential(
...     d = 2, 
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,1],
...     torchify = True)
>>> kernel(x,x).shape 
torch.Size([4, 3, 6, 5])
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
torch.Size([4, 3, 6, 5, 5])
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
torch.Size([4, 3, 6, 6, 5, 5])
Source code in qmcpy/kernel/abstract_kernel.py 
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
270
271
272
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.parse_assign_param(
        pname = "scale",
        param = scale, 
        shape_param = shape_scale,
        requires_grad_param = requires_grad_scale,
        tfs_param = tfs_scale,
        endsize_ops = [1],
        constraints = ["POSITIVE"])
    self.tfs_scale = tfs_scale
    self.batch_param_names.append("scale")
    self.raw_lengthscales = self.parse_assign_param(
        pname = "lengthscales",
        param = lengthscales, 
        shape_param = [self.d] if shape_lengthscales is None else shape_lengthscales,
        requires_grad_param = requires_grad_lengthscales,
        tfs_param = tfs_lengthscales,
        endsize_ops = [1,self.d],
        constraints = ["POSITIVE"])
    self.tfs_lengthscales = tfs_lengthscales
    self.batch_param_names.append("lengthscales")

KernelRationalQuadratic

Bases: AbstractKernelScaleLengthscales

Rational Quadratic kernel

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

Examples:

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

Multiple randomizations 

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

Batch hyperparameters 

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

PyTorch

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

Multiple randomizations 

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

Batch hyperparameters 

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

Parameters:

Name Type Description Default
d int

Dimension.

required
scale Union[ndarray, Tensor]

Scaling factor \(S\).

 1.0
lengthscales Union[ndarray, Tensor]

Lengthscales \(\boldsymbol{\gamma}\).

 1.0
alpha Union[ndarray, Tensor]

Scale mixture parameter \(\alpha\).

 1.0
shape_scale list

Shape of scale when np.isscalar(scale).

[1]
shape_lengthscales list

Shape of lengthscales when np.isscalar(lengthscales)

 None
shape_alpha list

Shape of alpha when np.isscalar(alpha)

 [1]
tfs_scale Tuple[callable, callable]

The first argument transforms to the raw value to be optimized; the second applies the inverse transform.

 (tf_exp_eps_inv, tf_exp_eps)
tfs_lengthscales Tuple[callable, callable]

The first argument transforms to the raw value to be optimized; the second applies the inverse transform.

 (tf_exp_eps_inv, tf_exp_eps)
tfs_alpha Tuple[callable, callable]

The first argument transforms to the raw value to be optimized; the second applies the inverse transform.

 (tf_exp_eps_inv, tf_exp_eps)
torchify bool

If True, use the torch backend. Set to True if computing gradients with respect to inputs and/or hyperparameters.

 False
requires_grad_scale bool

If True and torchify, set requires_grad=True for scale.

 True
requires_grad_lengthscales bool

If True and torchify, set requires_grad=True for lengthscales.

 True
requires_grad_alpha bool

If True and torchify, set requires_grad=True for alpha.

 True
device device

If torchify, put things onto this device.

 'cpu'
compile_call bool

If True, torch.compile the parsed___call__ method.

False
comiple_call_kwargs dict

When compile_call is True, pass these keyword arguments to torch.compile.

 {}
Source code in qmcpy/kernel/common_kernels.py 
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
507
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.parse_assign_param(
        pname = "alpha",
        param = alpha, 
        shape_param = shape_alpha,
        requires_grad_param = requires_grad_alpha,
        tfs_param = tfs_alpha,
        endsize_ops = [1],
        constraints = ["POSITIVE"])
    self.tfs_alpha = tfs_alpha
    self.batch_param_names.append("alpha")

KernelMatern12

Bases: AbstractKernelScaleLengthscales

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

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

Examples:

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

Multiple randomizations 

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

Batch hyperparameters 

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

PyTorch

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

Multiple randomizations 

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

Batch hyperparameters 

>>> kernel = KernelMatern12(
...     d = 2, 
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,1],
...     torchify = True)
>>> kernel(x,x).shape 
torch.Size([4, 3, 6, 5])
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
torch.Size([4, 3, 6, 5, 5])
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
torch.Size([4, 3, 6, 6, 5, 5])
Source code in qmcpy/kernel/abstract_kernel.py 
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
270
271
272
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.parse_assign_param(
        pname = "scale",
        param = scale, 
        shape_param = shape_scale,
        requires_grad_param = requires_grad_scale,
        tfs_param = tfs_scale,
        endsize_ops = [1],
        constraints = ["POSITIVE"])
    self.tfs_scale = tfs_scale
    self.batch_param_names.append("scale")
    self.raw_lengthscales = self.parse_assign_param(
        pname = "lengthscales",
        param = lengthscales, 
        shape_param = [self.d] if shape_lengthscales is None else shape_lengthscales,
        requires_grad_param = requires_grad_lengthscales,
        tfs_param = tfs_lengthscales,
        endsize_ops = [1,self.d],
        constraints = ["POSITIVE"])
    self.tfs_lengthscales = tfs_lengthscales
    self.batch_param_names.append("lengthscales")

KernelMatern32

Bases: AbstractKernelScaleLengthscales

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

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

Examples:

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

Multiple randomizations 

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

Batch hyperparameters 

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

PyTorch

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

Multiple randomizations 

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

Batch hyperparameters 

>>> kernel = KernelMatern32(
...     d = 2, 
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,1],
...     torchify = True)
>>> kernel(x,x).shape 
torch.Size([4, 3, 6, 5])
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
torch.Size([4, 3, 6, 5, 5])
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
torch.Size([4, 3, 6, 6, 5, 5])
Source code in qmcpy/kernel/abstract_kernel.py 
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
270
271
272
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.parse_assign_param(
        pname = "scale",
        param = scale, 
        shape_param = shape_scale,
        requires_grad_param = requires_grad_scale,
        tfs_param = tfs_scale,
        endsize_ops = [1],
        constraints = ["POSITIVE"])
    self.tfs_scale = tfs_scale
    self.batch_param_names.append("scale")
    self.raw_lengthscales = self.parse_assign_param(
        pname = "lengthscales",
        param = lengthscales, 
        shape_param = [self.d] if shape_lengthscales is None else shape_lengthscales,
        requires_grad_param = requires_grad_lengthscales,
        tfs_param = tfs_lengthscales,
        endsize_ops = [1,self.d],
        constraints = ["POSITIVE"])
    self.tfs_lengthscales = tfs_lengthscales
    self.batch_param_names.append("lengthscales")

KernelMatern52

Bases: AbstractKernelScaleLengthscales

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

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

Examples:

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

Multiple randomizations 

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

Batch hyperparameters 

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

PyTorch

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

Multiple randomizations 

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

Batch hyperparameters 

>>> kernel = KernelMatern52(
...     d = 2, 
...     shape_scale = [4,3,1],
...     shape_lengthscales = [3,1],
...     torchify = True)
>>> kernel(x,x).shape 
torch.Size([4, 3, 6, 5])
>>> kernel(x[:,:,None,:],x[:,None,:,:]).shape
torch.Size([4, 3, 6, 5, 5])
>>> kernel(x[:,None,:,None,:],x[None,:,None,:,:]).shape
torch.Size([4, 3, 6, 6, 5, 5])
Source code in qmcpy/kernel/abstract_kernel.py 
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
270
271
272
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.parse_assign_param(
        pname = "scale",
        param = scale, 
        shape_param = shape_scale,
        requires_grad_param = requires_grad_scale,
        tfs_param = tfs_scale,
        endsize_ops = [1],
        constraints = ["POSITIVE"])
    self.tfs_scale = tfs_scale
    self.batch_param_names.append("scale")
    self.raw_lengthscales = self.parse_assign_param(
        pname = "lengthscales",
        param = lengthscales, 
        shape_param = [self.d] if shape_lengthscales is None else shape_lengthscales,
        requires_grad_param = requires_grad_lengthscales,
        tfs_param = tfs_lengthscales,
        endsize_ops = [1,self.d],
        constraints = ["POSITIVE"])
    self.tfs_lengthscales = tfs_lengthscales
    self.batch_param_names.append("lengthscales")

KernelMultiTask

Bases: AbstractKernel

Multi-task kernel

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

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

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

Examples:

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

Batched inference

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

Integrals 

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

Cholesky Construction 

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

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

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

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

PyTorch

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

Batched inference

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

Integrals 

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

Cholesky Construction 

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

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

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

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

Parameters:

Name Type Description Default
base_kernel AbstractKernel

\(K_{\mathrm{base}}\).

required
num_tasks int

Number of tasks \(T>1\).

required
factor Union[ndarray, Tensor]

Factor \(\mathsf{F}\).

1.0
diag Union[ndarray, Tensor]

Diagonal parameter \(\boldsymbol{v}\).

1.0
shape_factor list

Shape of factor when np.isscalar(factor).

None
shape_diag list

Shape of diag when np.isscalar(diag).

None
tfs_factor Tuple[callable, callable]

The first argument transforms to the raw value to be optimized; the second applies the inverse transform.

 (tf_identity, tf_identity)
tfs_diag Tuple[callable, callable]

The first argument transforms to the raw value to be optimized; the second applies the inverse transform.

 (tf_exp_eps_inv, tf_exp_eps)
requires_grad_factor bool

If True and torchify, set requires_grad=True for factor.

 True
requires_grad_diag bool

If True and torchify, set requires_grad=True for diag.

 True
method str

"LOW RANK" or "CHOLESKY"

 'LOW RANK'
Source code in qmcpy/kernel/multitask_kernel.py 
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
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,
        method = "LOW RANK"
        ):
    r"""
    Args:
        base_kernel (AbstractKernel): $K_{\mathrm{base}}$. 
        num_tasks (int): Number of tasks $T>1$. 
        factor (Union[np.ndarray,torch.Tensor]): Factor $\mathsf{F}$. 
        diag (Union[np.ndarray,torch.Tensor]): Diagonal parameter $\boldsymbol{v}$. 
        shape_factor (list): Shape of `factor` when `np.isscalar(factor)`. 
        shape_diag (list): Shape of `diag` when `np.isscalar(diag)`. 
        tfs_factor (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        tfs_diag (Tuple[callable,callable]): The first argument transforms to the raw value to be optimized; the second applies the inverse transform.
        requires_grad_factor (bool): If `True` and `torchify`, set `requires_grad=True` for `factor`.
        requires_grad_diag (bool): If `True` and `torchify`, set `requires_grad=True` for `diag`.
        method (str): `"LOW RANK"` or "CHOLESKY"
    """
    assert isinstance(base_kernel,AbstractKernel)
    super().__init__(
        d = base_kernel.d,
        torchify = base_kernel.torchify,
        device = base_kernel.device,
        compile_call = False,
        comiple_call_kwargs = {})
    self.base_kernel = base_kernel
    self.AUTOGRADKERNEL = base_kernel.AUTOGRADKERNEL 
    assert np.isscalar(num_tasks) and num_tasks%1==0
    self.num_tasks = num_tasks
    assert np.isscalar(rank_factor) and rank_factor%1==0 and 0<=rank_factor<=self.num_tasks
    self.method = str(method).upper().replace("_"," ").strip()
    if self.method=="LOW RANK":
        if shape_factor is None:
            shape_factor = [self.num_tasks,rank_factor]
        factor_endsize_ops = list(range(self.num_tasks+1))
        diag_endsize_ops = [1,self.num_tasks]
    elif self.method=="CHOLESKY":
        if self.torchify:
            self.lti0,self.lti1 = self.npt.tril_indices(num_tasks,num_tasks,-1)
        else:
            self.lti0,self.lti1 = self.npt.tril_indices(num_tasks,-1)
        if shape_factor is None:
            shape_factor = [len(self.lti0)]
        factor_endsize_ops = [len(self.lti0)]
        diag_endsize_ops = [self.num_tasks]
    else:
        raise Exception("invalid method = %s, must be in ['LOW RANK','CHOLESKY']"%self.method)
    self.raw_factor = self.parse_assign_param(
        pname = "factor",
        param = factor,
        shape_param = shape_factor,
        requires_grad_param = requires_grad_factor,
        tfs_param = tfs_factor,
        endsize_ops = factor_endsize_ops,
        constraints = [])
    self.tfs_factor = tfs_factor
    if self.method=="LOW RANK":
        assert self.raw_factor.shape[-2]==self.num_tasks
    self.raw_diag = self.parse_assign_param(
        pname = "diag",
        param = diag, 
        shape_param = [self.num_tasks] if shape_diag is None else shape_diag,
        requires_grad_param = requires_grad_diag,
        tfs_param = tfs_diag,
        endsize_ops = diag_endsize_ops,
        constraints = ["NON-NEGATIVE"])
    self.tfs_diag = tfs_diag
    self.eye_num_tasks = self.npt.eye(self.num_tasks,**self.nptkwargs)
    self.batch_param_names.append("taskmat")
    self._nbdim_base = None

__call__ 

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

Evaluate the kernel with (optional) partial derivatives

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

Parameters:

Name Type Description Default
task0 Union[int, ndarray, Tensor]

First task indices \(i_0\).

required
task1 Union[int, ndarray, Tensor]

Second task indices \(i_1\).

 required
x0 Union[ndarray, Tensor]

Shape x0.shape=(...,d) first input to kernel.

 required
x1 Union[ndarray, Tensor]

Shape x1.shape=(...,d) second input to kernel.

required
beta0 Union[ndarray, Tensor]

Shape beta0.shape=(p,d) derivative orders with respect to first inputs, \(\boldsymbol{\beta}_0\).

 None
beta1 Union[ndarray, Tensor]

Shape beta1.shape=(p,d) derivative orders with respect to first inputs, \(\boldsymbol{\beta}_1\).

 None
c Union[ndarray, Tensor]

Shape c.shape=(p,) coefficients of derivatives.

 None

Returns:

Name Type Description
k Union[ndarray, Tensor]

Kernel evaluations with batched shape, see the doctests for examples.
Source code in qmcpy/kernel/multitask_kernel.py 
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
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 
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
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 
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
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