Gaussian Likelihood Models¶

misfit ¶

Misfit/Likelihood for usage in Hippylib.

This module implements the likelihood for Bayesian inverse problems in the Hippylib framework. In the optimization context, it can be interpreted as a misfit functional. We focus on additive, zero-centered Gaussian noise models. Given a PDE solution operator \(\mathcal{G}(\mathbf{m}) = \mathbf{u}\), we define for each output component \(u\) the projection to points of observation \(u_d = \mathcal{B}u \in \mathbb{R}^q\) with observation operator \(\mathcal{B}: \mathcal{U} \to \mathbb{R}^q\). We then write for the component data \(d\) that

\[ d = u_d + \eta, \quad \eta \sim \mathcal{N}(0, \Gamma). \]

Therefore, we can define the likelihood density for a single variable component u as

\[ \pi_{\text{like}}(d|u) \propto \exp\left(-\frac{1}{2} \| \mathcal{B}u - d \|_{\Gamma^{-1}}^2\right). \]

The below methods implement discrete versions of the projection operator, which we refer to as \(\mathbf{B}\), and a diagonal noise precision matrix \(\Gamma^{-1}\) with pointwise varying coefficients. The misfit functional is then given as the negative log-likelihood. We further compute gradients and Hessian-vector products of the misfit w.r.t. \(u\). The misfit functional for a single component is implemented in the DiscreteMisfit class. Multiple misfit objects can be combined in the VectorMisfit wrapper.

Classes:

Name	Description
`MisfitSettings`	Dataclass to store settings for misfit construction.
`DiscreteMisfit`	Misfit class for scalar misfit problems.
`VectorMisfit`	Misfit class for vector-valued misfit problems.
`TDMisfit`	Misfit class for time-dependent misfit problems (not implemented yet).
`Misfit`	Dataclass to store assembled misfit objects.
`MisfitBuilder`	Builder class to assemble misfit objects from settings.

Functions:

Name	Description
`assemble_pointwise_observation_operator`	Assemble pointwise observation operator.
`assemble_noise_precision_matrix`	Assemble noise precision matrix.

spin.hippylib.misfit.DiscreteMisfit ¶

Bases: hl.Misfit

Implementation of a single-component misfit functional, as described above.

Methods:

Name	Description
`cost`	Evaluate the cost of the misfit functional.
`grad`	Compute the gradient of the misfit functional.
`setLinearizationPoint`	Set point for Hessian evaluation.
`apply_ij`	Apply Hessian-vector product.

observation_matrix `property` ¶

observation_matrix: dl.Matrix

Return internally stored observation matrix.

init ¶

__init__(data: npt.NDArray[np.floating], observation_matrix: dl.Matrix, noise_precision_matrix: dl.PETScMatrix) -> None

Constructor, initializes internal data structures and buffers.

Parameters:

Name	Type	Description	Default
`data`	`npt.NDArray[np.floating]`	Observed data.	required
`observation_matrix`	`dl.Matrix`	Projection matrix from function to observation space.	required
`noise_precision_matrix`	`dl.PETScMatrix`	Precision of the noise model.	required

cost ¶

cost(state_list: Iterable[dl.Vector | dl.PETScVector, dl.Vector | dl.PETScVector | None, dl.Vector | dl.PETScVector | None]) -> Real

Evaluate the cost of the misfit functional.

For a given state \(u\), data \(d\), projection matrix \(\mathbf{B}\) and noise precision \(\Gamma^{-1}\), the cost is given as

\[ J_{\text{misfit}}(u) = \frac{1}{2} \| \mathbf{B} u - d \|_{\Gamma^{-1}}^2. \]

Parameters:

Name	Type	Description	Default
`state_list`	`Iterable`	List of forward, parameter and adjoint variables \((u,m,p)\), only \(u\) is required.	required

Returns:

Name	Type	Description
`Real`	`Real`	Cost value for given state.

grad ¶

grad(derivative_type: Annotated[int, IsEqual[hl.STATE] | IsEqual[hl.PARAMETER]], state_list: Iterable[dl.Vector | dl.PETScVector, dl.Vector | dl.PETScVector | None, dl.Vector | dl.PETScVector | None], output_vector: dl.Vector | dl.PETScVector) -> None

Gradient of the cost functional given in the cost method.

Parameters:

Name	Type	Description	Default
`derivative_type`	`Annotated[int, IsEqual[hl.STATE] \| IsEqual[hl.PARAMETER]]`	Variable with respect to which the gradient is computed. Only `STATE` yields a non-zero result.	required
`state_list`	`Iterable`	List of forward, parameter and adjoint variables \((u,m,p)\), only \(u\) is required.	required
`output_vector`	`dl.Vector \| dl.PETScVector`	Computed gradient.	required

setLinearizationPoint ¶

setLinearizationPoint(_state_list: Iterable[dl.Vector | dl.PETScVector | None, dl.Vector | dl.PETScVector | None, dl.Vector | dl.PETScVector | None], _gauss_newton_approx: bool) -> None

Set point for Hessian evaluation.

This method does nothing, as the misfit is a quadratic form and its second variation is thus constant. Only required for interface compatibility.

apply_ij ¶

apply_ij(first_derivative_type: Annotated[int, IsEqual[hl.STATE] | IsEqual[hl.PARAMETER]], second_derivative_type: Annotated[int, IsEqual[hl.STATE] | IsEqual[hl.PARAMETER]], hvp_direction: dl.Vector | dl.PETScVector, output_vector: dl.Vector | dl.PETScVector) -> None

Apply Hessian vector-product.

Second derivatives are only computed with respect to the forward variable \(u\).

Parameters:

Name	Type	Description	Default
`first_derivative_type`	`Annotated[int, IsEqual[hl.STATE] \| IsEqual[hl.PARAMETER]]`	Variable to compute first derivative with respect to. Only `STATE` yields a non-zero result.	required
`second_derivative_type`	`Annotated[int, IsEqual[hl.STATE] \| IsEqual[hl.PARAMETER]]`	Variable to compute second derivative with respect to. Only `STATE` yields a non-zero result.	required
`hvp_direction`	`dl.Vector \| dl.PETScVector`	Direction vector for which to compute Hessian-vector product.	required
`output_vector`	`dl.Vector \| dl.PETScVector`	Resulting vector.	required

spin.hippylib.misfit.VectorMisfit ¶

Bases: hl.Misfit

Misfit object for vector-valued functions.

The vector misfit is basically a wrapper to a collection of single-component misfiit objects, i.e. the DiscreteMisfit class. For vector-valued quantities, it distributes computations for cost, gradient, and Hessian-vector product to the individual components, and recombines the results accordingly. This makes the vector misfit completely independent of the underlying implementation for individual components.

Info

The only underlying assumption at the moment is that the different components are quadratic forms of the forward variable \(u\).

Methods:

Name	Description
`cost`	Evaluate the cost of the misfit functional.
`grad`	Compute the gradient of the misfit functional.
`setLinearizationPoint`	Set point for Hessian evaluation.
`apply_ij`	Apply Hessian-vector product.

init ¶

__init__(misfit_list: Iterable[hl.Misfit], function_space: dl.FunctionSpace) -> None

Constructor, initialize list of misfits and buffers.

Parameters:

Name	Type	Description	Default
`misfit_list`	`Iterable[hl.Misfit]`	List of misfit objects for each component.	required
`function_space`	`dl.FunctionSpace`	Vector function space of all components.	required

Raises:

Type	Description
`ValueError`	Checks that number of misfits matches number of components in function space.

cost ¶

cost(state_list: Iterable[dl.Vector | dl.PETScVector, dl.Vector | dl.PETScVector | None, dl.Vector | dl.PETScVector | None]) -> float

Evaluate the cost of the misfit functional.

The vector misfit wrapper computes the cost for each component individually and sums up the contributions.

Parameters:

Name	Type	Description	Default
`state_list`	`Iterable`	List of forward, parameter and adjoint variables \((u,m,p)\), only \(u\) is required.	required

Returns:

Name	Type	Description
`Real`	`float`	Cost value for given state.

grad ¶

grad(derivative_type: Annotated[int, IsEqual[hl.STATE] | IsEqual[hl.PARAMETER]], state_list: Iterable[dl.Vector | dl.PETScVector, dl.Vector | dl.PETScVector | None, dl.Vector | dl.PETScVector | None], output_vector: dl.Vector | dl.PETScVector) -> None

Gradient of the cost functional given in the cost method.

The wrapper computes individual gradient components \(g_i\) and assembles them into a vector \(\mathbf{g} = \text{vec}(g_i)\).

Parameters:

Name	Type	Description	Default
`derivative_type`	`Annotated[int, IsEqual[hl.STATE] \| IsEqual[hl.PARAMETER]]`	Variable with respect to which the gradient is computed. Only `STATE` yields a non-zero result.	required
`state_list`	`Iterable`	List of forward, parameter and adjoint variables \((u,m,p)\), only \(u\) is required.	required
`output_vector`	`dl.Vector \| dl.PETScVector`	Computed gradient.	required

setLinearizationPoint ¶

setLinearizationPoint(_state_list: Iterable[dl.Vector | dl.PETScVector | None, dl.Vector | dl.PETScVector | None, dl.Vector | dl.PETScVector | None], _gauss_newton_approx: bool) -> None

Set point for Hessian evaluation.

This method does nothing, as the misfit is a quadratic form and its second variation is thus constant. Only required for interface compatibility.

apply_ij ¶

apply_ij(first_derivative_type: Annotated[int, IsEqual[hl.STATE] | IsEqual[hl.PARAMETER]], second_derivative_type: Annotated[int, IsEqual[hl.STATE] | IsEqual[hl.PARAMETER]], hvp_direction: dl.Vector | dl.PETScVector, output_vector: dl.Vector | dl.PETScVector) -> None

Apply Hessian vector-product.

Second derivatives are only computed with respect to the forward variable \(u\). The vector wrapper distributes the Hessian-vector product computation to the individual components and recombines the results into a single vector.

Parameters:

Name	Type	Description	Default
`first_derivative_type`	`Annotated[int, IsEqual[hl.STATE] \| IsEqual[hl.PARAMETER]]`	Variable to compute first derivative with respect to. Only `STATE` yields a non-zero result.	required
`second_derivative_type`	`Annotated[int, IsEqual[hl.STATE] \| IsEqual[hl.PARAMETER]]`	Variable to compute second derivative with respect to. Only `STATE` yields a non-zero result.	required
`hvp_direction`	`dl.Vector \| dl.PETScVector`	Direction vector for which to compute Hessian-vector product.	required
`output_vector`	`dl.Vector \| dl.PETScVector`	Resulting vector.	required

_create_buffers `staticmethod` ¶

_create_buffers(function_space: dl.FunctionSpace) -> tuple[Iterable[dl.Vector], Iterable[dl.Vector]]

Create buffers for in-memory computations.

spin.hippylib.misfit.TDMisfit ¶

Bases: hl.Misfit

Time-dependent misfit.

Warning

This class is not implemented yet.

init ¶

__init__(stationary_misfits: Iterable[hl.Misfit], observation_times: Iterable[Real]) -> None

Constructor (raises error).

spin.hippylib.misfit.MisfitSettings `dataclass` ¶

Configuration and data for a misfit object.

Depending on the function space, the observation points, values, and noise variance given, the MisfitBuilder will assemble either a DiscreteMisfit or a VectorMisfit object.

Warning

Settings for time-dependent problems are only stubs, time-dependent misfits are not implemented yet.

Attributes:

Name	Type	Description
`function_space`	`dl.FunctionSpace`	Function space of the forward variable. Can be scalar or vector-valued.
`observation_points`	`npt.NDArray[np.floating] \| Iterable[npt.NDArray[np.floating]]`	Points of observation to project to. Either single array or collection of arrays for multiple components.
`observation_values`	`npt.NDArray \| Iterable[npt.NDArray] \| Iterable[Iterable[npt.NDArray]]`	Observed data. Either single array or collection of arrays for multiple components.
`noise_variance`	`npt.NDArray[np.floating] \| Iterable[npt.NDArray[np.floating]]`	Diagonal elements of the noise precision matrix. Either single array or collection of arrays for multiple components.
`stationary`	`bool`	Whether the misfit is time-dependent.
`observation_times`	`npt.NDArray[np.floating] \| None`	Observation times for time-dependent misfit.

Raises:

Type	Description
`ValueError`	Checks that only single arrays are provided for scalar function spaces.
`ValueError`	Checks that the correct number of arrays is provided for vector function spaces.
`ValueError`	Checks that observations times are provided for time-dependent misfits.

__post_init__ ¶

__post_init__() -> None

Post init.

spin.hippylib.misfit.Misfit `dataclass` ¶

Misfit object returned by the builder.

Wraps the corresponding Hippylib object and scipy representations of the noise precision and observation matrices.

spin.hippylib.misfit.MisfitBuilder ¶

Builder for misfit objects.

The builder takes a MisfitSettings object and builds a misfit object according to the provided configuration. The builder distinguishes between single-component and vector-valued function spaces, as well as stationary and time-dependent misfits (time-dependent misfit are not yet implemented though!).

Methods:

Name	Description
`build`	Main interface to build a misfit object.

init ¶

__init__(settings: MisfitSettings) -> None

Constructor, set data structures from settings.

Parameters:

Name	Type	Description	Default
`settings`	`MisfitSettings`	Configuration object for the builder	required

build ¶

build() -> Misfit

Main interface of the builder.

Creates Misfit object from settings.

Returns:

Name	Type	Description
`Misfit`	`Misfit`	SPIN misfit object wrapper.

_assemble_observation_matrices ¶

_assemble_observation_matrices() -> dl.Matrix | list[dl.Matrix]

Assemble observation matrices.

For a single component, this is just one matrix, for multiple components, it is a list of matrices. Simply calls the assemble_pointwise_observation_operator function.

Returns:

Type	Description
`dl.Matrix \| list[dl.Matrix]`	dl.Matrix \| list[dl.Matrix]: Discretizes observation operator(s).

_assemble_noise_precision_matrices ¶

_assemble_noise_precision_matrices() -> dl.PETScMatrix | list[dl.PETScMatrix]

Assemble noise precision matrices.

For a single component, this is just one matrix, for multiple components, it is a list of matrices. Simply calls the assemble_noise_precision_matrix function.

Returns:

Type	Description
`dl.PETScMatrix \| list[dl.PETScMatrix]`	dl.PETScMatrix \| list[dl.PETScMatrix]: Noise precision matrix/matrices.

_build_misfit ¶

_build_misfit() -> Misfit

Builds the misfit from the assembled noise precision and observation matrices.

The method distinguishes four different cases:

Single component, stationary: Assemble a single DiscreteMisfit object.
Single component, time-dependent (not implemented): Assemble a single DiscreteMisfit object. Use it as per-time-step misfit in a time-dependent misfit object.
Multiple components, stationary: Assemble a list of DiscreteMisfit objects and wrap them in a VectorMisfit object.
Multiple components, time-dependent (not implemented): Assemble a list of DiscreteMisfit objects and wrap them in a VectorMisfit object. Use it as per-time-step misfit in a time-dependent misfit object.

Returns:

Name	Type	Description
`Misfit`	`Misfit`	SPIN misfit object wrapper

_convert_matrices_to_scipy ¶

_convert_matrices_to_scipy(matrices: dl.Matrix | dl.PETScMatrix | Iterable[dl.Matrix | dl.PETScMatrix]) -> sp.sparse.coo_array | Iterable[sp.sparse.coo_array]

Convert sparse dolfin matrices to scipy COO arrays.

Uses the convert_matrix_to_scipy function in the Fenics converter module.

Parameters:

Name	Type	Description	Default
`matrices`	`dl.Matrix \| dl.PETScMatrix \| Iterable[dl.Matrix \| dl.PETScMatrix]`	Dolfin matrices to convert	required

Returns:

Type	Description
`sp.sparse.coo_array \| Iterable[sp.sparse.coo_array]`	sp.sparse.coo_array \| Iterable[sp.sparse.coo_array]: Resulting Scipy matrices

spin.hippylib.misfit.assemble_pointwise_observation_operator ¶

assemble_pointwise_observation_operator(function_space: dl.FunctionSpace, observation_points: npt.NDArray[np.floating]) -> dl.Matrix

Assemble pointwise observation matrix.

This is a simple PEP8-conformant wrapper to Hippylibs assemblePointwiseObservation function. The function takes an input function space and points of observation, to which a function shall be projected. Supports multidimensional domains. In a \(d\)-dimensional domain with \(q_d\) observation points per dimension, the observation_points array should have shape \(q_d\times d\). The resulting projection matrix \(\mathbf{B}\) has shape \(q_d d \times N\), where \(N\) is the number of degrees of freedom in the function space.

Warning

The observation operator works only component-wise, on scalar function spaces.

Parameters:

Name	Type	Description	Default
`function_space`	`dl.FunctionSpace`	Function space of vector to project.	required
`observation_points`	`npt.NDArray[np.floating]`	Points of observation to project to.	required

Returns:

Type	Description
`dl.Matrix`	dl.Matrix: Projection operator/matrix.

spin.hippylib.misfit.assemble_noise_precision_matrix ¶

assemble_noise_precision_matrix(noise_variance: npt.NDArray[np.floating]) -> dl.PETScMatrix

Assemble the precision matrix for a Gaussian noise model.

This function only supports diagonal precision matrices, but allows for pointwise varying values. For \(q\) overall observation points, the noise_variance array should have shape \(q\), and the resulting precision matrix \(\Gamma^{-1}\) will have shape \(q \times q\).

Warning

The precision matrix works only component-wise, on scalar function spaces.

Parameters:

Name	Type	Description	Default
`noise_variance`	`npt.NDArray[np.floating]`	Diagonal elements of the precision matrix.	required

Returns:

Type	Description
`dl.PETScMatrix`	dl.PETScMatrix: Sparse precision matrix in dolfin format.

Gaussian Likelihood Models¶

misfit ¶

spin.hippylib.misfit.DiscreteMisfit ¶

observation_matrix property ¶

__init__ ¶

cost ¶

grad ¶

setLinearizationPoint ¶

apply_ij ¶

spin.hippylib.misfit.VectorMisfit ¶

__init__ ¶

cost ¶

grad ¶

setLinearizationPoint ¶

apply_ij ¶

_create_buffers staticmethod ¶

spin.hippylib.misfit.TDMisfit ¶

__init__ ¶

spin.hippylib.misfit.MisfitSettings dataclass ¶

__post_init__ ¶

spin.hippylib.misfit.Misfit dataclass ¶

spin.hippylib.misfit.MisfitBuilder ¶

__init__ ¶

build ¶

_assemble_observation_matrices ¶

_assemble_noise_precision_matrices ¶

_build_misfit ¶

_convert_matrices_to_scipy ¶

spin.hippylib.misfit.assemble_pointwise_observation_operator ¶

spin.hippylib.misfit.assemble_noise_precision_matrix ¶

observation_matrix `property` ¶

init ¶

init ¶

_create_buffers `staticmethod` ¶

init ¶

spin.hippylib.misfit.MisfitSettings `dataclass` ¶

spin.hippylib.misfit.Misfit `dataclass` ¶

init ¶