glearn.Covariance.solve#

Covariance.solve(Y, sigma=None, sigma0=None, scale=None, p=1, derivative=[])#

Solve linear system involving the powers of covariance matrix or its derivatives.

Parameters:
Ynumpy.ndarray

The right-hand side matrix of the linear system of equations. The size of this matrix is \(n \times m\) where \(n\) is the size of the covariance matrix.

sigmafloat, default=None

The hyperparameter \(\sigma\) of the covariance model where \(\sigma^2\) represents the variance of the correlated errors of the model. \(\sigma\) should be positive and cannot be None.

sigma0float, default=None

The hyperparameter \(\varsigma\) of the covariance model where \(\varsigma^2\) represents the variance of the input noise to of the model. \(\sigma\) should be positive and cannot be None.

scalefloat or array_like[float], default=None

The scale hyperparameters \(\boldsymbol{\alpha} = (\alpha_1, \dots, \alpha_d)\) in scales the distance between data points in \(\mathbb{R}^d\). If an array of the size \(d\) is given, each \(\alpha_i\) scales the distance in the \(i\)-th dimension. If a scalar value \(\alpha\) is given, all dimensions are scaled isometrically. \(\boldsymbol{\alpha}\) cannot be None.

pfloat, default=1

The integer exponent \(p\) (negative or positive) of the covariance matrix \(\boldsymbol{\Sigma}^{p}\) (see Notes below).

derivativelist, default=[]

Specifies a list of derivatives of covariance matrix with respect to the hyperparameters \(\boldsymbol{\alpha} = (\alpha_1, \dots, \alpha_d)\). A list of the size \(q\) with the components [i, j, ..., k] corresponds to take the derivative

\[\left. \frac{\partial^q}{\partial \alpha_{i+1} \partial \alpha_{j+1} \dots \partial \alpha_{k+1}} \boldsymbol{\Sigma}^{p}(\boldsymbol{\alpha} \vert \sigma^2, \varsigma^2) \right|_{\boldsymbol{\alpha}}.\]

Note

The derivative with respect to each hyperparameter \(\alpha_i\) can be at most of the order two, \(\partial^2 / \partial \alpha_i^2\). That is, each index in the derivative list can appear at most twice. For instance derivative=[1, 1] (second order derivative with respect to \(\alpha_{2}\)) is a valid input argument, how ever derivative=[1, 1, 1] (third order derivative) is an invalid input.

Note

When the derivative order is non-zero (meaning that derivative is not []), the exponent \(p\) should be 1.

Returns:
Xnumpy.ndarray

The solved array with the same size as of Y.

Notes

This function solves the linear system

\[\boldsymbol{\Sigma}^{p}_{(i, j, \dots, k)} \mathbf{X} = \mathbf{Y},\]

where \(\boldsymbol{\Sigma}^{p}_{(i, j, \dots, k)}\) is defined as

\[\boldsymbol{\Sigma}^{p}_{(i, j, \dots, k)} = \frac{\partial^q}{\partial \alpha_{i} \partial \alpha_{j} \dots \partial \alpha_{k}} \boldsymbol{\Sigma}^{p}(\boldsymbol{\alpha} \vert \sigma, \varsigma).\]

In the above, \(p\) is the matrix exponent and \(q\) is the order of derivation. Also, the covariance matrix \(\boldsymbol{\Sigma}\) is defined by

\[\boldsymbol{\Sigma}(\boldsymbol{\alpha}, \sigma, \varsigma) = \sigma^2 \mathbf{K}(\boldsymbol{\alpha}) + \varsigma^2 \mathbf{I}.\]

In the above, \(\mathbf{I}\) is the identity matrix and \(\mathbf{K}\) is the correlation matrix that depends on a set of scale hyperparameters \(\boldsymbol{\alpha}=(\alpha_1, \dots, \alpha_d)\).

Derivatives:

Note that the indices in list derivative=[i, j, ..., k] are zero-indexed, meaning that the index i corresponds to take derivative with respect to the hyperparameter \(\alpha_{i+1}\). For instance:

  • [] corresponds to no derivative.

  • [0] corresponds to \(\partial / \partial \alpha_1\) and [1] corresponds to \(\partial / \partial \alpha_2\).

  • [0, 2] corresponds to \(\partial^2 / \partial \alpha_1 \partial \alpha_3\).

  • [0, 0] corresponds to \(\partial^2 / \partial \alpha_1^2\).

  • [0, 2, 2, 4] corresponds to \(\partial^4 / \partial \alpha_1 \partial \alpha_{3}^2 \partial \alpha_5\).

Configuring Computation Settings:

This function passes the computation of the log-determinant to the function imate.logdet(). To configure the latter function, create a dictionary of input arguments to this function and pass the dictionary with glearn.Covariance.set_imate_options(). See examples below for details.

Examples

Basic Usage:

Create a covariance matrix based on a set of sample data with four points in \(d=2\) dimensional space.

>>> # Generate a set of points
>>> from glearn.sample_data import generate_points
>>> x = generate_points(num_points=4, dimension=2)

>>> # Create a covariance object
>>> from glearn import Covariance
>>> cov = Covariance(x)

In the following, we create a sample right-hand side matrix \(\mathbf{Y}\) of the size \(n \times 2\). The size of the covariance, \(n\), is also the same as the size of the number of points generated in the above. We solve the linear system

\[\boldsymbol{\Sigma}^{2} \mathbf{X} = \mathbf{Y},\]

for the hyperparameters \(\sigma=2\), \(\varsigma = 3\), and \(\boldsymbol{\alpha} = (1, 2)\).

>>> import numpy
>>> numpy.random.seed(0)
>>> n = cov.get_size()
>>> m = 2
>>> Y = numpy.random.randn(n, m)

>>> # Solve linear system.
>>> cov.solve(Y, sigma=2.0, sigma0=3.0, scale=[1.0, 2.0], p=2)
array([[ 0.00661562,  0.00021078],
       [-0.00194021,  0.02031901],
       [ 0.00776767, -0.0137298 ],
       [-0.00239916, -0.00373247]])

Taking Derivatives:

Solve the linear system involving the second mixed derivative

\[\boldsymbol{\Sigma}_{(1, 2)} \mathbf{X} = \mathbf{Y},\]

where here \(\boldsymbol{\Sigma}_{(1, 2)}\) is

\[\boldsymbol{\Sigma}_{(1, 2)} = \frac{\partial^2}{\partial \alpha_1 \partial \alpha_2} \boldsymbol{\Sigma} (\boldsymbol{\alpha} \vert \sigma, \varsigma).\]
>>> # Compute second mixed derivative
>>> cov.solve(Y, sigma=2.0, sigma0=3.0, scale=[1.0, 2.0], p=1,
...           derivative=[0, 1])
array([[  -64.59508396,    46.20009763],
       [   80.00325735,   -49.33936488],
       [-1320.94755293,   817.88667267],
       [  314.55331235,  -172.5169869 ]])