First example with a diagonal covariance matrix

In the following, we define a (3 x 3) covariance matrix defining only its diagonal (i.e., all elements of the random vectors are independent). Using scipy, it is possible to compute the pdf.

d = [1, 2, 3]
A33 = np.diag(d)  # a diagonal covariance matrix
x = [4, -2, 5]  # a point of interest
dist = sp.stats.multivariate_normal(mean=[0, 0, 0], cov=A33)
dist.pdf(x)

np.float64(4.9595685102808205e-08)

It is possible to obtain a dense representation in a straightforward manner:

cov_diag33 = covmats.CovViaDiagonal(d)
dist = sp.stats.multivariate_normal(mean=[0, 0, 0], cov=cov_diag33)
dist.pdf(x)

np.float64(4.9595685102808205e-08)

It is compatible with the stats API from scipy since the base class inherit from Covariance.

Every representation also exposes todense, solve, precision and log_pdet, as well as sample_mvnormal, whiten and colorize for fast Monte-Carlo sampling and data whitening:

samples = cov_diag33.sample_mvnormal(shape=[10000], random_state=rng_seed)
np.round(np.cov(samples, rowvar=False), 1)

array([[ 1.,  0., -0.],
       [ 0.,  2., -0.],
       [-0., -0.,  3.]])

Cholesky and precision-based representations

CovViaCholesky wraps a dense Cholesky factor L such that A = L @ L.T, while CovViaPrecisionCholesky does the same for the precision matrix Q = A^-1 (the natural representation for Gaussian Markov Random Fields, where Q is sparse). Sparse counterparts, CovViaSparseCholesky and CovViaSparsePrecisionCholesky, build on a SparseCholeskyFactor (an LDL’ factorization) for large, sparse problems:

rng = np.random.default_rng(rng_seed)
n = 4
B = rng.random((n, n))
A = B @ B.T + n * np.eye(n)  # a random SPD matrix
cov_cho = covmats.CovViaCholesky(np.linalg.cholesky(A))
np.allclose(cov_cho.todense(), A)

True

Low-rank representations

For very large problems, CovViaEigenFactorization (a truncated eigen decomposition) and CovViaEnsemble (an ensemble of anomalies, as used in ensemble Kalman filtering) never require forming the dense covariance matrix:

ensemble = rng.multivariate_normal(np.zeros(6), np.eye(6) * 2 + 0.3, size=200)
cov_ens = covmats.CovViaEnsemble(ensemble)
cov_ens.shape

(6, 6)

get_linop_eigen_factorization and eigen_factorize_cov_mat build a randomized low-rank eigen factorization directly from any CovarianceMatrix or CovKernelAsLinop instance, and get_explained_var reports how much variance each retained mode captures – useful to pick a truncation rank.

Kernel-based linear operators

When the covariance is defined analytically through a kernel evaluated on a point cloud, CovKernelAsLinop (dense evaluation, any point cloud) and CovKernelAsLinopViaFFT (FFT-based, regular grids only) expose a LinearOperator without requiring the covariance matrix to be assembled up front:

pts = covmats.get_pts_coords_regular_grid(mesh_dim=1.0, shape=(6, 6))
cov_kernel = covmats.CovKernelAsLinop(
    pts, lambda d: np.exp(-d), len_scale=np.array([2.0, 2.0])
)
cov_kernel.shape

(36, 36)

Priors and drift matrices

Alongside covariance representations, covmats provides PriorTerm subclasses to describe the deterministic mean/trend of a field, from a simple NullPriorTerm to a LinearDriftMatrix expressing a spatially-varying trend as m = X @ beta:

drift = covmats.LinearDriftMatrix(pts)
drift.beta = np.array([1.0, 0.5, -0.5])
drift.get_values(np.zeros(pts.shape[0])).shape

(36,)

Conditional and unconditional simulations

See https://github.com/antoinecollet5/covmats/blob/main/notebooks/examples_simulations.ipynb

🏗️ Complete example with supporting paper coming Q1 2027.