Observation Models

Observation models define the relationship between observations y and the latent GMRF field x, typically through likelihood functions. They enable Bayesian inference by connecting your data to the underlying Gaussian process through flexible probabilistic models.

GaussianMarkovRandomFields.jl implements observation models using a factory pattern that separates model configuration from materialized evaluation instances. This design provides major performance benefits in optimization loops and cleaner automatic differentiation boundaries.

Core Concepts

The Factory Pattern

Observation models follow a two-stage pattern:

ObservationModel: A factory that defines the model structure and hyperparameters
ObservationLikelihood: A materialized instance with specific data and hyperparameters for fast evaluation

julia

# Step 1: Configure observation model (factory)
obs_model = ExponentialFamily(Normal)

# Step 2: Materialize with data and hyperparameters  
obs_lik = obs_model(y; σ=1.2)

# Step 3: Fast evaluation in hot loops
ll = loglik(x, obs_lik)      # Only x argument needed!
grad = loggrad(x, obs_lik)   # Fast x-only evaluation
hess = loghessian(x, obs_lik)

This pattern eliminates the need to repeatedly pass data and hyperparameters, providing significant performance benefits in optimization and sampling algorithms.

Evaluation Interface

All materialized observation likelihoods support a common interface:

loglik(x, obs_lik): Evaluate log-likelihood
loggrad(x, obs_lik): Compute gradient with respect to latent field
loghessian(x, obs_lik): Compute Hessian matrix

Exponential Family Models

The most common observation models are exponential family distributions connected to the latent field through link functions.

Basic Usage

julia

using GaussianMarkovRandomFields
using Distributions

# Poisson model for count data (canonical LogLink)
poisson_model = ExponentialFamily(Poisson)
x = [1.0, 2.0]  # Latent field (log-intensity due to LogLink)
y = [2, 7]      # Count observations
obs_lik = poisson_model(y)
ll = loglik(x, obs_lik)

# Normal model for continuous data (canonical IdentityLink)
normal_model = ExponentialFamily(Normal)
x = [1.5, 2.3]  # Latent field (direct mean due to IdentityLink)
y = [1.2, 2.8]  # Continuous observations
obs_lik = normal_model(y; σ=0.5)  # Normal requires σ hyperparameter
ll = loglik(x, obs_lik)

# Bernoulli model for binary data (canonical LogitLink)
bernoulli_model = ExponentialFamily(Bernoulli)
x = [0.0, 1.5]  # Latent field (logit-probability due to LogitLink)
y = [0, 1]      # Binary observations
obs_lik = bernoulli_model(y)
ll = loglik(x, obs_lik)

Supported Distributions and Links

Distribution	Canonical Link	Alternative Links	Hyperparameters
Normal	IdentityLink	LogLink	σ (std. dev.)
Poisson	LogLink	IdentityLink	none
Bernoulli	LogitLink	LogLink	none
Binomial	LogitLink	IdentityLink	none*

*For Binomial, the number of trials is provided through the data structure BinomialObservations, not as a hyperparameter.

Custom Link Functions

julia

# Non-canonical link function
poisson_identity = ExponentialFamily(Poisson, IdentityLink())
# Note: Requires positive latent field values for valid Poisson intensities

Custom Observation Models

For models not covered by exponential families, you can define custom log-likelihood functions using automatic differentiation.

Basic AutoDiff Models

julia

# Define custom log-likelihood function
function custom_loglik(x; y=[1.0, 2.0], σ=1.0)
    μ = sin.(x)  # Custom transformation
    return -0.5 * sum((y .- μ).^2) / σ^2 - length(y) * log(σ)
end

# Create observation model
obs_model = AutoDiffObservationModel(custom_loglik; n_latent=2, hyperparams=(:y, :σ))

# Materialize with data
obs_lik = obs_model(y=[1.2, 1.8], σ=0.5)

# Use normally - gradients and Hessians computed automatically!
x = [0.5, 1.0]
ll = loglik(x, obs_lik)
grad = loggrad(x, obs_lik)    # Automatic differentiation
hess = loghessian(x, obs_lik) # Potentially sparse!

Automatic Differentiation Requirements

AutoDiff observation models require an automatic differentiation backend. We support and recommend the following backends in order of preference:

Enzyme.jl (recommended for performance)
Mooncake.jl (good balance of performance and compatibility)
Zygote.jl (reliable fallback)
ForwardDiff.jl (for small problems)

julia

# Load an AD backend (required for AutoDiffObservationModel)
using Enzyme  # Recommended

# Or use another supported backend:
# using Mooncake
# using Zygote
# using ForwardDiff

# Now you can use AutoDiff models
obs_model = AutoDiffObservationModel(my_loglik; n_latent=10)
obs_lik = obs_model(y=data)
grad = loggrad(x, obs_lik)  # Uses your loaded AD backend

Sparse Hessian Computation

AutoDiff observation models can automatically detect and exploit sparsity in Hessian matrices using our package extensions. This requires loading both an AD backend and additional sparsity packages:

julia

# Load AD backend + sparse AD packages
using Enzyme  # Or your preferred AD backend
using SparseConnectivityTracer, SparseMatrixColorings

# The package extension is automatically activated
obs_model = AutoDiffObservationModel(my_loglik; n_latent=100)
obs_lik = obs_model(y=data)

# Hessian computation now automatically:
# - Detects sparsity pattern using TracerSparsityDetector  
# - Uses greedy coloring for efficient computation
# - Returns sparse matrix when beneficial
hess = loghessian(x, obs_lik)  # May be sparse!

The sparse Hessian features provide dramatic performance improvements for large-scale problems with structured sparsity.

Nonlinear Least Squares

Use when observations are Gaussian with mean given by an arbitrary (possibly nonlinear) function of the latent field: y | x ~ Normal(f(x), σ).

Key properties

Out-of-place f: define f(x)::AbstractVector with length equal to length(y).
Gauss–Newton: gradient and Hessian use the Gauss–Newton approximation (no exact Hessian term).
- ∇ℓ(x) = J(x)' (w ⊙ r), where r = y − f(x), w = 1 ./ σ.^2
- ∇²ℓ(x) ≈ − J(x)' Diagonal(w) J(x)
σ: accepts a scalar or vector (heteroskedastic), both interpreted as standard deviations.
Sparse autodiff: requires loading SparseConnectivityTracer and SparseMatrixColorings to activate the sparse Jacobian backend.

Example

julia

using GaussianMarkovRandomFields
using SparseConnectivityTracer, SparseMatrixColorings  # activate sparse Jacobian backend

# Nonlinear mapping f: R^2 -> R^3
f(x) = [x[1] + 2x[2], sin(x[1]), x[2]^2]

# Observations and noise
y = [1.0, 0.5, 2.0]
σ = [0.3, 0.4, 0.5]  # vector sigma allowed

# Build model and materialize likelihood
model = NonlinearLeastSquaresModel(f, 2)
lik = model(y; σ=σ)

# Evaluate
x = [0.1, 0.2]
ll = loglik(x, lik)
g  = loggrad(x, lik)     # uses sparse DI.jacobian under the hood
H  = loghessian(x, lik)  # Gauss–Newton: -J' W J

# Conditional distribution p(y | x)
dist = conditional_distribution(model, x; σ=0.3)

API

GaussianMarkovRandomFields.NonlinearLeastSquaresModel Type

julia

NonlinearLeastSquaresModel(f, n)

Observation model for nonlinear least squares with Gaussian noise: y | x ~ Normal(f(x), σ)

This model uses a Gauss–Newton approximation for the Hessian: ∇ℓ(x) = J(x)' (w ⊙ r), where r = y - f(x), w = 1 ./ σ.^2 ∇²ℓ(x) ≈ -J(x)' Diagonal(w) J(x)

Notes

Requires the sparse-AD extension (SparseConnectivityTracer + SparseMatrixColorings) to be loaded. If missing, construction or evaluation will error with a clear message.
f must be out-of-place with signature f(x)::AbstractVector of the same length as y.
σ can be a scalar or a vector matching length(y) (heteroskedastic case). It must be positive.

Observation Models ​

Core Concepts ​

The Factory Pattern ​

Evaluation Interface ​

Exponential Family Models ​

Basic Usage ​

Supported Distributions and Links ​

Custom Link Functions ​

Custom Observation Models ​

Basic AutoDiff Models ​

Automatic Differentiation Requirements ​

Sparse Hessian Computation ​

Nonlinear Least Squares ​

Advanced Features ​

Linear Transformations and Design Matrices ​

Binomial Observations ​

Composite Observations ​

API Reference ​

Core Types and Interface ​

Exponential Family Models ​

Link Functions ​

Custom AutoDiff Models ​

FEM Helper Functions ​

Advanced Features ​

Observation Models

Core Concepts

The Factory Pattern

Evaluation Interface

Exponential Family Models

Basic Usage

Supported Distributions and Links

Custom Link Functions

Custom Observation Models

Basic AutoDiff Models

Automatic Differentiation Requirements

Sparse Hessian Computation

Nonlinear Least Squares

Advanced Features

Linear Transformations and Design Matrices

Binomial Observations

Composite Observations

API Reference

Core Types and Interface

Exponential Family Models

Link Functions

Custom AutoDiff Models

FEM Helper Functions

Advanced Features