cofi.utils.SPDEMaternReg

cofi.utils.SPDEMaternReg#

class cofi.utils.SPDEMaternReg(model_shape: tuple, ell: float, sigma: float = 1.0, grid_spacing=1.0, reference_model: ndarray | None = None)[source]#

Sparse Matérn ν=1 regularization for 2D spatial fields via the SPDE approach.

Implements the regularization term

\[\mathcal{R}(\mathbf{m}) = \|\mathbf{R}(\mathbf{m} - \mathbf{m}_0)\|^2, \quad \mathbf{R} = \tau \mathbf{M}^{1/2}(\kappa^2 \mathbf{I} - \mathbf{L}_h), \quad \kappa = \sqrt{2} / \ell, \quad \tau = \frac{1}{2 \sqrt{\pi}\,\kappa\,\sigma}\]

where \(\mathbf{L}_h\) is the 2D grid-spacing-aware Neumann Laplacian, \(\mathbf{M}\) is the lumped mass matrix, \(\ell\) is the Matérn length scale for \(\nu = 1\), and \(\sigma\) is the marginal standard deviation in the continuous SPDE model. For the uniform rectangular grids implemented here, \(\mathbf{M} \approx h_x h_y \mathbf{I}\) so \(\mathbf{R} = \tau \sqrt{h_x h_y}(\kappa^2 \mathbf{I} - \mathbf{L}_h)\). The matrix \(\mathbf{R}\) is a sparse precision factor of the Matérn ν=1 precision matrix

\[\mathbf{Q} = \mathbf{R}^\top \mathbf{R} = \tau^2 \mathbf{B}_h^\top \mathbf{M} \mathbf{B}_h, \quad \mathbf{B}_h = \kappa^2 \mathbf{I} - \mathbf{L}_h\]

which reduces to \(\tau^2 (\kappa^2 \mathbf{I} - \mathbf{L})^2\) on a unit grid. (Lindgren, Rue & Lindström 2011). The matrix \(\mathbf{R}\) is sparse (\(O(n)\) non-zeros), unlike the dense precision matrix of a standard Gaussian prior.

The Matérn practical range (distance at which correlation drops to ~0.14) is \(\rho = 2\ell\) for \(\nu = 1\).

Parameters:

  • model_shape (tuple of (int, int)) – Shape of the 2D model grid (n_lon, n_lat) (or equivalently any two positive integers whose product is the number of model parameters).

  • ell (float) – Matérn length scale in the same physical units as grid_spacing. For \(\nu=1\), this is related to the SPDE wavenumber by \(\kappa = \sqrt{2}/\ell\). The practical range (correlation ≈ 0.14) is \(\rho = 2\ell\).

  • sigma (float, optional) – Prior marginal standard deviation of model perturbations. The internal SPDE scaling uses \(\tau = 1 / (2 \sqrt{\pi}\,\kappa\,\sigma)\). Default 1.0.

  • grid_spacing (float or tuple of (float, float), optional) – Physical spacing of the regular grid in the longitude/x and latitude/y directions. A scalar applies the same spacing to both axes. grid_spacing=1.0 reproduces the current unit-grid interpretation.

  • reference_model (np.ndarray, optional) – Background model \(\mathbf{m}_0\). If provided, the penalty is on \(\mathbf{m} - \mathbf{m}_0\); if omitted, the penalty is on \(\mathbf{m}\) directly.

Notes

Boundary conditions

The discrete Laplacian uses Neumann (zero-flux) boundary conditions. This means no smoothness penalty is imposed across the grid edge, so boundary cells are treated identically to interior cells. Dirichlet (zero-endpoint) boundaries would artificially anchor boundary cells to the reference model and are not appropriate for most geophysical grids. As in the SPDE literature, these boundary conditions inflate marginal variances near the edges of a finite grid.

Examples

>>> from cofi.utils import SPDEMaternReg
>>> import numpy as np
>>> reg = SPDEMaternReg(model_shape=(10, 8), ell=2.0, sigma=0.02)
>>> reg(np.zeros((10, 8)))
0.0

Reference Details

__call__(model: ndarray) → Number#

a class instance itself can also be called as a function

It works exactly the same as reg.

In other words, the following two usages are exactly the same:

>>> my_reg = QuadraticReg(factor=1, model_size=3)
>>> my_reg_value = my_reg(np.array([1,2,3]))            # usage 1
>>> my_reg_value = my_reg.reg(np.array([1,2,3]))        # usage 2
gradient(model: ndarray) → ndarray#

the gradient of regularization function with respect to model given a model

The usual size for the gradient is \((M,)\) where \(M\) is the number of model parameters

hessian(model: ndarray) → ndarray#

the hessian of regularization function with respect to model given a model

The usual size for the Hessian is \((M,M)\) where \(M\) is the number of model parameters

reg(model: ndarray) → Number#

the regularization function value given a model to evaluate

ell#

Matérn length scale.

grid_shape#

Original 2D grid shape (n_lon, n_lat).

grid_spacing#

Physical grid spacing (h_lon, h_lat).

kappa#

Wavenumber parameter κ = sqrt(2) / ell for ν = 1.

matrix#

the regularization matrix

This is either an identity matrix, or first/second order difference matrix (generated by Python package findiff), or a custom matrix brought on your own.

model_shape#
model_size#

the number of unknowns that current regularization function accepts

rho#

Matérn practical range (rho = 2 * ell for nu=1).

sigma#

Prior marginal standard deviation.

back to top

Contents