Fast Karhunen-Loève Expansions via FFT-Accelerated Toeplitz Operators

Abstract

Gaussian random fields are a versatile tool used in the fields of stochastic PDEs, uncertainty quantification, and geostatistical simulation. One way to obtain them is to use a truncated Karhunen-Loève expansion (KLE). Computing the expansion requires the leading eigenpairs of an $N\times N$ covariance matrix, where $N$ is the total number of grid cells. These are usually computed with a Krylov eigensolver, which relies on the covariance operator only within matrix-vector products. Stored densely, the matrix takes $\mathcal{O}(N^2)$ memory and each product $\mathcal{O}(N^2)$ time. For a stationary kernel on an equispaced grid, the covariance matrix becomes (block-) Toeplitz and the product evaluates in $\mathcal{O}(N\log N)$ time using FFT-based circulant embedding, without the need to assemble the dense matrix. In a matched single-threaded comparison, the median speedup of the eigensolve grows from $18\times$ at $N=4096$ to $183\times$ at $N=2^{15}$. This makes it possible to compute discretized fields that would otherwise be infeasible to compute in the standard formulation. We show that the same construction carries over to non-separable kernels as well as $d$ dimensions, using block-Toeplitz matrices. We extend it to piecewise-constant fields on arbitrary domains given as subsets of a tensor grid. Computational savings grow with problem size, and storage drops from $\mathcal{O}(N^2)$ to $\mathcal{O}(2^{d}N)$.