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

In deankoch/pkern: Fast Kriging with Product Kernels

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

Usage

pkern_var_mult(
  g_obs,
  pars,
  fac_method = "eigen",
  fac = NULL,
  quad = FALSE,
  p = -1
)

Arguments

g_obs	list of form returned by pkern_grid or numeric vector or matrix of non-NA data
pars	list of form returned by pkern_pars (with entries 'y', 'x', 'eps', 'psill')
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

Alternatively, out='quad' computes the quadratic form t(z) %*% W %*% z.

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_obs has complete data (ie no NA values). Note that the structure of any supplied fac overrules the fac_method argument, so if your g_obs is complete and you supply the Cholesky decomposition, the function will throw an error.

As factorization is the slow part of the computations, it can be pre-computed using pkern_var(..., scaled=TRUE) and passed to pkern_var_mult in argument fac. This must be the factorization of the covariance matrix after dividing by the partial sill (see ?pkern_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_obs can in this case be a numeric vector or matrix, containing one or more layers of observed data (with NAs omitted).

Value

numeric matrix

Examples

# 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)
n = prod(gdim)
z_all = rnorm(n)
g_obs = modifyList(pkern_grid(gdim), list(gval = z_all))

# define covariance parameters
pars = modifyList(pkern_pars(g_obs, 'gau'), list(psill=2, eps=0.5))

# COMPLETE CASE

# compute the full covariance matrix
V = pkern_var(g_obs, pars, sep=FALSE)
V_inv = chol2inv(chol(V))
out_reference = V_inv %*% z_all
out_reference_quad = t(z_all) %*% out_reference
max( rel_err(pkern_var_mult(g_obs, pars), out_reference) )
rel_err(pkern_var_mult(g_obs, pars, quad=TRUE), out_reference_quad)

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

# matrix powers
out_reference = V %*% z_all
max( rel_err(pkern_var_mult(g_obs, pars, fac_method='eigen', p=1), out_reference) )
rel_err(pkern_var_mult(g_obs, pars, fac_method='eigen', p=1, quad=T), t(z_all) %*% out_reference)

# INCOMPLETE CASE

n_sample = floor(n/10)
idx_sampled = sample.int(n, n_sample) |> sort()
z_obs = rep(NA, n)
z_obs[idx_sampled] = z_all[idx_sampled]
g_obs = modifyList(g_obs, list(gval = z_obs))
V = pkern_var(g_obs, pars)
pkern_plot(V)

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

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

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

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

# multiply g_obs twice by a square root of V
g_obs_sqrt = g_obs
g_obs_sqrt$gval[!is.na(g_obs$gval)] = pkern_var_mult(g_obs, pars, p=1/2)
max( rel_err(pkern_var_mult(g_obs_sqrt, pars, p=1/2), out_reference) )

deankoch/pkern documentation built on Oct. 26, 2023, 8:54 p.m.