Correlation modules refactor #12 – circulant embedding core

[raoanirudh]

Efficient approach for simulating spatially-cross-correlated ground motion fields

Following 11 PRs linked to #11230 and focused on establishing the correlation architecture, documentation, model implementations and verifications, this PR begins the performance phase of the correlation module improvements.

This PR is linked with #10833 and #11230 (comment).

It introduces a model-independent multivariate circulant-embedding core for efficiently simulating spatially and cross-IM-correlated residual fields on regular grids. It does not alter any calculator path yet and therefore no test numbers are changed. This PR establishes and verifies the numerical foundation before it is connected to scenario and event-based calculations in following PRs.

Why performance of spatial (and cross-IMT) correlation is crucial

Spatial correlation determines whether nearby sites tend to experience high or low within-event residuals together. Ignoring it can substantially distort the distribution of geographically aggregated damage and loss.

Cross-IM correlation is equally important when fragility, damage, or consequence calculations depend on several intensity measures. Simulating GMFs for every IMT independently can produce physically inconsistent combinations at a site and incoherent spatial patterns between IMTs.

Joint spatial–cross-IM simulation preserves both relationships simultaneously. Several studies have established the importance of including spatial-cross-correlation for risk assessment, Weatherill et al. (2015) being one such example.

While correlation can be switched off for unconditioned GMF simulations, calculations conditioned on seismic-station observations require correlation models as part of the conditioning algorithm.

Realistic target calculations include:

• scenario calculations involving anywhere from ~10,000 sites (eg. 2014 M6.0 South Napa earthquake) to up to ~250,000 sites (eg. 2023 M7.8 Kahramanmaraş earthquake), using 4–5 IMTs and 500–1,000 ground-motion field realizations per GMM
• event-based calculations requiring GMF generation for hundreds of thousands of events
• near-real-time scenario applications following significant earthquakes, potentially conditioned on available station data, where computational time is of crucial importance

The existing method in the engine constructs and factorizes a covariance matrix covering every site–IMT combination. For $$N$$ sites and $$M$$ IMTs, this matrix would have $$(MN)^2$$ elements. Its memory demand therefore grows quadratically with the number of sites, while its factorization cost grows even faster.

For the two example performance cases mentioned above, we have the following memory requirements with the existing method:

• the 2014 South Napa 500-field conditioned calculation is estimated to require 40 GB – still manageable with a 64 GB memory allocation
• the 2023 Kahramanmaraş conditioned calculation is estimated to require 19.4 TB!

In practice, these computational limits have often forced users to disable spatial correlation altogether, because realistically sized correlated calculations simply cannot run to completion with the existing approach, eg. see #11695. Without performance improvements and consideration of more efficient algorithms, these calculations are essentially infeasible on typical workstations.

Circulant embedding for regular grids

We begin by looking at the easier case of GMFs simulated on regularly spaced grids. A regular spatial grid contains a lot of repeated covariance information. Pairs of points separated by the same distance (or grid offset) have the same covariance. In the existing approach in the engine, we build dense matrices to store and process all of those repeated entries individually.

Circulant embedding (CE) places the required grid inside a somewhat larger periodic grid. The resulting covariance has a repeating structure that can be transformed into the frequency domain using fast Fourier transforms (FFTs). https://www.joshuahtouyz.com/blog/251210_fast_multivariate_grf_simulation provides a useful visual explanation.

In simplified terms, the CE method:

1. embeds the requested grid inside a larger periodic grid
2. represents the repeated covariance structure in the frequency domain
3. factorizes only a small $$M \times M$$ cross-IMT covariance matrix at each frequency
4. uses FFTs to generate complete correlated fields
5. extracts the requested sites from the larger grid

This avoids constructing the enormous dense site-by-site covariance matrix as is being done with the existing approach in the engine. The memory requirement grows only (approximately) linearly with the number of embedded grid cells rather than its square, and each additional realization can be generated using FFT operations.

Circulant embedding seems to be the most promising approach for our expected combination of large grids, a modest number of IMTs, and many realizations. It should also allow the expensive spectral factorization to be reused across batches of realizations and, where applicable, across events sharing the same grid and correlation setup.

Prior USGS work

Verros et al. (2017)¹ developed a memory-efficient iterative parallel algorithm based on the decay of the spatial correlation function. Their work demonstrated both the prohibitive cost of conventional large-grid simulation. Bailey et al. (2022)² subsequently adapted circulant embedding to conditional simulation with irregularly located observations. Because standard circulant embedding operates on regular grids, they investigated approximate extensions using local kriging and nearest-neighbor kriging. Their work was explicitly motivated by near-real-time ShakeMap applications and demonstrated that conditioned ensembles could be generated on an interactive timescale.

Pilot implementation and verification in this PR

The pilot implementation in this PR follows the circulant-embedding formulation of Dietrich and Newsam (1993)³ and its multivariate extension by Chan and Wood (1999)⁴, adding only a regular-grid unconditional core algorithm. The Bailey et al. (2022) approach would be an important reference when we add support for stations that do not coincide with target-grid cells.

This pilot method evaluates the spatial and cross-IM covariance over a rectangular periodic embedding, factorizes the small cross-IM spectral covariance matrix at each Fourier mode, and applies the resulting factors to independent residual samples using real FFTs.

The implementation supports rectangular grids and filtered subsets of grid cells while preserving the requested site and IMT ordering, selects FFT-efficient dimensions, and automatically enlarges an embedding when its initial spectrum is not positive semidefinite.

The verification tests apply the FFT factor to an identity basis and reconstruct the complete covariance matrix, avoiding the uncertainty of Monte Carlo testing.

The end goal is scalable joint spatial–cross-IM simulation for both routine event-based risk calculations and near-real-time earthquake scenarios, whether or not station observations are available for the scenarios. Conditioned GMF scenarios should be runnable on a workstation or a computer with a reasonable amount of available memory, say 36 GB.

References

Footnotes

1. Verros, S. A., Wald, D. J., Worden, C. B., Hearne, M. G., & Ganesh, M. (2017). Computing spatial correlation of ground motion intensities for ShakeMap. Computers & Geosciences, 99, 145–154. https://doi.org/10.1016/j.cageo.2016.11.004 ↩

2. Bailey, M. D., Bandyopadhyay, S., & Nychka, D. W. (2022). Adapting conditional simulation using circulant embedding for irregularly spaced spatial data. Stat, 11(1), e446. https://doi.org/10.1002/sta4.446 ↩

3. Dietrich, C. R., & Newsam, G. N. (1993). A fast and exact method for multidimensional Gaussian stochastic simulations. Water Resources Research, 29(8), 2861–2869. https://doi.org/10.1029/93WR01070 ↩

4. Chan, G., & Wood, A. T. A. (1999). Simulation of stationary Gaussian vector fields. Statistics and Computing, 9, 265–268. https://doi.org/10.1023/A:1008903804954 ↩

[micheles]

Looks fine to me, but I would be curious to see some performance numbers comparing this implementation with what we have now on a realistic case.

[raoanirudh]

Looks fine to me, but I would be curious to see some performance numbers comparing this implementation with what we have now on a realistic case.

Indeed. For that, we'll have to plug this in to the gmf computer, which needs to be done carefully, and it needs to be an option the user can select. Will do that in the next pull request(s).