GMRFs.jl

Appearance

GMRFs ​

GaussianMarkovRandomFields.AbstractGMRF Type
julia
AbstractGMRF

A Gaussian Markov Random Field (GMRF) is a special case of a multivariate normal distribution where the precision matrix is sparse. The zero entries in the precision correspond to conditional independencies.

source
GaussianMarkovRandomFields.GMRF Type
julia
GMRF(mean, precision, alg=LinearSolve.DefaultLinearSolver(); Q_sqrt=nothing, rbmc_strategy=RBMCStrategy(1000), linsolve_cache=nothing)

A Gaussian Markov Random Field with mean mean and precision matrix precision.

Arguments

  • mean::AbstractVector: The mean vector of the GMRF.

  • precision::Union{LinearMap, AbstractMatrix}: The precision matrix (inverse covariance) of the GMRF.

  • alg: LinearSolve algorithm to use for linear system solving. Defaults to LinearSolve.DefaultLinearSolver().

  • Q_sqrt::Union{Nothing, AbstractMatrix}: Square root of precision matrix Q, used for sampling when algorithm doesn't support backward solve.

  • rbmc_strategy: RBMC strategy for marginal variance computation when selected inversion is unavailable. Defaults to RBMCStrategy(1000).

  • linsolve_cache::Union{Nothing, LinearSolve.LinearCache}: Existing LinearSolve cache to reuse. If nothing, creates a new cache. Useful for iterative algorithms requiring factorization reuse.

Type Parameters

  • T<:Real: The numeric type (e.g., Float64).

  • PrecisionMap<:Union{LinearMap{T}, AbstractMatrix{T}}: The type of the precision matrix.

Fields

  • mean::Vector{T}: The mean vector.

  • precision::PrecisionMap: The precision matrix.

  • Q_sqrt::Union{Nothing, AbstractMatrix{T}}: Square root of precision matrix for sampling.

  • linsolve_cache::LinearSolve.LinearCache: The LinearSolve cache for efficient operations.

  • rbmc_strategy: RBMC strategy for variance computation fallback.

Notes

The LinearSolve cache is constructed automatically (if not provided) and is used to compute means, variances, samples, and other GMRF quantities efficiently. The algorithm choice determines which optimization strategies (selected inversion, backward solve) are available. When selected inversion is not supported, marginal variances are computed using the configured RBMC strategy. Providing an existing linsolve_cache enables factorization reuse in iterative algorithms.

source
GaussianMarkovRandomFields.ChordalGMRF Type
julia
ChordalGMRF{T, Hrm, Fac, Mea} <: AbstractGMRF{T, Hrm}

A GMRF backed by a chordal Cholesky factorization (via CliqueTrees.Multifrontal.ChordalCholesky) instead of CHOLMOD.

The pure-Julia chordal factorization composes naturally with Mooncake's reverse-mode AD through the rrules shipped by MooncakeSparse, so logpdf and gaussian_approximation give correct gradients with respect to the hyperparameters that produced Q. This is the recommended GMRF type for Mooncake-based hyperparameter optimization (e.g. L-BFGS / Adam on the marginal likelihood).

Fields

  • μ::AbstractVector: Mean.

  • Q::Hermitian: Precision matrix.

  • F::ChordalCholesky: Chordal Cholesky factorization of Q.

Construction

julia
ChordalGMRF(μ, Q)              # factorize Q via ChordalCholesky
ChordalGMRF(μ, Q, F)           # reuse a precomputed factorization
source
GaussianMarkovRandomFields.precision_map Function
julia
precision_map(::AbstractGMRF)

Return the precision (inverse covariance) map of the GMRF.

source
julia
precision_map(d::ConstrainedGMRF)

Return the precision map of the unconstrained base GMRF.

The precision of the constrained distribution x | Ax = e is rank-deficient (its null space is spanned by the columns of A'), so there is no canonical non-singular precision matrix to return. The base precision is what the unnormalised log-density on the constraint manifold uses,

log p(x | Ax = e) ∝ -0.5 (x - μ)' Q (x - μ)   for x with Ax = e,

so generic methods like sqmahal and gradlogpdf give the correct values when evaluated at points that satisfy the constraint.

source
GaussianMarkovRandomFields.precision_matrix Function
julia
precision_matrix(::AbstractGMRF)

Return the precision (inverse covariance) matrix of the GMRF.

source
julia
precision_matrix(model::LatentModel; θ...) -> AbstractMatrix

Precision matrix of the Gaussian latent prior at hyperparameters θ.

source
GaussianMarkovRandomFields.InformationVector Type
julia
InformationVector(data::AbstractVector)

Wrapper type for information vectors (Q * μ) used in GMRF construction. This allows distinguishing between constructors that take mean vectors vs information vectors.

source
GaussianMarkovRandomFields.information_vector Function
julia
information_vector(d::GMRF)

Return the information vector (Q * μ) for the GMRF. If stored, returns the cached value; otherwise computes it.

source

Metadata ​

GaussianMarkovRandomFields.MetaGMRF Type
julia
MetaGMRF{M <: GMRFMetadata, T, P, G <: AbstractGMRF{T, P}} <: AbstractGMRF{T, P}

A wrapper that combines a core GMRF with metadata of type M. This allows for specialized behavior based on the metadata type while preserving the computational efficiency of the underlying GMRF.

Fields

  • gmrf::G: The core computational GMRF (parametric type)

  • metadata::M: Domain-specific metadata

Usage

julia
# Define metadata types
struct SpatialMetadata <: GMRFMetadata
    coordinates::Matrix{Float64}
    boundary_info::Vector{Int}
end

# Create wrapped GMRF
meta_gmrf = MetaGMRF(my_gmrf, SpatialMetadata(coords, boundary))

# Dispatch on metadata type for specialized behavior
function some_spatial_operation(mgmrf::MetaGMRF{SpatialMetadata})
    # Access coordinates via mgmrf.metadata.coordinates
    # Access GMRF via mgmrf.gmrf
end
source
GaussianMarkovRandomFields.GMRFMetadata Type
julia
GMRFMetadata

Abstract base type for metadata that can be attached to GMRFs via MetaGMRF. Concrete subtypes should contain domain-specific information about the GMRF structure, coordinates, naming, etc.

source

Arithmetic ​

GaussianMarkovRandomFields.linear_condition Function
julia
linear_condition(gmrf::GMRF; A, Q_ϵ, y, b=zeros(size(A, 1)), obs_precision_contrib=nothing)

Condition a GMRF on linear observations y = A * x + b + ϵ where ϵ ~ N(0, Q_ϵ^(-1)).

Arguments

  • gmrf::GMRF: The prior GMRF

  • A::Union{AbstractMatrix, LinearMap}: Observation matrix

  • Q_ϵ::Union{AbstractMatrix, LinearMap}: Precision matrix of observation noise

  • y::AbstractVector: Observation values

  • b::AbstractVector: Offset vector (defaults to zeros)

  • obs_precision_contrib: Precomputed A' * Q_ϵ * A (optional optimization)

Returns

A new GMRF representing the posterior distribution with updated mean and precision.

Performance Optimization

When A has special structure (e.g., index selection), users can precompute A' * Q_ϵ * A more efficiently and pass it to avoid redundant computation.

Notes

Uses information vector arithmetic for efficient conditioning without intermediate solves.

source
GaussianMarkovRandomFields.joint_gmrf Function
julia
joint_gmrf(
    x1::AbstractGMRF,
    A::AbstractMatrix,
    Q_ϵ::AbstractMatrix,
    b::AbstractVector=spzeros(size(A, 1))
)

Return the joint GMRF of x1 and x2 = A * x1 + b + ϵ where ϵ ~ N(0, Q_ϵ⁻¹).

Arguments

  • x1::AbstractGMRF: The first GMRF.

  • A::AbstractMatrix: The matrix A.

  • Q_ϵ::AbstractMatrix: The precision matrix of the noise term ϵ.

  • b::AbstractVector=spzeros(size(A, 1)): Offset vector b; optional.

Returns

A GMRF object representing the joint GMRF of x1 and x2 = A * x1 + b + ϵ.

source

Spatiotemporal setting ​

Types ​

GaussianMarkovRandomFields.AbstractSpatiotemporalGMRF Type
julia
AbstractSpatiotemporalGMRF

A spatiotemporal GMRF is a GMRF that explicitly encodes the spatial and temporal structure of the underlying random field. All time points are modelled in one joint GMRF. It provides utilities to get statistics, draw samples and get the spatial discretization at a given time.

source
GaussianMarkovRandomFields.ImplicitEulerConstantMeshSTGMRF Type
julia
ImplicitEulerConstantMeshSTGMRF{D, SSM, T, P, G}

Implicit Euler spatiotemporal GMRF with constant spatial discretization, as a MetaGMRF over ImplicitEulerMetadata.

source
GaussianMarkovRandomFields.ConcreteConstantMeshSTGMRF Type
julia
ConcreteConstantMeshSTGMRF{D, T, P, G}

Concrete spatiotemporal GMRF with constant spatial discretization, as a MetaGMRF over ConcreteSTMetadata.

source

Quantities ​

GaussianMarkovRandomFields.N_t Function
julia
N_t(::AbstractSpatiotemporalGMRF)

Return the number of time points in the spatiotemporal GMRF.

source
GaussianMarkovRandomFields.time_means Function
julia
time_means(::AbstractSpatiotemporalGMRF)

Return the means of the spatiotemporal GMRF at each time point.

Returns

  • A vector of means of length Nₜ, one for each time point.
source
GaussianMarkovRandomFields.time_vars Function
julia
time_vars(::AbstractSpatiotemporalGMRF)

Return the marginal variances of the spatiotemporal GMRF at each time point.

Returns

  • A vector of marginal variances of length Nₜ, one for each time point.
source
GaussianMarkovRandomFields.time_stds Function
julia
time_stds(::AbstractSpatiotemporalGMRF)

Return the marginal standard deviations of the spatiotemporal GMRF at each time point.

Returns

  • A vector of marginal standard deviations of length Nₜ, one for each time point.
source
GaussianMarkovRandomFields.time_rands Function
julia
time_rands(::AbstractSpatiotemporalGMRF, rng::AbstractRNG)

Draw samples from the spatiotemporal GMRF at each time point.

Returns

  • A vector of sample values of length Nₜ, one sample value vector for each time point.
source
GaussianMarkovRandomFields.discretization_at_time Function
julia
discretization_at_time(::AbstractSpatiotemporalGMRF, t::Int)

Return the spatial discretization at time t.

source

Utilities ​

GaussianMarkovRandomFields.spatial_to_spatiotemporal Function
julia
spatial_to_spatiotemporal(
    spatial_matrix::AbstractMatrix,
    t_idx::Int,
    N_t::Int,
)

Make a spatial matrix applicable to a spatiotemporal system at time index t_idx. Results in a matrix that selects the spatial information exactly at time t_idx.

Arguments

  • spatial_matrix::AbstractMatrix: The spatial matrix.

  • t_idx::Integer: The time index.

  • N_t::Integer: The number of time points.

source
GaussianMarkovRandomFields.kronecker_product_spatiotemporal_model Function
julia
kronecker_product_spatiotemporal_model(
    Q_t::AbstractMatrix,
    Q_s::AbstractMatrix,
    spatial_disc::FEMDiscretization;
    algorithm = nothing,
)

Create a spatiotemporal GMRF through a Kronecker product of the temporal and spatial precision matrices.

Arguments

  • Q_t::AbstractMatrix: The temporal precision matrix.

  • Q_s::AbstractMatrix: The spatial precision matrix.

  • spatial_disc::FEMDiscretization: The spatial discretization.

Keyword arguments

  • algorithm: The LinearSolve algorithm to use.
source