sk_var_mult: Multiply a vector by a power of the covariance matrix

Description

Computes W %*% z, where z is the vector of non-NA data in g, and W is the pth power of V, the covariance matrix for z. By default, p=-1, so the function computes products with the inverse covariance matrix. Alternatively, out='quad' computes the quadratic form t(z) %*% W %*% z.

Usage

sk_var_mult(g, pars, fac_method = "eigen", fac = NULL, quad = FALSE, p = -1)

Value

numeric matrix

Arguments

g: a sk grid object or (if fac supplied) numeric vector or matrix of non-NA data
pars: list of form returned by sk_pars (with entries 'y', 'x', 'eps', 'psill')
fac_method: either 'eigen' (the default) or 'chol'
fac: factorization of scaled covariance matrix of z (V divided by psill)
quad: logical, if TRUE the function returns the quadratic form t(z) %*% V_inv %*% z
p: numeric, the matrix power of V^p to multiply (ignored when method=='chol')

Details

fac_method specifies the covariance matrix factorization to use: either 'chol' (Cholesky factorization), which only supports p=-1; or 'eigen' (eigen-decomposition), which supports any integer or numeric power for p. By default, the 'eigen' method is used, unless a Cholesky decomposition (matrix) is passed in fac.

The 'eigen' method is required when g has complete data (ie no NA values). Note that the structure of any supplied fac overrules the fac_method argument, so if your g is complete and you supply the Cholesky decomposition, the function will halt with an error even if you have set fac_method='eigen'.

factorization is often the slow part of the computations, so we allow it to be pre-computed using sk_var(..., scaled=TRUE) and passed in argument fac. This must be the factorization of the covariance matrix after dividing by the partial sill (see ?sk_var); either the lower Cholesky factor (eg the transposed output of base::chol), or a list of eigen-vectors and eigen-values (eg. the output of base::eigen). In the separable case, the eigen-decomposition is done on each of the x and y components separately, and fac should be a list with elements 'x' and 'y', each one a list of eigen-vectors and eigen-values.

When a factorization is supplied, all entries in pars, except for psill, are ignored, as they are baked into the factorization already. g can in this case be a numeric vector or matrix, containing one or more layers of observed data, with NAs omitted.

Examples

Run this code

# relative error comparing output x to reference y
rel_err = \(x, y) ifelse(y == 0, 0, abs( (x - y) / y ) )

# define example grid and data
gdim = c(10, 15)
g = sk(gdim)
n = length(g)
g[] = stats::rnorm(n)

# define covariance parameters
pars = utils::modifyList(sk_pars(g, 'gau'), list(psill=2, eps=0.5))

# COMPLETE CASE

# compute the full covariance matrix
V = sk_var(g, pars, sep=FALSE)
V_inv = chol2inv(chol(V))
out_reference = V_inv %*% g[]
out_reference_quad = t(g[]) %*% out_reference
max( rel_err(sk_var_mult(g, pars), out_reference) )
rel_err(sk_var_mult(g, pars, quad=TRUE), out_reference_quad)

# pre-computed factorization on separable components of correlation matrix
fac_corr = sk_var(utils::modifyList(g, list(gval=NULL)), pars, fac_method='eigen', sep=TRUE)
max( rel_err(sk_var_mult(g, pars, fac=fac_corr), out_reference) )
rel_err(sk_var_mult(g, pars, fac=fac_corr, quad=TRUE), out_reference_quad)

# matrix powers
out_reference = V %*% g[]
max( rel_err(sk_var_mult(g, pars, fac_method='eigen', p=1), out_reference) )
rel_err(sk_var_mult(g, pars, fac_method='eigen', p=1, quad=TRUE), t(g[]) %*% out_reference)

# INCOMPLETE CASE

n_sample = floor(n/10)
idx_sampled = sort(sample.int(n, n_sample))
g_miss = sk(gdim)
g_miss[idx_sampled] = g[idx_sampled]
V = sk_var(g_miss, pars)
sk_plot(V)

# correctness check (eigen used by default)
z = matrix(g[idx_sampled], ncol=1)
V_inv = chol2inv(chol(V))
out_reference = (V_inv %*% z)
out_reference_quad = t(z) %*% out_reference
max(rel_err(sk_var_mult(g_miss, pars), out_reference))
rel_err(sk_var_mult(g_miss, pars, quad=TRUE), out_reference_quad)

# check non-default Cholesky method
max( rel_err(sk_var_mult(g_miss, pars, fac_method='chol'), out_reference) )
rel_err(sk_var_mult(g_miss, pars, quad=TRUE, fac_method='chol'), out_reference_quad)

# supply data as a vector instead of list by pre-computing factorization
fac_chol = sk_var(g_miss, pars, scaled=TRUE, fac_method='chol')
fac_eigen = sk_var(g_miss, pars, scaled=TRUE, fac_method='eigen')
max(rel_err(sk_var_mult(z, pars, fac=fac_chol), out_reference))
max(rel_err(sk_var_mult(g_miss, pars, fac=fac_eigen), out_reference))
rel_err(sk_var_mult(z, pars, fac=fac_chol, quad=TRUE), out_reference_quad)
rel_err(sk_var_mult(g_miss, pars, fac=fac_eigen, quad=TRUE), out_reference_quad)

# matrix powers in eigen mode
out_reference = V %*% z
max(rel_err(sk_var_mult(g_miss, pars, p=1), out_reference))
rel_err(sk_var_mult(g_miss, pars, p=1, quad=TRUE), t(z) %*% out_reference)
max(rel_err(sk_var_mult(g_miss, pars, p=2), V %*% out_reference))

# verify that multiplying g_miss twice by a square root of V is same as multiplying by V
g_miss_sqrt = g_miss
g_miss_sqrt[!is.na(g_miss)] = sk_var_mult(g_miss, pars, p=1/2)
max( rel_err(sk_var_mult(g_miss_sqrt, pars, p=1/2), out_reference) )