The LU decomposition approach to simulation

§7.1 motivated conditional simulation. §7.2 develops the mathematical FOUNDATION: how do you actually generate a random field with a specified covariance structure? The CHOLESKY (LU) decomposition gives the exact answer for modest grid sizes.

The algorithm

To generate a random vector $y \in \mathbb{R}^N$ with covariance matrix C:

• Build the covariance matrix C from the variogram model. C[i, j] = sill × exp(-3 d_ij / range).
• Cholesky-factor C = LL' where L is lower-triangular.
• Generate independent standard-Normal noise $z = (z_1, \ldots, z_N) \sim N(0, I)$.
• Transform: $y = Lz$.

The result has $\text{Cov}(y) = E[(Lz)(Lz)'] = L E[zz'] L' = L I L' = LL' = C$. Exact covariance match.

Cholesky decomposition

For a positive-definite matrix C, the Cholesky algorithm constructs lower-triangular L row by row:

Time: O(N³). Memory: O(N²) for C + O(N²) for L (total 2N²).

The scalability problem

For N = 1000 grid points: memory ≈ 8 MB (two N² doubles); compute time ≈ 1-2 seconds. Feasible.

For N = 100,000 grid points: memory ≈ 80 GB; compute time ≈ days. Infeasible on standard hardware.

So Cholesky/LU simulation is restricted to MODEST grids. Real geostatistical applications routinely involve grids with 10^4 to 10^6 nodes. The Cholesky approach doesn't scale — hence §7.3 SEQUENTIAL Gaussian simulation, which trades exactness for scalability via sequential node-by-node sampling using a moving search neighbourhood.

When Cholesky/LU IS used

• Educational implementations and small benchmark problems.
• Validating SGS — check that SGS realisations match the variogram structure that Cholesky simulations would.
• Small dense grids (≤ 1000 nodes) where exact simulation is feasible.
• Spectral methods variants for very specialised regular grids.

For most modern geostat applications, SGS is the workhorse. But Cholesky/LU is the conceptual foundation.

Try it

• Default N = 30, range = 2.0. The widget computes Cholesky of the 30×30 covariance matrix, transforms 5 independent z-vectors, displays 5 realisations.
• The empirical lag-1 correlation across realisations matches the expected correlation from the variogram. Verifies Cholesky correctness.
• Drag range to 0.5 (short correlation). Realisations look noisy — each grid point is nearly independent of its neighbours. Empirical correlation is small.
• Drag range to 5.0 (long correlation). Realisations are smooth, slowly-varying functions. Adjacent points are highly correlated.
• Crank N to 60. Computation is still fast (60×60 Cholesky is trivial). Realisations look similar in structure but at finer resolution. Scaling to N=1000 would be ~30 seconds; N=10,000 would be ~hours; N=100,000 would be days — hence the need for SGS.

A geostatistician needs to simulate permeability on a 200 × 200 × 50 reservoir grid (2 million nodes). Why is Cholesky/LU not feasible, and what alternative method should they use?

What you now know

Cholesky/LU decomposition is the mathematical foundation of simulation: factor C = LL', transform y = Lz where z is iid standard Normal. Exact covariance match. O(N²) memory + O(N³) time — restricts to N ≤ ~1000. For larger grids, use SGS (§7.3). The Cholesky approach remains pedagogically important and benchmark-relevant.

Mark section complete →

References

• Davis, M.W. (1987). "Production of conditional simulations via the LU triangular decomposition of the covariance matrix." Mathematical Geology 19(2), 91–98. (Davis original.)
• Alabert, F. (1987). "The practice of fast conditional simulations through the LU decomposition of the covariance matrix." Mathematical Geology 19(5), 369–386.
• Goovaerts, P. (1997). Geostatistics for Natural Resources Evaluation, Chapter 8. Oxford.
• Chilès, J.-P., Delfiner, P. (2012). Geostatistics: Modeling Spatial Uncertainty, 2nd ed. Wiley. Section 7.3.
• Cressie, N. (1993). Statistics for Spatial Data. Wiley. (Comprehensive treatment of simulation methods.)