Gaussian Approximation

Provided by GaussianMarkovRandomFields.jl

The gaussian_approximation function is implemented in GaussianMarkovRandomFields.jl v0.4+ and re-exported by Latte.jl for user convenience. This documentation describes the re-exported functionality.

The Gaussian approximation functionality provides efficient Newton-Raphson optimization for finding posterior modes in INLA. This is the core computational engine for constructing Gaussian approximations to non-Gaussian posteriors.

Overview

The Gaussian approximation process finds the mode of the posterior distribution p(x|y) and constructs a Gaussian approximation around it using Fisher scoring (Newton-Raphson with Fisher information matrix). This approximation forms the foundation of the INLA methodology.

Main Function

GaussianMarkovRandomFields.gaussian_approximation Function

julia

gaussian_approximation(prior_gmrf, obs_lik; kwargs...) -> AbstractGMRF

Find Gaussian approximation to the posterior using Fisher scoring.

This function finds the mode of the posterior distribution and constructs a Gaussian approximation around it using Fisher scoring (Newton-Raphson with Fisher information matrix).

Works for GMRF, ConstrainedGMRF, and ChordalGMRF priors, automatically handling constraint projection when needed.

Arguments

prior_gmrf: Prior GMRF distribution for the latent field (GMRF, ConstrainedGMRF, or ChordalGMRF)
obs_lik: Materialized observation likelihood (contains data and hyperparameters)

Keyword Arguments

max_iter::Int=50: Maximum number of Fisher scoring iterations
mean_change_tol::Real=1e-4: Convergence tolerance for mean change
newton_dec_tol::Real=1e-5: Newton decrement convergence tolerance
adaptive_stepsize::Bool=true: Enable adaptive stepsize with backtracking line search
max_linesearch_iter::Int=10: Maximum line search iterations per Newton step
verbose::Bool=false: Print iteration information

Returns

Gaussian approximation to the posterior p(x | θ, y) (same type as input prior)

Example

julia

# Set up components
prior_gmrf = GMRF(μ_prior, Q_prior)
obs_model = ExponentialFamily(Poisson)
obs_lik = obs_model(y)

# Find Gaussian approximation - uses adaptive stepsize by default
posterior_gmrf = gaussian_approximation(prior_gmrf, obs_lik)

# For well-conditioned problems, disable adaptive stepsize for speed
posterior_gmrf = gaussian_approximation(prior_gmrf, obs_lik; adaptive_stepsize=false)

source

julia

gaussian_approximation(prior::WorkspaceGMRF, obs_lik::ObservationLikelihood; kwargs...)

Workspace-aware Gaussian approximation via Fisher scoring. Uses the workspace's factorisation engine for numeric-only refactorisation on each Newton step.

source

julia

gaussian_approximation(prior::LatentModel, obs_lik; θ, ws=nothing, kwargs...)

Materialise the Gaussian latent prior at θ and delegate to the (::AbstractGMRF, obs_lik) / (::WorkspaceGMRF, obs_lik) dispatch for fixed-Q Newton.

source

julia

gaussian_approximation(prior::NonGaussianLatentPrior, obs_lik; θ, ws=nothing, x0, kwargs...)

Iterated-linearisation Gaussian approximation. The Newton loop calls local_quadratic(prior, x_k; θ...) per iterate to re-linearise the prior at the current iterate; the line-search merit uses the exact log p(x | θ) carried in LocalLatentQuadratic.logp_ref.

Hyperparameter values may be passed splatted (τ = 1.0) or as a θ::NamedTuple keyword; both forms merge before dispatch. Non-Gaussian priors don't have a canonical mean, so x0 defaults to zeros(length(prior)) — pass it explicitly for problems where zero is a poor starting point (e.g. priors with reflection symmetry through zero).

When a ws::GMRFWorkspace is supplied, it must be seeded with the full structural sparsity pattern that local_quadratic(prior, x; θ...) produces across all Newton iterates (the workspace reuses one symbolic factorisation, and the per-iterate precision values are copied onto it positionally). Seeding from a generic, non-degenerate x is the safe way to capture every structural coupling.

source

Implementation Details

The gaussian_approximation function uses NonlinearSolve.jl with Newton-Raphson optimization and CHOLMOD factorization for efficient sparse linear algebra.

Basic Usage

julia

using Latte
using GaussianMarkovRandomFields
using Distributions

# Set up prior GMRF
μ_prior = zeros(10)
Q_prior = spdiagm(0 => ones(10))
prior_gmrf = GMRF(μ_prior, Q_prior)

# Set up observation model
obs_model = ExponentialFamily(Poisson)
θ_named = NamedTuple()  # No hyperparameters for Poisson

# Generate synthetic data
x_true = rand(prior_gmrf)
data_dist = conditional_distribution(obs_model, x_true; θ_named...)
y_obs = rand(data_dist)

# Find Gaussian approximation
obs_lik = obs_model(y_obs; θ_named...)
posterior_gmrf = gaussian_approximation(prior_gmrf, obs_lik)

# Extract posterior statistics
posterior_mean = mean(posterior_gmrf)
posterior_precision = precision_matrix(posterior_gmrf)
println("Gaussian approximation computed successfully")

Mathematical Properties

The Gaussian approximation finds the mode of the posterior distribution and constructs a Gaussian around it. For observation models with Gaussian likelihood, this approximation is exact:

julia

# For Gaussian observations, the approximation is exact
obs_model = ExponentialFamily(Normal)
θ_named = (σ = 0.5,)

# The posterior precision combines prior and observation precision
# Q_posterior = Q_prior + Q_obs
# μ_posterior = Q_posterior^(-1) * (Q_prior * μ_prior + Q_obs * y)

Performance Considerations

Fisher Information vs Hessian

The implementation uses Fisher scoring, which approximates the Hessian with the Fisher information matrix. This has several advantages:

Stability: Fisher information is always positive semi-definite
Efficiency: Often faster convergence than pure Newton-Raphson
Numerical robustness: Less sensitive to poor conditioning

Sparse Linear Algebra

The implementation leverages sparse linear algebra throughout:

Sparse precision matrices are preserved during optimization
Cholesky factorizations use sparse solvers when available
Memory usage scales with the sparsity pattern of the prior

Optimization Details

The implementation uses NonlinearSolve.jl with Newton-Raphson method:

Jacobian computation: Uses automatic differentiation for precise gradients and Hessians
Linear solver: CHOLMOD factorization for sparse precision matrices
Convergence: Built-in tolerances for abstol=1e-6 and reltol=1e-6
Permutation handling: Preserves sparsity patterns from the prior GMRF

Mathematical Background

Fisher Scoring Update

At each iteration, the algorithm performs the update:

Q_new = Q_prior - Hessian_obs(μ_current)
gradient = Q_prior * (μ_current - μ_prior) - grad_obs(μ_current)
μ_new = μ_current - Q_new⁻¹ * gradient

Where:

Q_prior is the prior precision matrix
μ_prior is the prior mean
Hessian_obs is the Hessian of the observation log-likelihood
grad_obs is the gradient of the observation log-likelihood