Modeling Large Nonstationary Spatial Data with the Full-Scale Basis Graphical Lasso

Technical Summary: Modeling Large Nonstationary Spatial Data with the Full-Scale Basis Graphical Lasso

Problem Statement
Spatially distributed datasets have grown increasingly large and complex, necessitating statistical models that can handle nonstationary processes while remaining computationally feasible. Traditional spatial statistics often rely on sparsity or rank-dependent structures in the covariance matrix. While basis function expansions (e.g., Fixed Rank Kriging, LatticeKrig, Basis Graphical Lasso) offer interpretability and flexibility, they face specific limitations:

• Fixed Rank Kriging often fails to capture small-scale dependence.
• LatticeKrig enforces approximate stationarity through spatial autoregressive structures on coefficients, limiting its ability to model complex nonstationary behaviors.
• Basis Graphical Lasso (BGL) relaxes stationarity assumptions by estimating a sparse precision matrix for coefficients but typically assumes the residual process is a simple nugget effect (white noise), failing to capture small-scale spatial correlation.
• Full-Scale Approximations (FSA) combine low-rank and small-scale components but often rely on empirical orthogonal functions (EOFs) for the low-rank portion or require estimating dense covariance matrices, which becomes computationally prohibitive for large numbers of basis functions.

The core challenge is to synthesize these approaches into a single framework that captures both large-scale nonstationary structures and small-scale locally correlated behavior without incurring prohibitive computational costs.

Methodology
The authors propose the Full-Scale Basis Graphical Lasso (FSBGL), a model that decomposes a spatial process $Z(\vec{x})$ into a low-rank component $Z_1(\vec{x})$ and a small-scale component $Z_2(\vec{x})$:
$$Z(\vec{x}) = Z_1(\vec{x}) + Z_2(\vec{x})$$

1. Low-Rank Component ($Z_1$): Modeled as a basis expansion $Z_1(\vec{x}) = \sum \beta_j \psi_j(\vec{x})$. Unlike previous methods that assume independence or stationarity among coefficients, the vector of coefficients $\vec{\beta}$ is endowed with a sparse precision matrix $Q$ (inverse covariance). The sparsity pattern of $Q$ is estimated using a Graphical Lasso penalty, allowing for flexible, nonstationary conditional dependencies between basis coefficients.
2. Small-Scale Component ($Z_2$): Modeled as a stationary, isotropic Gaussian process with a compactly supported covariance function $C_2$. This ensures the covariance matrix $D$ of the residuals is sparse, distinguishing FSBGL from BGL (which sets $Z_2=0$) and allowing for the modeling of small-scale spatial correlation.
3. Estimation: The estimation employs a penalized likelihood approach. The authors derive a regularized likelihood function (Eq. 9) that incorporates the sparsity of $D$ and the graphical lasso penalty on $Q$.
   • The optimization for $Q$ is solved using a Difference-of-Convex (DC) scheme, leveraging the QUIC algorithm for efficiency.
   • The small-scale parameters (defining $D$) and the low-rank precision parameters are estimated in a two-step or joint optimization process, utilizing the sparsity of $D$ to reduce computational complexity from $O(n^3)$ to more manageable levels (e.g., $O(k^2n)$ for banded matrices).
4. Model Selection: The optimal penalty parameter $\lambda$ is selected using the Conditional Akaike Information Criterion (cAIC), which accounts for the effective degrees of freedom in the model.

Key Contributions

• Synthesis of Methods: The FSBGL unifies the Full-Scale Approximation (low-rank + small-scale) and the Basis Graphical Lasso (sparse precision for coefficients) into a single framework.
• Flexible Nonstationarity: By applying the graphical lasso to the basis coefficients, the model captures complex, nonstationary dependencies across different scales of resolution, a feature not possible in LatticeKrig where resolution levels are often assumed independent.
• Small-Scale Modeling: Unlike the BGL, FSBGL explicitly models small-scale spatial correlation via a compactly supported covariance function ($Z_2$), rather than relying solely on a nugget effect.
• Computational Efficiency: The method leverages the sparsity of both the precision matrix $Q$ and the small-scale covariance $D$ to make estimation feasible for large datasets (tens of thousands of locations) and large numbers of basis functions.
• Basis Agnosticism: The framework allows for the use of various basis functions (e.g., Needlets, Radial Basis Functions) chosen by the user based on domain expertise, rather than being restricted to specific fixed bases.

Results
The paper evaluates the FSBGL through simulation studies and an application to thermospheric neutral temperature fields derived from the WACCM-X Earth system model.

• Simulation Studies: Using synthetic fields with known precision matrix structures (block-diagonal, hub-and-spoke, random graphs), the FSBGL successfully recovered the structure of $Q$ and the parameters of the residual process. The accuracy of recovering $Q$ improved with the number of independent realizations, and the method performed comparably to the BGL in recovering known structures.
• Thermospheric Application:
  • Data: The model was applied to high-resolution (downsampled to $2^\circ \times 2^\circ$) neutral temperature data (41 replicates) representing the lower thermosphere (80–150 km).
  • Model Fit: Among various small-scale covariance options (Gaspari-Cohn, Wendland, Tapered Matérn), the Tapered Matérn model yielded the best fit (lowest cAIC) and the most physically interpretable parameters (e.g., a smoothness parameter $\nu \approx 0.63$ indicating a non-differentiable field, consistent with atmospheric turbulence).
  • Structure Recovery: The estimated precision matrix $Q$ revealed significant conditional dependencies between the lowest resolution level (global scale) and higher resolution levels, a feature LatticeKrig cannot capture. The neighborhood structures showed physical interpretations, such as solar heating effects and latitude-dependent dynamics.
  • Comparative Performance:
    • Likelihood: FSBGL achieved the lowest negative log-likelihood on out-of-sample data compared to LatticeKrig, BGL, and standard Full-Scale Approximation (FSA).
    • Prediction: In subgrid-scale predictions (predicting at $0.25^\circ$ resolution from $2^\circ$ training data), FSBGL significantly outperformed competitors. The Mean and Median Continuous Rank Probability Scores (CRPS) and Root Mean Square Error (RMSE) for FSBGL were less than half those of the BGL and LatticeKrig.
    • Generalization: Unlike the standard FSA (which relies on EOFs and tends to overfit training data), FSBGL demonstrated superior generalization to held-out data.

Significance and Claims
The paper claims that the FSBGL provides a flexible and computationally tractable approach for modeling large, nonstationary spatial datasets. Its primary significance lies in its ability to:

1. Capture Salient Features: It better captures the complex, multi-scale, and nonstationary features of thermospheric temperature fields than state-of-the-art alternatives.
2. Enable Subgrid Prediction: The model's ability to predict subgrid-scale variations suggests its potential utility as a statistical downscaler. This could allow engineers to run computationally expensive climate models at lower resolutions and use FSBGL to statistically estimate fine-scale structures and uncertainties, which is critical for modeling aerodynamic forces on reentry vehicles.
3. Domain Adaptability: The framework is adaptable to various applications by allowing users to select appropriate basis functions and residual covariance structures based on domain knowledge.

The authors conclude that while the method shows promise for operational use in augmenting ensemble simulations for space vehicle reentry, further investigation is needed to confirm its operational accuracy and efficiency in real-time scenarios.