Matérn Prior Fields¶

prior ¶

Hippylib Bilaplacian prior fields for vector-valued functions.

This module implements Gaussian prior fields for non-parametric Bayesian inference with Hippylib. It extends the Hippylib SqrtPrecisionPDEPrior class to vector-valued fields, where the individual component fields are non-stationary, but statistically independent. We define the prior for a component with a mean function \(\bar{m}\) and a precision operator \(\mathcal{R}\) as

\[ \pi_{\text{prior}} = \mathcal{N}(\bar{m}, \mathcal{R}^{-1}). \]

For optimization-based inference, we define the negative log distribution of the prior as a cost functional. With a given precision matrix \(\mathbf{R}\) and mean vector \(\bar{\mathbf{m}}\), the discretized cost functional reads

\[ J_{\text{prior}}(\mathbf{m}) = \frac{1}{2} ||(\mathbf{m} - \bar{\mathbf{m}})||_\mathbf{R}^2. \]

The prior object provides methods for evaluating this functional, as well as its gradient and Hessian-vector product. In addition, we can draw samples from the distribution. We focus on a special class of priors with Matérn or Matérn-like covariance structure. We exploit the correspondence between such fields and the solution of SPDEs with specific left-hand-side operators, as proposed here. This allows for the definition of a sparsity-promoting precision operator, in our case one reminiscent of a bilaplacian operator. For a domain \(\Omega\subset\mathbb{R}^d\), we write

\[ \mathcal{R} = \left(\delta - \gamma \Delta\right)^2,\ x\in\Omega. \]

The parameters \(\gamma\) and \(\delta\) can be spatially varying, and have direct correspondence (neglecting boundary effects) to the variance and correlation length of the prior field,

\[ \sigma^2 = \frac{\Gamma(\nu)}{(4\pi^{d/2})}\frac{1}{\delta^\nu\gamma^{d/2}},\quad \rho = \sqrt{\frac{8\nu\gamma}{\delta}},\quad \nu = 2 - \frac{d}{2}. \]

Lastly, to mitigate boundary effects, we can apply Robin boundary conditions to the prior field,

\[ \mathcal{R} = \gamma \nabla m \cdot \mathbf{n} + \frac{\sqrt{\delta\gamma}}{c} m,\ x\in\partial\Omega, \]

with an empirically optimized constant \(c=1.42\).

Classes:

Name	Description
`SqrtPrecisionPDEPrior`	Re-implementation of Hippylib's `SqrtPrecisionPDEPrior` prior.
`PriorSettings`	Settings for the setup of a bilaplacian vector prior.
`Prior`	SPIN prior object, wrapping the Hippylib object with addiaional data and functionality.
`BilaplacianVectorPriorBuilder`	Builder for vector-valued Bilaplacian priors.

spin.hippylib.prior.SqrtPrecisionPDEPrior ¶

Bases: hl.prior._Prior

Re-implementation of Hippylib's SqrtPrecisionPDEPrior prior.

This class is an almost one-to-one re-implementation of the SqrtPrecisionPDEPrior class from Hippylib. It introduces some minor improvements to the original design, and, more importantly, has better support for vector-valued prior fields.

On the component level, this class implements Gaussian prior fields, whose precision operator is given as a bilaplacian-like operator, reading for a domain \(\Omega\subset\mathbb{R}^d\):

\[ R_i= \big(\delta_i(x) - \gamma_I(x) \Delta\big)^2,\ \mathbf{x}\in\Omega, \]

where the parameters \(\delta_i(x)\) and \(\gamma_i(x)\) can be spatially varying. Vector-valued priors treat components as statistically independent, meaning that the precision is block- structured.

According to the Hippylib _Prior nominal subtyping, this class has methods for evaluation of the prior cost functional (the negative log distribution), its parametric gradient, and the action of the Hessian (which is the precision operator). It further allows for sampling from the prior field via a specialized, sparse Cholesky decomposition.

Attributes:

Name	Type	Description
`mean`	`dl.Vector \| dl.PETScVector`	Mean vector of the prior field (Hippylib interface)
`M`	`dl.Matrix`	Mass matrix of the prior field (Hippylib interface)
`Msolver`	`dl.PETScKrylovSolver`	Solver for the mass matrix (Hippylib interface)
`R`	`hl.prior._BilaplacianR`	Precision operator action of the prior field (Hippylib interface)
`Rsolver`	`hl.prior._BilaplacianRsolver`	Covariance operator action of the prior field (Hippylib interface)

Methods:

Name	Description
`init_vector`	Initialize a vector for prior computations
`sample`	Draw a sample from the prior field

init ¶

__init__(function_space: dl.FunctionSpace, variational_form_handler: Callable[[ufl.Argument, ufl.Argument], ufl.Form], mean: dl.Vector | dl.PETScVector, cg_solver_relative_tolerance: Annotated[float, Is[lambda x: 0 < x < 1]] = 1e-12, cg_solver_max_iter: Annotated[int, Is[lambda x: x > 0]] = 1000) -> None

Constructor, building the prior internally.

The Prior works on vector function spaces, with the corresponding build-up in the BilaplacianVectorPriorBuilder class.

Strictly speaking, the size of the constructor is a violation of good design principle, in the sense that this class is its own builder. However, we stick to this pattern according to the Hippylib interface, and sub-divide the constructor into a sequence of smaller methods for clarity.

Parameters:

Name	Type	Description	Default
`function_space`	`dl.FunctionSpace`	Scalar function space.	required
`variational_form_handler`	`Callable[[ufl.Argument, ufl.Argument], ufl.Form]`	variational form describing the underlying prior field as an SPDE.	required
`mean`	`dl.Vector \| dl.PETScVector`	Mean vector.	required
`cg_solver_relative_tolerance`	`int`	Relative tolerance for CG solver for matrix-free inversion. Defaults to 1e-12.	`1e-12`
`cg_solver_max_iter`	`int`	Maximum number of iterations for CG solver for matrix-free inversion. Defaults to 1000.	`1000`

Raises:

Type	Description
`ValueError`	Checks if the mean vector has the same dimension as the function space.

_assemble_system_matrices `staticmethod` ¶

_assemble_system_matrices(trial_function: ufl.Argument, test_function: ufl.Argument, variational_form_handler: Callable) -> tuple[dl.Matrix, dl.Matrix]

Assemble the mass and bilaplacian operator matrix for the FEM system.

Parameters:

Name	Type	Description	Default
`trial_function`	`ufl.Argument`	FEM trial function.	required
`test_function`	`ufl.Argument`	FEM test function.	required
`variational_form_handler`	`Callable`	UFL form callable for the bilaplacian operator.	required

Returns:

Type	Description
`tuple[dl.Matrix, dl.Matrix]`	tuple[dl.Matrix, dl.Matrix]: Mass and operator matrices.

_initialize_solvers ¶

_initialize_solvers(cg_solver_max_iter: int, cg_solver_relative_tolerance: Real) -> tuple[dl.PETScKrylovSolver, dl.PETScKrylovSolver]

Initialize solvers for matrix-free inversion of the mass and bilaplacian operator matrix.

Both matrices are inverted matrix-free using the CG method. For the mass matrix, we employ a simple Jacobi preconditioner, while for the bilaplacian operator matrix, we use the algebraic multigrid method.

Parameters:

Name	Type	Description	Default
`cg_solver_max_iter`	`int`	Maximum number of CG iterations.	required
`cg_solver_relative_tolerance`	`Real`	Relative tolerance for CG termination.	required

Returns:

Type	Description
`tuple[dl.PETScKrylovSolver, dl.PETScKrylovSolver]`	tuple[dl.PETScKrylovSolver, dl.PETScKrylovSolver]: Initialized solver objects, acting as application of the inverse matrices to vectors.

_modify_quadrature_representation `staticmethod` ¶

_modify_quadrature_representation() -> tuple[object, object]

Change UFL form representation for quadrature space.

The change in settings is utlized for the assembly of the mass matrix cholesky factor.

Returns:

Type	Description
`tuple[object, object]`	tuple[object, object]: Buffers holding the default representation settings.

_restore_quadrature_representation `staticmethod` ¶

_restore_quadrature_representation(representation_buffers: tuple[object, object]) -> None

Change UFL representation back to default.

Reverts effects of _modify_quadrature_representation.

Parameters:

Name	Type	Description	Default
`representation_buffers`	`tuple[object, object]`	Buffers holding the default representation settings.	required

_set_up_quadrature_space ¶

_set_up_quadrature_space(quadrature_degree: int) -> tuple[dl.FunctionSpace, ufl.Argument, ufl.Argument]

Set up quadrature space, trial and test function.

Used for the assembly of the mass matrix cholseky factor.

Parameters:

Name	Type	Description	Default
`quadrature_degree`	`int`	Degree of the quadrature representation.	required

Returns:

Type	Description
`tuple[dl.FunctionSpace, ufl.Argument, ufl.Argument]`	tuple[dl.FunctionSpace, ufl.Argument, ufl.Argument]: Quadrature space, trial and test functions.

_assemble_mass_matrix_cholesky_factor ¶

_assemble_mass_matrix_cholesky_factor(quadrature_degree: int, quadrature_space: dl.FunctionSpace, quadrature_trial_function: ufl.Argument, quadrature_test_function: ufl.Argument, test_function: ufl.Argument) -> dl.Matrix

Assemble mass matrix Cholesky factor.

As in Hippylib, we provide a special form of the mass matrix decomposition for sampling from the prior field. The decomposition is rectangular and sparse, based on the quadrature space functionality in Fenics. For more details, we refer to the appendix of the Hippylib paper.

Parameters:

Name	Type	Description	Default
`quadrature_degree`	`int`	Degree of the quadrature representation.	required
`quadrature_space`	`dl.FunctionSpace`	Quadrature function space.	required
`quadrature_trial_function`	`ufl.Argument`	Quadrature space trial function.	required
`quadrature_test_function`	`ufl.Argument`	Quadrature space test function.	required
`test_function`	`ufl.Argument`	Test function on original space.	required

Returns:

Type	Description
`dl.Matrix`	dl.Matrix: Sparse Cholesky factor of the mass matrix

_set_up_hippylib_interface ¶

_set_up_hippylib_interface() -> None

Assign internal attributes to public properties, to adhere to Hippylib interface.

init_vector ¶

init_vector(vector_to_init: dl.Vector, matrix_dim: int | str) -> None

Initialize a vector to be used for prior computations.

Parameters:

Name	Type	Description	Default
`vector_to_init`	`dl.Vector`	Dolfin vector to initialize.	required
`matrix_dim`	`int \| str`	Matrix "dimension" to initialize the vector for.	required

sample ¶

sample(noise_vector: dl.Vector, matern_field_vector: dl.Vector, add_mean: bool = True) -> None

Draw a sample from the prior field.

Parameters:

Name	Type	Description	Default
`noise_vector`	`dl.Vector`	Noise vector for the sample.	required
`matern_field_vector`	`dl.Vector`	Resulting sample vector.	required
`add_mean`	`bool`	Whether to add prior mean vector to sample. Defaults to True.	`True`

spin.hippylib.prior.PriorSettings `dataclass` ¶

Settings for the setup of a bilaplacian vector prior.

Attributes:

Name	Type	Description
`function_space`	`dl.FunctionSpace`	Function space the prior is defined on.
`mean`	`Iterable[str]`	Mean functions, given as a sequence of strings in dolfin syntax, will be compiled to dolfin expressions.
`variance`	`Iterable[str]`	Variance functions, given as a sequence of strings in dolfin syntax, will be compiled to dolfin expressions.
`correlation_length`	`Iterable[str]`	Correlation length functions, given as a sequence of strings in dolfin syntax, will be compiled to dolfin expressions.
`anisotropy_tensor`	`Iterable[hl.ExpressionModule.AnisTensor2D]`	Anisitropy tensor for anisotropic covariance structure. Defaults to None.
`cg_solver_relative_tolerance`	`float`	Tolerance of CG solver for matrix-free application of inverse operators. Defaults to 1e-12.
`cg_solver_max_iter`	`int`	Maximum iteration number of CG solver for matrix-free application of inverse operators. Defaults to 1000.
`robin_bc`	`bool`	Whether to apply Robin boundary conditions. Defaults to False.
`robin_bc_const`	`Real`	Constant for the Robin boundary condition. Defaults to 1.42.

spin.hippylib.prior.Prior `dataclass` ¶

SPIN prior object, wrapping the Hippylib object with additional data and functionality.

Attributes:

Name	Type	Description
`hippylib_prior`	`hl.prior._Prior`	Hippylib prior object
`function_space`	`dl.FunctionSpace`	Underlying function space
`mean_array`	`npt.NDArray[np.floating]`	Mean as numpy array
`variance_array`	`npt.NDArray[np.floating]`	Ideal pointwise variance (without boundary effects) as numpy array
`correlation_length_array`	`npt.NDArray[np.floating]`	Ideal correlation length (without boundary effects) as numpy array
`spde_matern_matrix`	`sp.sparse.coo_array`	Matern SPDE matrix in COO format
`mass_matrix`	`sp.sparse.coo_array`	Mass matrix in COO format

Methods:

Name	Description
`compute_variance_with_boundaries`	Approximate the actual variance field with boundary effects.
`compute_precision_with_boundaries`	Compute the actual precision matrix with boundary effects (only for testing and development).
`evaluate_gradient`	Evaluate the parametric gradient of the prior cost functional.

compute_variance_with_boundaries ¶

compute_variance_with_boundaries(method: Annotated[str, Is[lambda x: x in Exact, Randomized]], num_eigenvalues_randomized: Annotated[int, Is[lambda x: x > 0]] | None = None) -> npt.NDArray[np.floating]

Compute the actual variance of the prior field.

Opposed to the prescribed, theoretical prior variance, this method evaluates the actual variance for the discretized prior field on a bounded domain. The variance can be computed exactly or approximated by a randomized method. The exact method should only be used for small test problems. The randomized method is based on a truncated eigendecomposition of the covariance matrix, and can therefore be used for larger problems.

Parameters:

Name	Type	Description	Default
`method`	`str`	Method for variance computation, either "`Exact`" or "`Randomized`".	required
`num_eigenvalues_randomized`	`int`	Number of eigenvalues to use for randomized diagonal estimation. Only necessary for option "Randomized". Defaults to None.	`None`

Raises:

Type	Description
`ValueError`	Checks that a number of eigenvalues is provided for the "Randomized" method.

Returns:

Type	Description
`npt.NDArray[np.floating]`	npt.NDArray[np.floating]: Pointwise variance of the prior field..

compute_precision_with_boundaries ¶

compute_precision_with_boundaries() -> npt.NDArray[np.floating]

Compute the entire precision matrix of the prior field.

Should only be used for small test problems.

Returns:

Type	Description
`npt.NDArray[np.floating]`	npt.NDArray[np.floating]: Prior field precision matrix.

evaluate_gradient ¶

evaluate_gradient(parameter_array: npt.NDArray[np.floating]) -> npt.NDArray[np.floating]

Evaluate parametric gradient of the prior cost functional (negative log distribution).

Parameters:

Name	Type	Description	Default
`parameter_array`	`npt.NDArray[np.floating]`	Parameter value for which to compute the gradient.	required

Returns:

Type	Description
`npt.NDArray[np.floating]`	npt.NDArray[np.floating]: Gradient array

spin.hippylib.prior.BilaplacianVectorPriorBuilder ¶

Builder for vector-valued Bilaplacian priors.

This builder assembles a vector-valued prior field, with bilaplacian-like precision operator. It constructs component-wise variational forms and initializes an SqrtPrecisionPDEPrior object with the resulting data. Different variance and correlation structures can be supplied for each component, but the components themselves are statistically independent of each other.

Methods:

Name	Description
`build`	Main interface of the builder.

init ¶

__init__(prior_settings: PriorSettings) -> None

Constructor, internally initializes data structures.

Parameters:

Name	Type	Description	Default
`prior_settings`	`PriorSettings`	Configuration and data for the prior field.	required

build ¶

build() -> Prior

Main interface of the builder.

Returns a Prior wrapper object, containing the Hippylib prior object alongside additional data and functionalities.

Returns:

Name	Type	Description
`Prior`	`Prior`	SPIN Prior object.

_convert_prior_coefficients ¶

_convert_prior_coefficients() -> tuple[dl.Function, dl.Function]

Convert variance and correlation length fields to prior coefficients.

For the Bilaplacian prior, the variance \(\sigma^2\) and correlation length \(\rho\) can be converted to the prior SPDE parameters \(\gamma\) and \(\delta\) as

\[ \begin{gather*} \nu = 2- \frac{d}{2}, \quad \kappa = \frac{\sqrt{8\nu}}{\rho}, \\ s = \sigma \kappa^\nu \sqrt{\frac{4\pi^{d/2}}{\Gamma(\nu)}}, \\ \gamma = \frac{1}{s}, \quad \delta = \frac{\kappa^2}{s}. \end{gather*} \]

Returns:

Type	Description
`tuple[dl.Function, dl.Function]`	tuple[dl.Function, dl.Function]: Coefficients \(\gamma\) and \(\delta\) for the prior field.

_generate_scalar_form ¶

_generate_scalar_form(trial_function: ufl.Argument, test_function: ufl.Argument, index: int) -> ufl.Form

Generate the compomnent-wise UFL form for the Bilaplacian prior.

The component form has three different contributions. Given a trial function \(u\) and test function \(v\), this contributions are:

The mass matrix term, simply the discretization of \(\int_{\Omega}\delta u v d\mathbf{x}\)
The stiffness matrix term, discretizing \(\int_{\Omega}\gamma \nabla u \cdot \nabla v d\mathbf{x}\)
The Robin boundary term, discretizing \(\int_{\partial\Omega} \frac{\sqrt{\delta\gamma}}{c} u v dx\)

The boundary term is only applied if Robin boundary conditions are set to true with a boundary constant \(c\). The stiffness matrix term can additionally contain an anisotropy tensor, if the user provides one.

Parameters:

Name	Type	Description	Default
`trial_function`	`ufl.Argument`	FE Trial function.	required
`test_function`	`ufl.Argument`	FE Test function.	required
`index`	`int`	Index of the component being discretized in a vector context.	required

Returns:

Type	Description
`ufl.Form`	ufl.Form: UFL form for the component-wise discretization.

Matérn Prior Fields¶

prior ¶

spin.hippylib.prior.SqrtPrecisionPDEPrior ¶

__init__ ¶

_assemble_system_matrices staticmethod ¶

_initialize_solvers ¶

_modify_quadrature_representation staticmethod ¶

_restore_quadrature_representation staticmethod ¶

_set_up_quadrature_space ¶

_assemble_mass_matrix_cholesky_factor ¶

_set_up_hippylib_interface ¶

init_vector ¶

sample ¶

spin.hippylib.prior.PriorSettings dataclass ¶

spin.hippylib.prior.Prior dataclass ¶

compute_variance_with_boundaries ¶

compute_precision_with_boundaries ¶

evaluate_gradient ¶

spin.hippylib.prior.BilaplacianVectorPriorBuilder ¶

__init__ ¶

build ¶

_convert_prior_coefficients ¶

_generate_scalar_form ¶

init ¶

_assemble_system_matrices `staticmethod` ¶

_modify_quadrature_representation `staticmethod` ¶

_restore_quadrature_representation `staticmethod` ¶

spin.hippylib.prior.PriorSettings `dataclass` ¶

spin.hippylib.prior.Prior `dataclass` ¶

init ¶