pkern_toep_mult: Efficiently compute yzx for symmetric Toeplitz matrices y and...

In deankoch/pkern: Fast Kriging with Product Kernels

Efficiently compute yzx for symmetric Toeplitz matrices y and x

Description

Computes the product y %*% z or y %*% z %*% x for symmetric Toeplitz matrices y and x and any numeric matrix z.

Usage

pkern_toep_mult(y, z = NULL, x = NULL, idx_obs = NULL, gdim = NULL)

Arguments

y	numeric matrix or vector, the symmetric Toeplitz matrix y or its first row
z	numeric matrix or vector with dimensionality conforming with y (and x)
x	numeric matrix or vector, the symmetric Toeplitz matrix x or its first row
idx_obs	integer vector, indices of the observed grid points

Details

Argument(s) y (and x) can be vector(s) supplying the first row of the matrix. By default, z is the identity matrix, so for matrix y, pkern_toep_mult(y) returns its argument, and for vector y, it returns the Toeplitz matrix generated by y, the same as base::toeplitz(y).

Fast Fourier transforms are used to reduce the memory footprint of computations, The first row(s) of y (and x) are embedded in a zero-padded vector representing a circulant matrix, whose action on the zero-padded version of z is equivalent to element-wise product in Fourier space. This allows the desired matrix product to be computed without explicitly creating matrices y or x in memory.

The function is optimized for grid data z that are sparse (many zeros). Before computing any transformations it first scans for and removes columns and rows of z which are all zero, replacing them afterwards.

To avoid unnecessarily copying large sparse matrices, z can be the vector of non-zero matrix entries only, where gdim specifies the full matrix dimensions and idx_obs the indices of the non-zero entries.

Value

numeric matrix, the product of yzx or yz (if x is NULL)

Examples

# define example matrix from 1D exponential variogram
n = 10
y = exp(1-seq(n))
y_mat = pkern_toep_mult(y)
max( abs(y_mat - stats::toeplitz(y))  )

# multiply by random matrix and compare with default matrix multiply
z = matrix(rnorm(n^2), n)
result_default = y_mat %*% z
max( abs( result_default - pkern_toep_mult(y_mat, z) ) )

# save memory by passing only the first row of the Toeplitz matrix
max( abs( result_default - pkern_toep_mult(y, z) ) )

# sparsify z and repeat
idx_sparse = sample.int(n^2, n^2 - n)
z[idx_sparse] = 0
result_default = y_mat %*% z
max( abs( result_default - pkern_toep_mult(y, z) ) )

# right-multiply with another kernel
x = exp( 2 *( 1-seq(n) ) )
x_mat = pkern_toep_mult(x)
result_default = result_default %*% x_mat
max( abs( result_default - pkern_toep_mult(y, z, x) ) )

# z can also be supplied as vector of nonzero grid values
idx_obs = which(z != 0)
gdim = c(y=n, x=n)
max( abs( result_default - pkern_toep_mult(y, z=z[idx_obs], x, idx_obs, gdim) ) )

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