Likelihood Evaluations Based on Fast Fourier Transforms

The module implements Fourier Methods for one-dimensional signals and two-dimensional signals. The functions of immediate use to practitioners are gp_rfft_lpdf to evaluate the likelihood of a one-dimensional signal and gp_inv_rfft to construct a non-centered parameterization using white noise. The two-dimensional analogues are gp_rfft2_lpdf and gp_inv_rfft2. The module also provides functions to evaluate kernels directly in the Fourier domain 1, including gp_periodic_exp_quad_cov_rfft and gp_periodic_matern_cov_rfft.

Utility Functions

The following functions, and their two-dimensional analogues, primarily exist as utility functions but may be useful for more complex models.

• gp_rfft transforms Gaussian process realizations to the Fourier domain and scales the coefficients such that they are white noise under the Gaussian process prior.

• gp_rfft_log_abs_det_jacobian evaluates the log absolute determinant of the Jacobian associated with the transformations gp_rfft.

Together, these two functions are used by gp_rfft_lpdf to evaluate the likelihood.

• gp_unpack_rfft unpacks complex Fourier coefficients of size n %/% 2 + 1 to a real vector of size n for easier manipulation.

• gp_pack_rfft is, unsurprisingly, the inverse of gp_unpack_rfft and packs a real vector of size n into a complex vector of size n %/% 2 + 1 such that the inverse RFFT can be applied.

• gp_evaluate_rfft_scale evaluates the expected standard deviation of Fourier coefficients obtained by transforming a Gaussian process. The values are arranged to match the output of gp_unpack_rfft.

Function Reference

One-dimensional Signals

vector gp_evaluate_rfft_scale(vector cov_rfft, int n)

Evaluate the scale of Fourier coefficients.

Parameters

• cov_rfft – Precomputed real fast Fourier transform of the kernel with shape (..., n %/% 2 + 1).

• n – Size of the real signal. Necessary because the size cannot be inferred from cov_rfft.

Returns

Scale of Fourier coefficients with shape (..., n %/% 2 + 1).

vector gp_unpack_rfft(complex_vector x, int n)

Unpack the Fourier coefficients of a real Fourier transform with n %/% 2 + 1 elements to a vector of n elements.

Parameters

• z – Real Fourier transform coefficients.

• n – Size of the real signal. Necessary because the size cannot be inferred from rfft.

Returns

Unpacked vector of n elements comprising the n %/% 2 + 1 real parts of the zero frequency term, complex terms, and Nyqvist frequency term (for even n). The subsequent (n - 1) %/% 2 elements are the imaginary parts of complex coefficients.

vector gp_rfft(vector y, vector loc, vector cov_rfft)

Transform a Gaussian process realization to white noise in the Fourier domain.

Parameters

• y – Realization of the Gaussian process with shape (..., n).

• loc – Mean of the Gaussian process with shape (..., n).

• cov_rfft – Precomputed real fast Fourier transform of the kernel with shape (..., n // 2 + 1).

Returns

Fourier-domain white noise with shape (..., n). See gp_unpack_rfft for details on the data structure.

real gp_rfft_log_abs_det_jacobian(vector cov_rfft, int n)

Evaluate the log absolute determinant of the Jacobian associated with gp_rfft.

Parameters

• cov_rfft – Precomputed real fast Fourier transform of the kernel with shape (..., n %/% 2 + 1).

• n – Size of the real signal. Necessary because the size cannot be inferred from rfft.

Returns

Log absolute determinant of the Jacobian.

real gp_rfft_lpdf(vector y, vector loc, vector cov_rfft)

Evaluate the log probability of a one-dimensional Gaussian process realization in Fourier space.

Parameters

• y – Realization of a Gaussian process with n grid points.

• loc – Mean of the Gaussian process of size n.

• cov_rfft – Precomputed real fast Fourier transform of the kernel of size n %/% 2 + 1.

Returns

Log probability of the Gaussian process realization.

complex_vector gp_pack_rfft(vector z)

Transform a real vector with n elements to a vector of complex Fourier coefficients with n elements ready for inverse real fast Fourier transformation.

vector gp_inv_rfft(vector z, vector loc, vector cov_rfft)

Transform white noise in the Fourier domain to a Gaussian process realization, i.e., a non-centered* parameterization in the Fourier domain.

The n real white noise variables must be assembled into a length-n %/% 2 + 1 complex vector with structure expected by the fast Fourier transform. The input vector z comprises

• the real zero frequency term,

• (n - 1) %/% 2 real parts of positive frequency terms,

• the real Nyqvist frequency term if n is even,

• and (n - 1) %/% 2 imaginary parts of positive frequency terms.

Parameters

• z – Fourier-domain white noise comprising n elements.

• loc – Mean of the Gaussian process.

• cov_rfft – Real fast Fourier transform of the covariance kernel.

Returns

Realization of the Gaussian process with n elements.

vector gp_periodic_exp_quad_cov_rfft(int n, real sigma, real length_scale, real period)

Evaluate the real fast Fourier transform of the periodic squared exponential kernel.

Parameters

• n – Number of grid points.

• sigma – Scale of the covariance.

• length_scale – Correlation length.

• period – Period for circular boundary conditions.

Returns

Fourier transform of the squared exponential kernel of size n %/% 2 + 1.

vector gp_periodic_matern_cov_rfft(real nu, int n, real sigma, real length_scale, real period)

Evaluate the real fast Fourier transform of the periodic Matern kernel.

Parameters

• nu – Smoothness parameter (1 / 2, 3 / 2, and 5 / 2 are typical values).

• n – Number of grid points.

• sigma – Scale of the covariance.

• length_scale – Correlation length.

• period – Period for circular boundary conditions.

Returns

Fourier transform of the squared exponential kernel of size n %/% 2 + 1.

Two-dimensional Signals

matrix gp_evaluate_rfft2_scale(matrix cov_rfft2, int width)

Evaluate the scale of Fourier coefficients.

Parameters

• cov_rfft2 – Precomputed real fast Fourier transform of the kernel with shape (height, width %/% 2 + 1).

• width – Number of columns of the signal (cannot be inferred from the Fourier coefficients).

Returns

Scale of Fourier coefficients with shape (height, width %/% 2 + 1).

matrix gp_unpack_rfft2(complex_matrix z, int m)

Unpack the complex Fourier coefficients of a two-dimensional real Fourier transform with shape to a real matrix with shape (height, width).

TODO: add details on packing structure.

Parameters

• z – Two-dimensional real Fourier transform coefficients with shape (height, width %/% 2 + 1).

• width – Number of columns of the signal (cannot be inferred from the Fourier coefficients z).

Returns

Unpacked matrix with shape (height, width).

complex_matrix gp_pack_rfft2(matrix z)

Transform a real matrix with shape (height, width) to a matrix of complex Fourier coefficients with shape (height, width %/% 2 + 1) ready for inverse real fast Fourier transformation in two dimensions.

Parameters

z – Unpacked matrices with shape (height, width).

Returns

Two-dimensional real Fourier transform coefficients.

matrix gp_rfft2(matrix y, matrix loc, matrix cov_rfft2)

Transform a Gaussian process realization to white noise in the Fourier domain.

Parameters

• y – Realization of the Gaussian process with shape (height, width).

• loc – Mean of the Gaussian process with shape (height, width).

• cov_rfft2 – Precomputed real fast Fourier transform of the kernel with shape (height, width %/% 2 + 1).

Returns

Unpacked matrix with shape (height, width).

matrix gp_inv_rfft2(matrix z, matrix loc, matrix cov_rfft2)

Transform white noise in the Fourier domain to a Gaussian process realization.

Parameters

• z – Unpacked matrix with shape (height, width).

• loc – Mean of the Gaussian process with shape (height, width).

• cov_rfft2 – Precomputed real fast Fourier transform of the kernel with shape (height, width %/% 2 + 1).

Returns

Realization of the Gaussian process.

real gp_rfft2_log_abs_det_jacobian(matrix cov_rfft2, int width)

Evaluate the log absolute determinant of the Jacobian associated with gp_rfft2.

Parameters

• cov_rfft2 – Precomputed real fast Fourier transform of the kernel with shape (height, width %/% 2 + 1).

• width – Number of columns of the signal (cannot be inferred from the precomputed kernel Fourier transform cov_rfft2).

Returns

Log absolute determinant of the Jacobian.

real gp_rfft2_lpdf(matrix y, matrix loc, matrix cov_rfft2)

Evaluate the log probability of a two-dimensional Gaussian process realization in Fourier space.

Parameters

• y – Realization of a Gaussian process with shape (height, width), where height is the number of rows, and width is the number of columns.

• loc – Mean of the Gaussian process with shape (height, width).

• cov_rfft2 – Precomputed real fast Fourier transform of the kernel with shape (height, width %/% 2 + 1).

Returns

Log probability of the Gaussian process realization.

matrix gp_periodic_exp_quad_cov_rfft2(int height, int width, real sigma, vector length_scale, vector period)

Evaluate the two-dimensional real fast Fourier transform of the periodic squared exponential kernel.

Parameters

• height – Number of grid rows.

• width – Number of grid columns.

• sigma – Scale of the covariance.

• length_scale – Correlation lengths in each dimension.

• period – Periods for circular boundary conditions in each dimension.

Returns

Fourier transform of the periodic squared exponential kernel with shape (height, width %/% 2 + 1).

matrix gp_periodic_matern_cov_rfft2(real nu, int height, int width, real sigma, vector length_scale, vector period)

Evaluate the two-dimensional real fast Fourier transform of the periodic Matern kernel.

Parameters

• nu – Smoothness parameter (1 / 2, 3 / 2, and 5 / 2 are typical values).

• height – Number of grid rows.

• width – Number of grid columns.

• sigma – Scale of the covariance.

• length_scale – Correlation lengths in each dimension.

• period – Periods for circular boundary conditions in each dimension.

Returns

Fourier transform of the periodic Matern kernel with shape (height, width %/% 2 + 1).

1

Fourier-domain kernels are implemented by discretizing their power spectrum naively. This approach works well if the number of grid points is large and the correlation length is small compared with the size of the domain. More sophisticated techniques may be required otherwise.