Fast nonparametric spectral density estimation from irregularly sampled data

TL;DR

The paper tackles nonparametric spectral density estimation for processes observed at fully irregular locations by introducing a weighted nonuniform Fourier estimator hat{S}(ξ) whose weights are designed so that H_{α}(ω) closely matches the Fourier transform G(ω) of a window g. This quadrature-based design reduces aliasing, avoids ill-conditioning inherent in nonuniform inverse problems, and remains scalable via NUFFT-based computations; the authors provide a rigorous bias analysis and demonstrate effectiveness on large-scale 1D and 2D data, including up to a million observation points. A key contribution is the explicit framework for window selection (favoring highly concentrated prolate/Slepian functions with Kaiser as a practical alternative), together with a method to compute weights α by solving a linear system with Chebyshev-based sampling, supported by stability and rate guarantees under oversampling. The approach yields substantial improvements over Lomb-Scargle and related methods, as evidenced by theoretical aliasing bounds and numerical demonstrations, highlighting its potential for accurate, scalable spectral analysis in irregularly sampled, high-dimensional settings.

Abstract

We introduce a nonparametric spectral density estimator for continuous-time and continuous-space processes measured at fully irregular locations. Our estimator is constructed using a weighted nonuniform Fourier sum whose weights yield a high-accuracy quadrature rule with respect to a user-specified window function. The resulting estimator significantly reduces the aliasing seen in periodogram approaches and least squares spectral analysis, sidesteps the dangers of ill-conditioning of the nonuniform Fourier inverse problem, and can be adapted to a wide variety of irregular sampling settings. We describe methods for rapidly computing the necessary weights in various settings, making the estimator scalable to large datasets. We then provide a theoretical analysis of sources of bias, and close with demonstrations of the method's efficacy, including for processes that exhibit very slow spectral decay and are observed at up to a million locations in multiple dimensions.

Figures (10)

• Figure 1: A comparison of two existing estimators and the estimator proposed in this work, each estimating the spectral density $S$ of a process $Y(x)$ on $[0,1]$ with Matérn covariance and two faint spectral lines sampled on a jittered grid. Individual estimates for many samples are shown in grey, their average is shown in black, and the true SDF $S$ is shown in red. The forward estimator is as defined in (\ref{['eq:shat']}) and $g(x)$ is a standard Kaiser window.
• Figure 2: Accuracy of various weighting schemes in recovering $G$ using a weighted Fourier sum with $\left\{ x_j\right\}_{j=1}^n \mathrel{\overset{\hbox{i.i.d.}}{\hbox{[}1]{$∼$}}} \text{Unif}([-1, 1])$. Weights for our method are computed with $\Omega = 75$. The window $g$ is chosen to be a Kaiser function, discussed in the next subsection, and the weights are $\alpha_j = g(x_j)/n$ for no correction and $\alpha_j = g(x_j) \gamma_j$ as in \ref{['eq:irtrap']} for the irregular trapezoidal scheme.
• Figure 3: A comparison of the recovered window $\left|H_{\alpha}\right|^2$ in a sampling setting with $\left\{ x_j\right\}_{j=1}^n \mathrel{\overset{\hbox{i.i.d.}}{\hbox{[}1]{$∼$}}} \text{Unif}([-1, 1]\setminus[-0.15,0])$ and $n = 3000$ using a generic Kaiser window supported on $[-1, 1]$ and a prolate window supported on the union of disjoint intervals. Blue dots in the left subplot show a representative subset of sampling locations.
• Figure 4: A visual demonstration of the phenomenon described in Theorem \ref{['thm:aliasing']}. Left panel: the expected value of the estimator $\hat{S}(\omega)$ at various frequencies compared with the true $S(\omega)$ and the aliasing bias $\varepsilon(\omega)$. Right panel: $|H_{\alpha}(\omega - \xi)|^2$ for two values of $n$ but fixed $\Omega$ visualized against $S(\omega)$, where the integral of the product of those two functions is $\mathbb{E} \hat{S}(\omega)$. The two faint dotted lines give $\left\Vert\bm{\alpha}\right\Vert_{2}^2$ for each data size.
• Figure 5: Runtime measurements for weight computations in two sampling regimes: uniformly random points and gappy gridded points. For random points, a hierarchical matrix preconditioner is used to accelerate the Krylov solver. In the gappy gridded case, no preconditioner is used.

Theorems & Definitions (10)

• Theorem 1
• Corollary 1
• proof
• Theorem 2
• Corollary 2
• proof
• Lemma 1
• proof
• proof : Proof of Theorem \ref{['thm:aliasing']}
• proof