hIPPYlib - Inverse Problem PYthon library

Inverse Problem PYthon library

 __        ______  _______   _______   __      __  __  __  __
/  |      /      |/       /       /     /  |/  |/  |/  |
$$ |____  $$$$$$/ $$$$$$$  |$$$$$$$  |$$   /$$/ $$ |$$/ $$ |____
$$        $$ |  $$ |__$$ |$$ |__$$ | $$  \/$$/  $$ |/  |$$
$$$$$$$  |  $$ |  $$    $$/ $$    $$/   $$  $$/   $$ |$$ |$$$$$$$  |
$$ |  $$ |  $$ |  $$$$$$$/  $$$$$$$/     $$$$/    $$ |$$ |$$ |  $$ |
$$ |  $$ | _$$ |_ $$ |      $$ |          $$ |    $$ |$$ |$$ |__$$ |
$$ |  $$ |/ $$   |$$ |      $$ |          $$ |    $$ |$$ |$$    $$/
$$/   $$/ $$$$$$/ $$/       $$/           $$/     $$/ $$/ $$$$$$$/

https://hippylib.github.io

hIPPYlib implements state-of-the-art scalable algorithms for PDE-based deterministic and Bayesian inverse problems. It builds on FEniCS (a parallel finite element element library) for the discretization of the PDE and on PETSc for scalable and efficient linear algebra operations and solvers.

For building instructions, see the file INSTALL.md. Copyright information and licensing restrictions can be found in the file COPYRIGHT.

The best starting point for new users interested in hIPPYlib’s features are the interactive tutorials in the tutorial folder.

Conceptually, hIPPYlib can be viewed as a toolbox that provides the building blocks for experimenting new ideas and developing scalable algorithms for PDE-based deterministic and Bayesian inverse problems.

In hIPPYlib the user can express the forward PDE and the likelihood in weak form using the friendly, compact, near-mathematical notation of FEniCS, which will then automatically generate efficient code for the discretization. Linear and nonlinear, and stationary and time-dependent PDEs are supported in hIPPYlib. Currently, gradient and Hessian information needs to be provided by the user; our future plan includes automatic generation of this information by using FEniCS automatic differentiation capabilities for users who do not wish to derive the relevant weak forms.

Noise and prior covariance operators are modeled as inverses of elliptic differential operators allowing us to build on existing fast multigrid solvers for elliptic operators without explicitly constructing the dense covariance operator.

hIPPYlib provides a robust implementation of the inexact Newton-conjugate gradient algorithm to compute the maximum a posterior (MAP) point. The reduced gradient and Hessian actions are automatically computed via their weak form specification in FEniCS by constraining the state and adjoint variables to satisfy the forward and adjoint problem. The Newton system is solved inexactly by early termination of CG iterations via Eisenstat-Walker (to prevent oversolving) and Steihaug (to avoid negative curvature) criteria. Globalization is achieved with an Armijo back-tracking line search.

hIPPYlib offers different scalable methods to sample from the prior distribution: Krylov methods to approximate the action of a matrix square root on a vector, the conjugate gradient sampler, and finally samplers that exploit the finite element assembly procedure of the inverse covariance operator to obtain a symmetric decomposition. To sample from a local Gaussian approximation (such as at the MAP point) hIPPYlib exploits the low rank factorization of the Hessian misfit to correct samples from the prior distribution.

Finally, randomized and probing algorithms are available to compute the variance of the prior/posterior distribution.

Installation

Inverse Problem PYthon library

 __        ______  _______   _______   __      __  __  __  __
/  |      /      |/       /       /     /  |/  |/  |/  |
$$ |____  $$$$$$/ $$$$$$$  |$$$$$$$  |$$   /$$/ $$ |$$/ $$ |____
$$        $$ |  $$ |__$$ |$$ |__$$ | $$  \/$$/  $$ |/  |$$
$$$$$$$  |  $$ |  $$    $$/ $$    $$/   $$  $$/   $$ |$$ |$$$$$$$  |
$$ |  $$ |  $$ |  $$$$$$$/  $$$$$$$/     $$$$/    $$ |$$ |$$ |  $$ |
$$ |  $$ | _$$ |_ $$ |      $$ |          $$ |    $$ |$$ |$$ |__$$ |
$$ |  $$ |/ $$   |$$ |      $$ |          $$ |    $$ |$$ |$$    $$/
$$/   $$/ $$$$$$/ $$/       $$/           $$/     $$/ $$/ $$$$$$$/

https://hippylib.github.io

hIPPYlib depends on FEniCS version 1.6 or above. The suggested version of FEniCS to use with hIPPYlib is 2017.2.

FEniCS needs to be built with the following dependecies:

• numpy, scipy, matplotlib
• PETSc and petsc4py (version 3.7.0 or above)
• SLEPc and slepc4py (version 3.7.0 or above)
• PETSc dependencies: parmetis, scotch, suitesparse, superlu_dist, ml
• (optional): mshr, jupyter

Run FEniCS from Docker (Linux, MacOS, Windows)

An easy way to run FEniCS is to use their prebuilt Docker images.

First you will need to install Docker on your system. MacOS and Windows users should preferably use Docker for Mac or Docker for Windows — if it is compatible with their system — instead of the legacy version Docker Toolbox.

Among the many docker’s workflow discussed here, we suggest using the Jupyter notebookone.

Docker for Mac, Docker for Windows and Linux users (Setup and first use instructions)

We first create a new Docker container to run the jupyter-notebook command and to expose port 8888. From a command line shell, go to the hippylib folder and type:

docker run --name hippylib-nb -w /home/fenics/hippylib -v $(pwd):/home/fenics/hippylib -d -p 127.0.0.1:8888:8888 quay.io/fenicsproject/stable:2017.2.0.r4 'jupyter-notebook --ip=0.0.0.0'
docker logs hippylib-nb

The notebook will be available at http://localhost:8888/?token=<security_token_for_first_time_connection> in your web browser. From there you can run the interactive notebooks or create a new shell (directly from your browser) to run python scripts.

Docker Toolbox users on Mac/Windows (Setup and first use instructions)

Docker Toolbox is for older Mac and Windows systems that do not meet the requirements of Docker for Mac or Docker for Windows. Docker Toolbox will first create a lightweight linux virtual machine on your system and run docker from the virtual machine. This has implications on the workflow presented above.

We first create a new Docker container to run the jupyter-notebook command and to expose port 8888 on the virtual machine. From a command line shell, go to the hippylib folder and type:

docker run --name hippylib-nb -w /home/fenics/hippylib -v $(pwd):/home/fenics/hippylib -d -p $(docker-machine ip $(docker-machine active)):8888:8888 quay.io/fenicsproject/stable:2017.2.0.r4 'jupyter-notebook --ip=0.0.0.0'
docker logs hippylib-nb

To find out the IP of the virtual machine we type:

docker-machine ip $(docker-machine active)

The notebook will be available at http://<ip-of-virtual-machine>:8888/?token=<security_token_for_first_time_connection> in your web browser. From there you can run the interactive notebooks or create a new shell (directly from your browser) to run python scripts.

Subsequent uses

The docker container will continue to run in the background until we stop it:

docker stop hippylib-nb

To start it again just run:

docker start hippylib-nb

If you would like to see the log output from the Jupyter notebook server (e.g. if you need the security token) type:

docker logs hippylib-nb

Install FEniCS from Conda (Linux or MacOS)

To use the prebuilt Anaconda Python packages (Linux and Mac only), first install Anaconda, then run following commands in your terminal:

conda create -n fenicsproject -c conda-forge fenics=2017.2.0
source activate fenicsproject

Build FEniCS from source using hashdist (Linux and MacOS 10.12 or below)

To build FEniCS from source we suggest using the scripts and profile files in fenics-hashdist. These scripts and profile files contain small modifications with respect to the ones provided by the FEniCS community to ensure that all the dependencies needed by hIPPYlib are installed.

See fenics-hashdist/README.md for further details.

Other ways to build FEniCS

For instructions on other ways to build FEniCS, we refer to the FEniCS project download page. Note that this instructions always refer to the latest version of FEniCS which may or may not be yet supported by hIPPYlib. Always check the hIPPYlib website for supported FEniCS versions.

Build the hIPPYlib documentation using Sphinx

To build the documentation you need to have sphinx (tested on v.1.7.5), m2r and sphinx_rtd_theme - all of these can be installed via pip.

To build simply run make html from doc folder.

Changelog

Inverse Problem PYthon library

 __        ______  _______   _______   __      __  __  __  __
/  |      /      |/       /       /     /  |/  |/  |/  |
$$ |____  $$$$$$/ $$$$$$$  |$$$$$$$  |$$   /$$/ $$ |$$/ $$ |____
$$        $$ |  $$ |__$$ |$$ |__$$ | $$  \/$$/  $$ |/  |$$
$$$$$$$  |  $$ |  $$    $$/ $$    $$/   $$  $$/   $$ |$$ |$$$$$$$  |
$$ |  $$ |  $$ |  $$$$$$$/  $$$$$$$/     $$$$/    $$ |$$ |$$ |  $$ |
$$ |  $$ | _$$ |_ $$ |      $$ |          $$ |    $$ |$$ |$$ |__$$ |
$$ |  $$ |/ $$   |$$ |      $$ |          $$ |    $$ |$$ |$$    $$/
$$/   $$/ $$$$$$/ $$/       $$/           $$/     $$/ $$/ $$$$$$$/

https://hippylib.github.io

Version 2.0.0, released on June 15, 2018

• Introduce capabilities for non-Gaussian Bayesian inference using Mark Chain Monte Carlo methods. Kernels: mMALA, pCN, gpCN, IS. Note: API subject to change
• Support domain-decomposition parallelization (new parallel random number generator, and new randomized eigensolvers)
• The parameter, usually labeled a, throughout the library, has been renamed to m, for model parameter. Interface changes:
  • PDEProblem.eval_da –> PDEProblem.evalGradientParameter
  • Model.applyWua –> Model.applyWum
  • Model.applyWau –> Model.applyWmu
  • Model.applyRaa –> Model.applyWmm
  • gda_tolerance –> gdm_tolerance in the parameter list for Newton and QuasiNewton optimizers
  • gn_approx –> gass_newton_approx as parameter in function to compute Hessian/linearization point in classes Model, PDEProblem, Misfit, Qoi, ReducedQoi
• Organize hippylib in subpackages
• Add sphinx documentation (thanks to E. Khattatov and I. Ambartsumyan)

Version 1.6.0, released on May 16, 2018

• Bugfix in PDEVariationalProblem.solveIncremental for non self-adjoint models
• Add new estimator for the trace and diagonal of the prior covariance using randomized eigendecomposition
• In all examples and tutorial, use enviromental variable HIPPYLIB_BASE_DIR (if defined) to add hIPPYlib to PYTHONPATH

Version 1.5.0, released on Jan 24, 2018

• Add support for FEniCS 2017.2

Version 1.4.0, released on Nov 8, 2017

• Add support for Python 3
• Enchantments in PDEVariationalProblem: it now supports multiple Dirichlet condition and vectorial/mixed function spaces
• Bugfix: Set the correct number of global rows, when targets points fall outside the computational domain
• More extensive testing with Travis Integration

Version 1.3.0, released on June 28, 2017

• Improve hashdist installation support
• Switch license to GPL-2
• Add support for FEniCS 2017.1

Version 1.2.0, released on April 24, 2017

• Update instruction to build FEniCS: hashdist and docker
• Update notebook to nbformat 4
• Let FEniCS 2016.2 be the preferred version of FEniCS
• Add Travis integration

Version 1.1.0, released on Nov 28, 2016

• Add partial support for FEniCS 2016.1 (Applications and Tutorial)
• Improve performance of the randomized eigensolvers

Version 1.0.2, released on Sep 30, 2016

• Use vector2Function to safely convert dolfin.Vector to dolfin.Function
• Optimize the PDEVariationalProblem to exploit the case when the forward problem is linear
• Update notebook 1_FEniCS101.ipynb

Version 1.0.1, released on Aug 25, 2016

• Add support in hippylib.Model and hippylib.Misfit for misfit functional with explicit dependence on the parameter

Version 1.0.0, released on Aug 8, 2016

• Uploaded to https://hippylib.github.io.
• Initial release under GPL-3.

hippylib

hippylib package

Subpackages

hippylib.modeling package

Submodules

hippylib.modeling.PDEProblem module

class hippylib.modeling.PDEProblem.PDEProblem

Bases: object

Consider the PDE problem: Given \(m\), find \(u\) such that

\[F(u, m, p) = ( f(u, m), p) = 0, \quad \forall p.\]

Here \(F\) is linear in \(p\), but it may be non linear in \(u\) and \(m\).

apply_ij(i, j, dir, out)

Given \(u, m, p\); compute \(\delta_{ij} F(u, m, p; \hat{i}, \tilde{j})\) in the direction \(\tilde{j} =\) dir, \(\forall \hat{i}.\)

apply_ijk(i, j, k, x, jdir, kdir, out)

Given x = [u,a,p]; compute \(\delta_{ijk} F(u,a,p; \hat{i}, \tilde{j}, \tilde{k})\) in the direction \((\tilde{j},\tilde{k}) = (\)jdir,kdir), \(\forall \hat{i}.\)

evalGradientParameter(x, out)

Given \(u, m, p\); evaluate \(\delta_m F(u, m, p; \hat{m}),\, \forall \hat{m}.\)

generate_parameter()

Return a vector in the shape of the parameter.

generate_state()

Return a vector in the shape of the state.

init_parameter(m)

Initialize the parameter.

setLinearizationPoint(x, gauss_newton_approx)

Set the values of the state and parameter for the incremental forward and adjoint solvers. Set whether Gauss Newton approximation of the Hessian should be used.

solveAdj(state, x, adj_rhs, tol)

Solve the linear adjoint problem: Given \(m\), \(u\); find \(p\) such that

\[\delta_u F(u, m, p;\hat{u}) = 0, \quad \forall \hat{u}.\]

solveFwd(state, x, tol)

Solve the possibly nonlinear forward problem: Given \(m\), find \(u\) such that

\[\delta_p F(u, m, p;\hat{p}) = 0, \quad \forall \hat{p}.\]

solveIncremental(out, rhs, is_adj, mytol)

If is_adj = False:

Solve the forward incremental system: Given \(u, m\), find \(\tilde{u}\) such that

\[\delta_{pu} F(u, m, p; \hat{p}, \tilde{u}) = \mbox{rhs}, \quad \forall \hat{p}.\]

If is_adj = True:

Solve the adjoint incremental system: Given \(u, m\), find \(\tilde{p}\) such that

\[\delta_{up} F(u, m, p; \hat{u}, \tilde{p}) = \mbox{rhs}, \quad \forall \hat{u}.\]

class hippylib.modeling.PDEProblem.PDEVariationalProblem(Vh, varf_handler, bc, bc0, is_fwd_linear=False)

Bases: hippylib.modeling.PDEProblem.PDEProblem

_createLUSolver()

apply_ij(i, j, dir, out)

Given \(u, m, p\); compute \(\delta_{ij} F(u, m, p; \hat{i}, \tilde{j})\) in the direction \(\tilde{j} =\) dir, \(\forall \hat{i}\).

apply_ijk(i, j, k, x, jdir, kdir, out)

Given x = [u,a,p]; compute \(\delta_{ijk} F(u,a,p; \hat{i}, \tilde{j}, \tilde{k})\) in the direction \((\tilde{j},\tilde{k}) = (\)jdir,kdir), \(\forall \hat{i}.\)

evalGradientParameter(x, out)

Given \(u, m, p\); evaluate \(\delta_m F(u, m, p; \hat{m}),\, \forall \hat{m}.\)

generate_parameter()

Return a vector in the shape of the parameter.

generate_state()

Return a vector in the shape of the state.

init_parameter(m)

Initialize the parameter.

setLinearizationPoint(x, gauss_newton_approx)

Set the values of the state and parameter for the incremental forward and adjoint solvers.

solveAdj(adj, x, adj_rhs, tol)

Solve the linear adjoint problem: Given \(m, u\); find \(p\) such that

\[\delta_u F(u, m, p;\hat{u}) = 0, \quad \forall \hat{u}.\]

solveFwd(state, x, tol)

Solve the possibly nonlinear forward problem: Given \(m\), find \(u\) such that

\[\delta_p F(u, m, p;\hat{p}) = 0,\quad \forall \hat{p}.\]

solveIncremental(out, rhs, is_adj, mytol)

If is_adj == False:

Solve the forward incremental system: Given \(u, m\), find \(\tilde{u}\) such that

\[\delta_{pu} F(u, m, p; \hat{p}, \tilde{u}) = \mbox{rhs},\quad \forall \hat{p}.\]

If is_adj == True:

Solve the adjoint incremental system: Given \(u, m\), find \(\tilde{p}\) such that

\[\delta_{up} F(u, m, p; \hat{u}, \tilde{p}) = \mbox{rhs},\quad \forall \hat{u}.\]

hippylib.modeling.expression module

hippylib.modeling.misfit module

class hippylib.modeling.misfit.ContinuousStateObservation(Vh, dX, bcs, form=None)

Bases: hippylib.modeling.misfit.Misfit

This class implements continuous state observations in a subdomain \(X \subset \Omega\) or \(X \subset \partial \Omega\).

Constructor:

Vh: the finite element space for the state variable.

dX: the integrator on subdomain X where observation are presents. E.g. dX = dl.dx means observation on all \(\Omega\) and dX = dl.ds means observations on all \(\partial \Omega\).

bcs: If the forward problem imposes Dirichlet boundary conditions \(u = u_D \mbox{ on } \Gamma_D\); bcs is a list of dolfin.DirichletBC object that prescribes homogeneuos Dirichlet conditions \(u = 0 \mbox{ on } \Gamma_D\).

form: if form = None we compute the \(L^2(X)\) misfit: \(\int_X (u - u_d)^2 dX,\) otherwise the integrand specified in the given form will be used.

apply_ij(i, j, dir, out)

Apply the second variation \(\delta_{ij}\) (i,j = STATE,PARAMETER) of the cost in direction dir.

cost(x)

Given x evaluate the cost functional. Only the state u and (possibly) the parameter m are accessed.

grad(i, x, out)

Given the state and the paramter in x, compute the partial gradient of the misfit functional in with respect to the state (i == STATE) or with respect to the parameter (i == PARAMETER).

setLinearizationPoint(x, gauss_newton_approx=False)

Set the point for linearization.

Inputs:

x=[u, m, p] - linearization point

gauss_newton_approx (bool) - whether to use Gauss Newton approximation

class hippylib.modeling.misfit.Misfit

Bases: object

Abstract class to model the misfit componenet of the cost functional. In the following x will denote the variable [u, m, p], denoting respectively the state u, the parameter m, and the adjoint variable p.

The methods in the class misfit will usually access the state u and possibly the parameter m. The adjoint variables will never be accessed.

apply_ij(i, j, dir, out)

Apply the second variation \(\delta_{ij}\) (i,j = STATE,PARAMETER) of the cost in direction dir.

cost(x)

Given x evaluate the cost functional. Only the state u and (possibly) the parameter m are accessed.

grad(i, x, out)

Given the state and the paramter in x, compute the partial gradient of the misfit functional in with respect to the state (i == STATE) or with respect to the parameter (i == PARAMETER).

setLinearizationPoint(x, gauss_newton_approx=False)

Set the point for linearization.

Inputs:

x=[u, m, p] - linearization point

gauss_newton_approx (bool) - whether to use Gauss Newton approximation

class hippylib.modeling.misfit.MultiStateMisfit(misfits)

Bases: hippylib.modeling.misfit.Misfit

append(misfit)

apply_ij(i, j, dir, out)

Apply the second variation \(\delta_{ij}\) (i,j = STATE,PARAMETER) of the cost in direction dir.

cost(x)

Given x evaluate the cost functional. Only the state u and (possibly) the parameter m are accessed.

grad(i, x, out)

Given the state and the paramter in x, compute the partial gradient of the misfit functional in with respect to the state (i == STATE) or with respect to the parameter (i == PARAMETER).

setLinearizationPoint(x, gauss_newton_approx=False)

Set the point for linearization.

Inputs:

x=[u, m, p] - linearization point

gauss_newton_approx (bool) - whether to use Gauss Newton approximation

class hippylib.modeling.misfit.PointwiseStateObservation(Vh, obs_points)

Bases: hippylib.modeling.misfit.Misfit

This class implements pointwise state observations at given locations. It assumes that the state variable is a scalar function.

Constructor:

Vh is the finite element space for the state variable

obs_points is a 2D array number of points by geometric dimensions that stores the location of the observations.

apply_ij(i, j, dir, out)

Apply the second variation \(\delta_{ij}\) (i,j = STATE,PARAMETER) of the cost in direction dir.

cost(x)

Given x evaluate the cost functional. Only the state u and (possibly) the parameter m are accessed.

grad(i, x, out)

Given the state and the paramter in x, compute the partial gradient of the misfit functional in with respect to the state (i == STATE) or with respect to the parameter (i == PARAMETER).

setLinearizationPoint(x, gauss_newton_approx=False)

Set the point for linearization.

Inputs:

x=[u, m, p] - linearization point

gauss_newton_approx (bool) - whether to use Gauss Newton approximation

hippylib.modeling.model module

class hippylib.modeling.model.Model(problem, prior, misfit)

This class contains the full description of the inverse problem. As inputs it takes a PDEProblem object, a Prior object, and a Misfit object.

In the following we will denote with

• u the state variable
• m the (model) parameter variable
• p the adjoint variable

Create a model given:

• problem: the description of the forward/adjoint problem and all the sensitivities
• prior: the prior component of the cost functional
• misfit: the misfit componenent of the cost functional

Rsolver()

Return an object Rsovler that is a suitable solver for the regularization operator \(R\).

The solver object should implement the method Rsolver.solve(z,r) such that \(Rz pprox r\).

applyC(dm, out)

Apply the \(C\) block of the Hessian to a (incremental) parameter variable, i.e. out = \(C dm\)

Parameters:

• dm the (incremental) parameter variable
• out the action of the \(C\) block on dm

Note

This routine assumes that out has the correct shape.

applyCt(dp, out)

Apply the transpose of the \(C\) block of the Hessian to a (incremental) adjoint variable. out = \(C^t dp\)

Parameters:

• dp the (incremental) adjoint variable
• out the action of the \(C^T\) block on dp

..note:: This routine assumes that out has the correct shape.

applyR(dm, out)

Apply the regularization \(R\) to a (incremental) parameter variable. out = \(R dm\)

Parameters:

• dm the (incremental) parameter variable
• out the action of \(R\) on dm

Note

This routine assumes that out has the correct shape.

applyWmm(dm, out)

Apply the \(W_{mm}\) block of the Hessian to a (incremental) parameter variable. out = \(W_{mm} dm\)

Parameters:

• dm the (incremental) parameter variable
• out the action of the \(W_{mm}\) on block dm

Note

This routine assumes that out has the correct shape.

applyWmu(du, out)

Apply the \(W_{mu}\) block of the Hessian to a (incremental) state variable. out = \(W_{mu} du\)

Parameters:

• du the (incremental) state variable
• out the action of the \(W_{mu}\) block on du

Note

This routine assumes that out has the correct shape.

applyWum(dm, out)

Apply the \(W_{um}\) block of the Hessian to a (incremental) parameter variable. out = \(W_{um} dm\)

Parameters:

• dm the (incremental) parameter variable
• out the action of the \(W_{um}\) block on du

Note

This routine assumes that out has the correct shape.

applyWuu(du, out)

Apply the \(W_{uu}\) block of the Hessian to a (incremental) state variable. out = \(W_{uu} du\)

Parameters:

• du the (incremental) state variable
• out the action of the \(W_{uu}\) block on du

Note

This routine assumes that out has the correct shape.

apply_ij(i, j, d, out)

cost(x)

Given the list x = [u,m,p] which describes the state, parameter, and adjoint variable compute the cost functional as the sum of the misfit functional and the regularization functional.

Return the list [cost functional, regularization functional, misfit functional]

Note

p is not needed to compute the cost functional

evalGradientParameter(x, mg, misfit_only=False)

Evaluate the gradient for the variational parameter equation at the point x=[u,m,p].

Parameters:

• x = [u,m,p] the point at which to evaluate the gradient.
• mg the variational gradient \((g, mtest)\), mtest being a test function in the parameter space (Output parameter)

Returns the norm of the gradient in the correct inner product \(g_norm = sqrt(g,g)\)

generate_vector(component='ALL')

By default, return the list [u,m,p] where:

• u is any object that describes the state variable
• m is a dolfin.Vector object that describes the parameter variable. (Needs to support linear algebra operations)
• p is any object that describes the adjoint variable

If component = STATE return only u

If component = PARAMETER return only m

If component = ADJOINT return only p

init_parameter(m)

Reshape m so that it is compatible with the parameter variable

setPointForHessianEvaluations(x, gauss_newton_approx=False)

Specify the point x = [u,m,p] at which the Hessian operator (or the Gauss-Newton approximation) needs to be evaluated.

Parameters:

• x = [u,m,p]: the point at which the Hessian or its Gauss-Newton approximation needs to be evaluated.
• gauss_newton_approx (bool): whether to use Gauss-Newton approximation (default: use Newton)

Note

This routine should either:

• simply store a copy of x and evaluate action of blocks of the Hessian on the fly
• or partially precompute the block of the hessian (if feasible)

solveAdj(out, x, tol=1e-09)

Solve the linear adjoint problem.

Parameters:

• out: is the solution of the adjoint problem (i.e. the adjoint p) (Output parameter)

• x = [u, m, p] provides

  1. the parameter variable m for assembling the adjoint operator
  2. the state variable u for assembling the adjoint right hand side

  Note

  p is not accessed

• tol is the relative tolerance for the solution of the adjoint problem. [Default 1e-9].

solveAdjIncremental(sol, rhs, tol)

Solve the incremental adjoint problem for a given right-hand side

Parameters:

• sol the solution of the incremental adjoint problem (Output)
• rhs the right hand side of the linear system
• tol the relative tolerance for the linear system

solveFwd(out, x, tol=1e-09)

Solve the (possibly non-linear) forward problem.

Parameters:

• out: is the solution of the forward problem (i.e. the state) (Output parameters)

• x = [u,m,p] provides

  1. the parameter variable m for the solution of the forward problem
  2. the initial guess u if the forward problem is non-linear

  Note

  p is not accessed

• tol is the relative tolerance for the solution of the forward problem. [Default 1e-9].

solveFwdIncremental(sol, rhs, tol)

Solve the linearized (incremental) forward problem for a given right-hand side

Parameters:

• sol the solution of the linearized forward problem (Output)
• rhs the right hand side of the linear system
• tol the relative tolerance for the linear system

hippylib.modeling.modelVerify module

hippylib.modeling.modelVerify.modelVerify(model, m0, innerTol, is_quadratic=False, misfit_only=False, verbose=True, eps=None)

Verify the reduced Gradient and the Hessian of a model. It will produce two loglog plots of the finite difference checks for the gradient and for the Hessian. It will also check for symmetry of the Hessian.

hippylib.modeling.modelVerify.modelVerifyPlotErrors(is_quadratic, eps, err_grad, err_H)

hippylib.modeling.pointwiseObservation module

hippylib.modeling.pointwiseObservation.assemblePointwiseObservation(Vh, targets)

Assemble the pointwise observation matrix:

Inputs

• Vh: FEniCS finite element space.
• targets: observation points (numpy array).

hippylib.modeling.pointwiseObservation.exportPointwiseObservation(Vh, B, data, fname, varname='observation')

This function writes a VTK PolyData file to visualize pointwise data.

Inputs:

• B: observation operator.
• mesh: mesh.
• data: dolfin.Vector containing the data.
• fname: filename for the file to export (without extension).
• varname: name of the variable for the .vtk file.

hippylib.modeling.pointwiseObservation.write_vtk(points, data, fname, varname='observation')

This function writes a VTK PolyData file to visualize pointwise data.

Inputs:

• points: locations of the points (numpy array of size equal to number of points times space dimension).
• data: pointwise values (numpy array of size equal to number of points).
• fname: filename for the .vtk file to export.
• varname: name of the variable for the .vtk file.

hippylib.modeling.posterior module

class hippylib.modeling.posterior.GaussianLRPosterior(prior, d, U, mean=None)

Class for the low rank Gaussian Approximation of the Posterior. This class provides functionality for approximate Hessian apply, solve, and Gaussian sampling based on the low rank factorization of the Hessian.

In particular if \(d\) and \(U\) are the dominant eigenpairs of \(H_{\mbox{misfit}} U[:,i] = d[i] R U[:,i]\) then we have:

• low rank Hessian apply: \(y = ( R + RU D U^{T}) x.\)
• low rank Hessian solve: \(y = (R^-1 - U (I + D^{-1})^{-1} U^T) x.\)
• low rank Hessian Gaussian sampling: \(y = ( I - U S U^{T}) x\), where \(S = I - (I + D)^{-1/2}\) and \(x \sim \mathcal{N}(0, R^{-1}).\)

Construct the Gaussian approximation of the posterior. Input: - prior: the prior mode. - d: the dominant generalized eigenvalues of the Hessian misfit. - U: the dominant generalized eigenvector of the Hessian misfit \(U^T R U = I.\) - mean: the MAP point.

_sample_given_prior(s_prior, s_post)

_sample_given_white_noise(noise, s_prior, s_post)

cost(m)

init_vector(x, dim)

Inizialize a vector x to be compatible with the range/domain of \(H\). If dim == "noise" inizialize x to be compatible with the size of white noise used for sampling.

klDistanceFromPrior(sub_comp=False)

pointwise_variance(**kwargs)

Compute/estimate the pointwise variance of the posterior, prior distribution and the pointwise variance reduction informed by the data.

See _Prior.pointwise_variance for more details.

sample(*args, **kwargs)

possible calls:

1. sample(s_prior, s_post, add_mean=True)

   Given a prior sample s_prior compute a sample s_post from the posterior.

   • s_prior is a sample from the prior centered at 0 (input).
   • s_post is a sample from the posterior (output).
   • if add_mean=True (default) then the samples will be centered at the map point.

2. sample(noise, s_prior, s_post, add_mean=True)

   Given noise \(\sim \mathcal{N}(0, I)\) compute a sample s_prior from the prior and s_post from the posterior.

   • noise is a realization of white noise (input).
   • s_prior is a sample from the prior (output).
   • s_post is a sample from the posterior.
   • if add_mean=True (default) then the prior and posterior samples will be centered at the respective means.

trace(**kwargs)

Compute/estimate the trace of the posterior, prior distribution and the trace of the data informed correction.

See _Prior.trace for more details.

trace_update()

class hippylib.modeling.posterior.LowRankHessian(prior, d, U)

Operator that represents the action of the low rank approximation of the Hessian and of its inverse.

init_vector(x, dim)

inner(x, y)

mult(x, y)

solve(sol, rhs)

class hippylib.modeling.posterior.LowRankHessianMisfit(prior, d, U)

Operator that represents the action of the low rank approximation of the Hessian misfit and of its inverse.

init_vector(x, dim)

inner(x, y)

mult(x, y)

class hippylib.modeling.posterior.LowRankPosteriorSampler(prior, d, U)

Object to sample from the low rank approximation of the posterior.

\[y = ( I - U S U^{T}) x,\]

where

\(S = I - (I + D)^{-1/2}, x \sim \mathcal{N}(0, R^{-1}).\)

init_vector(x, dim)

sample(noise, s)

hippylib.modeling.prior module

class hippylib.modeling.prior.BiLaplacianPrior(Vh, gamma, delta, Theta=None, mean=None, rel_tol=1e-12, max_iter=1000)

Bases: hippylib.modeling.prior._Prior

This class implement a prior model with covariance matrix \(C = (\delta I + \gamma \mbox{div } \Theta \nabla) ^ {-2}\).

The magnitude of \(\delta\gamma\) governs the variance of the samples, while the ratio \(\frac{\gamma}{\delta}\) governs the correlation lenght.

Here \(\Theta\) is a SPD tensor that models anisotropy in the covariance kernel.

Construct the prior model. Input:

• Vh: the finite element space for the parameter
• gamma and delta: the coefficient in the PDE
• Theta: the SPD tensor for anisotropic diffusion of the PDE
• mean: the prior mean

init_vector(x, dim)

Inizialize a vector x to be compatible with the range/domain of \(R\).

If dim == "noise" inizialize x to be compatible with the size of white noise used for sampling.

sample(noise, s, add_mean=True)

Given noise \(\sim \mathcal{N}(0, I)\) compute a sample s from the prior.

If add_mean == True add the prior mean value to s.

class hippylib.modeling.prior.LaplacianPrior(Vh, gamma, delta, mean=None, rel_tol=1e-12, max_iter=100)

Bases: hippylib.modeling.prior._Prior

This class implements a prior model with covariance matrix \(C = (\delta I - \gamma \Delta) ^ {-1}\).

The magnitude of \(\gamma\) governs the variance of the samples, while the ratio \(\frac{\gamma}{\delta}\) governs the correlation length.

Note

\(C\) is a trace class operator only in 1D while it is not a valid prior in 2D and 3D.

Construct the prior model. Input:

• Vh: the finite element space for the parameter
• gamma and delta: the coefficient in the PDE
• Theta: the SPD tensor for anisotropic diffusion of the PDE
• mean: the prior mean

init_vector(x, dim)

Inizialize a vector x to be compatible with the range/domain of \(R\).

If dim == "noise" inizialize x to be compatible with the size of white noise used for sampling.

sample(noise, s, add_mean=True)

Given noise \(\sim \mathcal{N}(0, I)\) compute a sample s from the prior.

If add_mean == True add the prior mean value to s.

class hippylib.modeling.prior.MollifiedBiLaplacianPrior(Vh, gamma, delta, locations, m_true, Theta=None, pen=10.0, order=2, rel_tol=1e-12, max_iter=1000)

Bases: hippylib.modeling.prior._Prior

This class implement a prior model with covariance matrix \(C = \left( [\delta + \mbox{pen} \sum_i m(x - x_i) ] I + \gamma \mbox{div } \Theta \nabla\right) ^ {-2}\),

where

• \(\Theta\) is a SPD tensor that models anisotropy in the covariance kernel.
• \(x_i (i=1,...,n)\) are points were we assume to know exactly the value of the parameter (i.e., \(m(x_i) = m_{\mbox{true}}( x_i) \mbox{ for } i=1,...,n).\)
• \(m\) is the mollifier function: \(m(x - x_i) = \exp\left( - \left[\frac{\gamma}{\delta}\| x - x_i \|_{\Theta^{-1}}\right]^{\mbox{order}} \right).\)
• pen is a penalization parameter.

The magnitude of \(\delta \gamma\) governs the variance of the samples, while the ratio \(\frac{\gamma}{\delta}\) governs the correlation length.

The prior mean is computed by solving

\[\left( [\delta + \sum_i m(x - x_i) ] I + \gamma \mbox{div } \Theta \nabla \right) m = \sum_i m(x - x_i) m_{\mbox{true}}.\]

Construct the prior model. Input:

• Vh: the finite element space for the parameter
• gamma and delta: the coefficients in the PDE
• locations: the points \(x_i\) at which we assume to know the true value of the parameter
• m_true: the true model
• Theta: the SPD tensor for anisotropic diffusion of the PDE
• pen: a penalization parameter for the mollifier

init_vector(x, dim)

Inizialize a vector x to be compatible with the range/domain of \(R\).

If dim == "noise" inizialize x to be compatible with the size of white noise used for sampling.

sample(noise, s, add_mean=True)

Given noise \(\sim \mathcal{N}(0, I)\) compute a sample s from the prior.

If add_mean == True add the prior mean value to s.

class hippylib.modeling.prior._BilaplacianR(A, Msolver)

Operator that represent the action of the regularization/precision matrix for the Bilaplacian prior.

init_vector(x, dim)

inner(x, y)

mpi_comm()

mult(x, y)

class hippylib.modeling.prior._BilaplacianRsolver(Asolver, M)

Operator that represent the action of the inverse the regularization/precision matrix for the Bilaplacian prior.

init_vector(x, dim)

solve(x, b)

class hippylib.modeling.prior._Prior

Abstract class to describe the prior model. Concrete instances of a _Prior class should expose the following attributes and methods.

Attributes:

• R: an operator to apply the regularization/precision operator.
• Rsolver: an operator to apply the inverse of the regularization/precision operator.
• M: the mass matrix in the control space.
• mean: the prior mean.

Methods:

• init_vector(self,x,dim): Inizialize a vector x to be compatible with the range/domain of R If dim == "noise" inizialize x to be compatible with the size of white noise used for sampling.
• sample(self, noise, s, add_mean=True): Given noise \(\sim \mathcal{N}(0, I)\) compute a sample s from the prior. If add_mean==True add the prior mean value to s.

cost(m)

grad(m, out)

init_vector(x, dim)

pointwise_variance(method, k=1000000, r=200)

Compute/estimate the prior pointwise variance.

• If method=="Exact" we compute the diagonal entries of \(R^{-1}\) entry by entry. This requires to solve \(n\) linear system in \(R\) (not scalable, but ok for illustration purposes).

sample(noise, s, add_mean=True)

trace(method='Exact', tol=0.1, min_iter=20, max_iter=100, r=200)

Compute/estimate the trace of the prior covariance operator.

• If method=="Exact" we compute the trace exactly by summing the diagonal entries of \(R^{-1}M\). This requires to solve \(n\) linear system in \(R\) (not scalable, but ok for illustration purposes).
• If method=="Estimator" use the trace estimator algorithms implemeted in the class TraceEstimator. tol is a relative bound on the estimator standard deviation. In particular, we used enough samples in the Estimator such that the standard deviation of the estimator is less then tol\(tr(\mbox{Prior})\). min_iter and max_iter are the lower and upper bound on the number of samples to be used for the estimation of the trace.

class hippylib.modeling.prior._RinvM(Rsolver, M)

Operator that models the action of \(R^{-1}M\). It is used in the randomized trace estimator.

init_vector(x, dim)

mult(x, y)

hippylib.modeling.reducedHessian module

class hippylib.modeling.reducedHessian.FDHessian(model, m0, h, innerTol, misfit_only=False)

This class implements matrix free application of the reduced Hessian operator. The constructor takes the following parameters:

• model: the object which contains the description of the problem.
• m0: the value of the parameter at which the Hessian needs to be evaluated.
• h: the mesh size for FD.
• innerTol: the relative tolerance for the solution of the forward and adjoint problems.
• misfit_only: a boolean flag that describes whenever the full Hessian or only the misfit component of the Hessian is used.

Type help(Template) for more information on which methods model should implement.

Construct the reduced Hessian Operator

init_vector(x, dim)

Reshape the Vector x so that it is compatible with the reduced Hessian operator.

Parameters:

• x: the vector to reshape
• dim: if 0 then x will be reshaped to be compatible with the range of the reduced Hessian, if 1 then x will be reshaped to be compatible with the domain of the reduced Hessian

Note

Since the reduced Hessian is a self adjoint operator, the range and the domain is the same. Either way, we choosed to add the parameter dim for consistency with the interface of Matrix in dolfin.

inner(x, y)

Perform the inner product between x and y in the norm induced by the reduced Hessian \(H,\, (x, y)_H = x' H y\).

mult(x, y)

Apply the reduced Hessian (or the Gauss-Newton approximation) to the vector x. Return the result in y.

class hippylib.modeling.reducedHessian.ReducedHessian(model, innerTol, misfit_only=False)

This class implements matrix free application of the reduced Hessian operator. The constructor takes the following parameters:

• model: the object which contains the description of the problem.
• innerTol: the relative tolerance for the solution of the incremental forward and adjoint problems.
• misfit_only: a boolean flag that describes whenever the full Hessian or only the misfit component of the Hessian is used.

Type help(modelTemplate) for more information on which methods model should implement.

Construct the reduced Hessian Operator

GNHessian(x, y)

Apply the Gauss-Newton approximation of the reduced Hessian to the vector x. Return the result in y.

TrueHessian(x, y)

Apply the the reduced Hessian to the vector x. Return the result in y.

init_vector(x, dim)

Reshape the Vector x so that it is compatible with the reduced Hessian operator.

Parameters:

• x: the vector to reshape.
• dim: if 0 then x will be reshaped to be compatible with the range of the reduced Hessian, if 1 then x will be reshaped to be compatible with the domain of the reduced Hessian.

Note

Since the reduced Hessian is a self adjoint operator, the range and the domain is the same. Either way, we choosed to add the parameter dim for consistency with the interface of Matrix in dolfin.

inner(x, y)

Perform the inner product between x and y in the norm induced by the reduced Hessian \(H,\,(x, y)_H = x' H y\).

mult(x, y)

Apply the reduced Hessian (or the Gauss-Newton approximation) to the vector x. Return the result in y.

hippylib.modeling.timeDependentVector module

class hippylib.modeling.timeDependentVector.TimeDependentVector(times, tol=1e-10, mpi_comm=<sphinx.ext.autodoc.importer._MockObject object>)

Bases: object

A class to store time dependent vectors. Snapshots are stored/retrieved by specifying the time of the snapshot. Times at which the snapshot are taken must be specified in the constructor.

Constructor:

• times: time frame at which snapshots are stored.
• tol : tolerance to identify the frame of the snapshot.

axpy(a, other)

Compute \(x = x + \mbox{a*other}\) snapshot per snapshot.

copy(other)

Copy all the time frames and snapshot from other to self.

initialize(M, dim)

Initialize all the snapshot to be compatible with the range/domain of an operator M.

inner(other)

Compute the inner products: \(a+= (\mbox{self[i]},\mbox{other[i]})\) for each snapshot.

norm(time_norm, space_norm)

Compute the space-time norm of the snapshot.

retrieve(u, t)

Retrieve snapshot u relative to time t. If t does not belong to the list of time frame an error is raised.

store(u, t)

Store snapshot u relative to time t. If t does not belong to the list of time frame an error is raised.

zero()

Zero out each snapshot.

hippylib.modeling.variables module

Module contents

hippylib.algorithms package

Submodules

hippylib.algorithms.NewtonCG module

hippylib.algorithms.NewtonCG.LS_ParameterList()

Generate a ParameterList for line search globalization. type: LS_ParameterList().showMe() for default values and their descriptions

class hippylib.algorithms.NewtonCG.ReducedSpaceNewtonCG(model, parameters=<hippylib.utils.parameterList.ParameterList object>)

Inexact Newton-CG method to solve constrained optimization problems in the reduced parameter space. The Newton system is solved inexactly by early termination of CG iterations via Eisenstat-Walker (to prevent oversolving) and Steihaug (to avoid negative curvature) criteria. Globalization is performed using one of the following methods:

• line search (LS) based on the armijo sufficient reduction condition; or
• trust region (TR) based on the prior preconditioned norm of the update direction.

The stopping criterion is based on a control on the norm of the gradient and a control of the inner product between the gradient and the Newton direction.

The user must provide a model that describes the forward problem, cost functionals, and all the derivatives for the gradient and the Hessian.

More specifically the model object should implement following methods:

• generate_vector() -> generate the object containing state, parameter, adjoint
• cost(x) -> evaluate the cost functional, report regularization part and misfit separately
• solveFwd(out, x,tol) -> solve the possibly non linear forward problem up a tolerance tol
• solveAdj(out, x,tol) -> solve the linear adjoint problem
• evalGradientParameter(x, out) -> evaluate the gradient of the parameter and compute its norm
• setPointForHessianEvaluations(x) -> set the state to perform hessian evaluations
• solveFwdIncremental(out, rhs, tol) -> solve the linearized forward problem for a given rhs
• solveAdjIncremental(out, rhs, tol) -> solve the linear adjoint problem for a given rhs
• applyC(dm, out) –> Compute out \(= C_x dm\)
• applyCt(dp, out) –> Compute out = \(C_x dp\)
• applyWuu(du,out) –> Compute out = \((W_{uu})_x du\)
• applyWmu(dm, out) –> Compute out = \((W_{um})_x dm\)
• applyWmu(du, out) –> Compute out = \(W_{mu} du\)
• applyR(dm, out) –> Compute out = \(R dm\)
• applyWmm(dm,out) –> Compute out = \(W_{mm} dm\)
• Rsolver() –> A solver for the regularization term

Type help(Model) for additional information

Initialize the ReducedSpaceNewtonCG. Type ReducedSpaceNewtonCG_ParameterList().showMe() for list of default parameters and their descriptions.

_solve_ls(x)

Solve the constrained optimization problem with initial guess x.

_solve_tr(x)

solve(x)

Input:

x = [u, m, p] represents the initial guess (u and p may be None). x will be overwritten on return.

termination_reasons = ['Maximum number of Iteration reached', 'Norm of the gradient less than tolerance', 'Maximum number of backtracking reached', 'Norm of (g, dm) less than tolerance']

hippylib.algorithms.NewtonCG.ReducedSpaceNewtonCG_ParameterList()

Generate a ParameterList for ReducedSpaceNewtonCG. type: ReducedSpaceNewtonCG_ParameterList().showMe() for default values and their descriptions

hippylib.algorithms.NewtonCG.TR_ParameterList()

Generate a ParameterList for Trust Region globalization. type: RT_ParameterList().showMe() for default values and their descriptions

hippylib.algorithms.bfgs module

class hippylib.algorithms.bfgs.BFGS(model, parameters=<hippylib.utils.parameterList.ParameterList object>)

Implement BFGS technique with backtracking inexact line search and damped updating See Nocedal & Wright (06), ch.6.2, ch.7.3, ch.18.3

The user must provide a model that describes the forward problem, cost functionals, and all the derivatives for the gradient and the Hessian.

More specifically the model object should implement following methods:

• generate_vector() -> generate the object containing state, parameter, adjoint
• cost(x) -> evaluate the cost functional, report regularization part and misfit separately
• solveFwd(out, x,tol) -> solve the possibly non-linear forward problem up a tolerance tol
• solveAdj(out, x,tol) -> solve the linear adjoint problem
• evalGradientParameter(x, out) -> evaluate the gradient of the parameter and compute its norm
• applyR(dm, out) –> Compute out = \(R dm\)
• Rsolver() –> A solver for the regularization term

Type help(Model) for additional information

Initialize the BFGS solver. Type BFGS_ParameterList().showMe() for default parameters and their description

solve(x, H0inv, bounds_xPARAM=None)

Solve the constrained optimization problem with initial guess x = [u, m0, p].

Note

u and p may be None.

x will be overwritten.

H0inv: the initial approximated inverse of the Hessian for the BFGS operator. It has an optional method update(x) that will update the operator based on x = [u,m,p].

bounds_xPARAM: Bound constraint (list with two entries: min and max). Can be either a scalar value or a dolfin.Vector.

Return the solution [u,m,p]

termination_reasons = ['Maximum number of Iteration reached', 'Norm of the gradient less than tolerance', 'Maximum number of backtracking reached', 'Norm of (g, da) less than tolerance']

hippylib.algorithms.bfgs.BFGS_ParameterList()

class hippylib.algorithms.bfgs.BFGS_operator(parameters=<hippylib.utils.parameterList.ParameterList object>)

set_H0inv(H0inv)

Set user-defined operator corresponding to H0inv

Input:

H0inv: Fenics operator with method solve()

solve(x, b)

Solve system: \(H_{bfgs} x = b\) where \(H_{bfgs}\) is the approximation to the Hessian build by BFGS. That is, we apply

\[x = (H_{bfgs})^{-1} b = H_k b\]

where \(H_k\) matrix is BFGS approximation to the inverse of the Hessian. Computation done via double-loop algorithm.

Inputs:

x = dolfin.Vector - [out]

b = dolfin.Vector - [in]

update(s, y)

Update BFGS operator with most recent gradient update.

To handle potential break from secant condition, update done via damping

Inputs:

s = dolfin.Vector [in] - corresponds to update in medium parameters.

y = dolfin.Vector [in] - corresponds to update in gradient.

hippylib.algorithms.bfgs.BFGSoperator_ParameterList()

class hippylib.algorithms.bfgs.RescaledIdentity(init_vector=None)

Bases: object

Default operator for H0inv, corresponds to applying \(d0 I\)

init_vector(x, dim)

solve(x, b)

hippylib.algorithms.cgsampler module

class hippylib.algorithms.cgsampler.CGSampler

This class implements the CG sampler algorithm to generate samples from \(\mathcal{N}(0, A^{-1})\).

Reference:

Albert Parker and Colin Fox Sampling Gaussian Distributions in Krylov Spaces with Conjugate Gradient SIAM J SCI COMPUT, Vol 34, No. 3 pp. B312-B334

Construct the solver with default parameters tolerance = 1e-4

print_level = 0

verbose = 0

sample(noise, s)

Generate a sample \(s ~ N(0, A^{-1})\).

noise is a numpy.array of i.i.d. normal variables used as input. For a fixed realization of noise the algorithm is fully deterministic. The size of noise determine the maximum number of CG iterations.

set_operator(A)

Set the operator A, such that \(x \sim \mathcal{N}(0, A^{-1})\).

Note

A is any object that provides the methods init_vector() and mult()

hippylib.algorithms.cgsolverSteihaug module

class hippylib.algorithms.cgsolverSteihaug.CGSolverSteihaug(parameters=<hippylib.utils.parameterList.ParameterList object>, comm=<sphinx.ext.autodoc.importer._MockObject object>)

Solve the linear system \(A x = b\) using preconditioned conjugate gradient ( \(B\) preconditioner) and the Steihaug stopping criterion:

• reason of termination 0: we reached the maximum number of iterations (no convergence)
• reason of termination 1: we reduced the residual up to the given tolerance (convergence)
• reason of termination 2: we reached a negative direction (premature termination due to not spd matrix)
• reason of termination 3: we reached the boundary of the trust region

The stopping criterion is based on either

• the absolute preconditioned residual norm check: \(|| r^* ||_{B^{-1}} < atol\)
• the relative preconditioned residual norm check: \(|| r^* ||_{B^{-1}}/|| r^0 ||_{B^{-1}} < rtol,\)

where \(r^* = b - Ax^*\) is the residual at convergence and \(r^0 = b - Ax^0\) is the initial residual.

The operator A is set using the method set_operator(A). A must provide the following two methods:

• A.mult(x,y): y = Ax
• A.init_vector(x, dim): initialize the vector x so that it is compatible with the range (dim = 0) or the domain (dim = 1) of A.

The preconditioner B is set using the method set_preconditioner(B). B must provide the following method: - B.solve(z,r): z is the action of the preconditioner B on the vector r

To solve the linear system \(Ax = b\) call self.solve(x,b). Here x and b are assumed to be dolfin.Vector objects.

Type CGSolverSteihaug_ParameterList().showMe() for default parameters and their descriptions

reason = ['Maximum Number of Iterations Reached', 'Relative/Absolute residual less than tol', 'Reached a negative direction', 'Reached trust region boundary']

set_TR(radius, B_op)

set_operator(A)

Set the operator \(A\).

set_preconditioner(B_solver)

Set the preconditioner \(B\).

solve(x, b)

Solve the linear system \(Ax = b\)

update_x_with_TR(x, alpha, d)

update_x_without_TR(x, alpha, d)

hippylib.algorithms.cgsolverSteihaug.CGSolverSteihaug_ParameterList()

Generate a ParameterList for CGSolverSteihaug. Type CGSolverSteihaug_ParameterList().showMe() for default values and their descriptions

hippylib.algorithms.linalg module

class hippylib.algorithms.linalg.DiagonalOperator(d)

init_vector(x, dim)

inner(x, y)

mult(x, y)

hippylib.algorithms.linalg.GetFromOwnedGid(v, gid)

hippylib.algorithms.linalg.MatAtB(A, B)

Compute the matrix-matrix product \(A^T B\).

hippylib.algorithms.linalg.MatMatMult(A, B)

Compute the matrix-matrix product \(AB\).

hippylib.algorithms.linalg.MatPtAP(A, P)

Compute the triple matrix product \(P^T A P\).

class hippylib.algorithms.linalg.Operator2Solver(op, mpi_comm=<sphinx.ext.autodoc.importer._MockObject object>)

init_vector(x, dim)

inner(x, y)

solve(y, x)

hippylib.algorithms.linalg.SetToOwnedGid(v, gid, val)

class hippylib.algorithms.linalg.Solver2Operator(S, mpi_comm=<sphinx.ext.autodoc.importer._MockObject object>)

init_vector(x, dim)

inner(x, y)

mult(x, y)

hippylib.algorithms.linalg.Transpose(A)

Compute the matrix transpose

hippylib.algorithms.linalg.amg_method()

Determine which AMG preconditioner to use. If avaialable use ML, which is faster than the PETSc one.

hippylib.algorithms.linalg.estimate_diagonal_inv2(Asolver, k, d)

An unbiased stochastic estimator for the diagonal of \(A^{-1}\). \(d = [ \sum_{j=1}^k v_j .* A^{-1} v_j ] ./ [ \sum_{j=1}^k v_j .* v_j ]\) where

• \(v_j\) are i.i.d. \(\mathcal{N}(0, I)\)
• \(.*\) and \(./\) represent the element-wise multiplication and division of vectors, respectively.

Reference:

Costas Bekas, Effrosyni Kokiopoulou, and Yousef Saad, An estimator for the diagonal of a matrix, Applied Numerical Mathematics, 57 (2007), pp. 1214-1229.

hippylib.algorithms.linalg.get_diagonal(A, d)

Compute the diagonal of the square operator \(A\). Use Solver2Operator if \(A^{-1}\) is needed.

hippylib.algorithms.linalg.to_dense(A, mpi_comm=<sphinx.ext.autodoc.importer._MockObject object>)

Convert a sparse matrix A to dense. For debugging only.

hippylib.algorithms.linalg.trace(A, mpi_comm=<sphinx.ext.autodoc.importer._MockObject object>)

Compute the trace of a sparse matrix \(A\).

hippylib.algorithms.lowRankOperator module

class hippylib.algorithms.lowRankOperator.LowRankOperator(d, U, my_init_vector=None)

This class model the action of a low rank operator \(A = U D U^T\). Here \(D\) is a diagonal matrix, and the columns of are orthonormal in some weighted inner-product.

Note

This class only works in serial!

Construct the low rank operator given d and U.

get_diagonal(diag)

Compute the diagonal of A.

init_vector(x, dim)

Initialize x to be compatible with the range (dim=0) or domain (dim=1) of A.

inner(x, y)

mult(x, y)

Compute \(y = Ax = U D U^T x\)

solve(sol, rhs)

Compute \(\mbox{sol} = U D^-1 U^T x\)

trace(W=None)

Compute the trace of A. If the weight W is given, compute the trace of \(W^{1/2} A W^{1/2}\). This is equivalent to \(\mbox{tr}_W(A) = \sum_i \lambda_i\), where \(\lambda_i\) are the generalized eigenvalues of \(A x = \lambda W^{-1} x\).

Note

If \(U\) is a \(W\)-orthogonal matrix then \(\mbox{tr}_W(A) = \sum_i D(i,i)\).

trace2(W=None)

Compute the trace of \(A A\) (Note this is the square of Frobenius norm, since \(A\) is symmetic). If the weight W is provided, it will compute the trace of \((AW)^2\).

This is equivalent to \(\mbox{tr}_W(A) = \sum_i \lambda_i^2\), where \(\lambda_i\) are the generalized eigenvalues of \(A x = \lambda W^{-1} x\).

Note

If \(U\) is a \(W\)-orthogonal matrix then \(\mbox{tr}_W(A) = \sum_i D(i,i)^2\).

hippylib.algorithms.multivector module

hippylib.algorithms.multivector.MatMvMult(A, x, y)

class hippylib.algorithms.multivector.MultiVector(*args, **kwargs)

Bases: sphinx.ext.autodoc.importer._MockObject

Borthogonalize(B)

Returns \(QR\) decomposition of self. \(Q\) and \(R\) satisfy the following relations in exact arithmetic

\[ \begin{align}\begin{aligned}R \,= \,Z, && (1),\\Q^*BQ\, = \, I, && (2),\\Q^*BZ \, = \,R, && (3),\\ZR^{-1} \, = \, Q, && (4). \end{aligned}\end{align} \]

Returns:

Bq of type MultiVector -> B\(^{-1}\)-orthogonal vectors r of type ndarray -> The \(r\) of the QR decomposition.

Note

self is overwritten by \(Q\).

_mgs_reortho()

_mgs_stable(B)

Returns \(QR\) decomposition of self, which satisfies conditions (1)–(4). Uses Modified Gram-Schmidt with re-orthogonalization (Rutishauser variant) for computing the \(B\)-orthogonal \(QR\) factorization.

References:

1. A.K. Saibaba, J. Lee and P.K. Kitanidis, Randomized algorithms for Generalized Hermitian Eigenvalue Problems with application to computing Karhunen-Loe’ve expansion http://arxiv.org/abs/1307.6885
2. W. Gander, Algorithms for the QR decomposition. Res. Rep, 80(02), 1980

https://github.com/arvindks/kle

dot_mv(mv)

dot_v(v)

export(Vh, filename, varname='mv', normalize=False)

Export in paraview this multivector.

Inputs:

• Vh: the parameter finite element space.
• filename: the name of the paraview output file.
• varname: the name of the paraview variable.
• normalize: if True the vector is rescaled such that \(\| u \|_{\infty} = 1.\)

norm(norm_type)

orthogonalize()

Returns \(QR\) decomposition of self. \(Q\) and \(R\) satisfy the following relations in exact arithmetic

\[ \begin{align}\begin{aligned}QR \, = \, Z, && (1),\\Q^*Q \, = \, I, && (2),\\Q^*Z \, = \, R, && (3),\\ZR^{-1} \, = \, Q, && (4).\end{aligned}\end{align} \]

Returns:

r of type ndarray -> The r of the QR decomposition

Note

self is overwritten by \(Q\).

hippylib.algorithms.multivector.MvDSmatMult(X, A, Y)

hippylib.algorithms.randomizedEigensolver module

hippylib.algorithms.randomizedEigensolver.check_g(A, B, U, d)

Test the frobenious norm of \(U^TBU - I_k\).

Test the frobenious norm of \((V^TAV) - I_k\), with \(V = U D^{-1/2}\).

Test the \(l_2\) norm of the residual: \(r[i] = A U[i] - d[i] B U[i]\).

hippylib.algorithms.randomizedEigensolver.check_std(A, U, d)

Test the frobenious norm of \(U^TU - I_k\).

Test the frobenious norm of \((V^TAV) - I_k\), with \(V = U D^{-1/2}\).

Test the \(l_2\) norm of the residual: \(r[i] = A U[i] - d[i] U[i]\).

hippylib.algorithms.randomizedEigensolver.doublePass(A, Omega, k, s, check=False)

The double pass algorithm for the HEP as presented in [1].

Inputs:

• A: the operator for which we need to estimate the dominant eigenpairs.
• Omega: a random gassian matrix with \(m \geq k\) columns.
• k: the number of eigenpairs to extract.
• s: the number of power iterations for selecting the subspace.

Outputs:

• d: the estimate of the \(k\) dominant eigenvalues of \(A\).
• U: the estimate of the \(k\) dominant eigenvectors of \(A,\, U^T U = I_k\).

hippylib.algorithms.randomizedEigensolver.doublePassG(A, B, Binv, Omega, k, s=1, check=False)

The double pass algorithm for the GHEP as presented in [2].

Inputs:

• A: the operator for which we need to estimate the dominant generalized eigenpairs.
• B: the right-hand side operator.
• Binv: the inverse of the right-hand side operator.
• Omega: a random gassian matrix with \(m \geq k\) columns.
• k: the number of eigenpairs to extract.
• s: the number of power iterations for selecting the subspace.

Outputs:

• d: the estimate of the \(k\) dominant eigenvalues of \(A\).
• U: the estimate of the \(k\) dominant eigenvectors of \(A,\, U^T B U = I_k\).

hippylib.algorithms.randomizedEigensolver.singlePass(A, Omega, k, s=1, check=False)

The single pass algorithm for the Hermitian Eigenvalues Problems (HEP) as presented in [1].

Inputs:

• A: the operator for which we need to estimate the dominant eigenpairs.
• Omega: a random gassian matrix with \(m \geq k\) columns.
• k: the number of eigenpairs to extract.

Outputs:

• d: the estimate of the \(k\) dominant eigenvalues of \(A\).
• U: the estimate of the \(k\) dominant eigenvectors of \(A,\, U^T U = I_k\).

hippylib.algorithms.randomizedEigensolver.singlePassG(A, B, Binv, Omega, k, s=1, check=False)

The single pass algorithm for the Generalized Hermitian Eigenvalues Problems (GHEP) as presented in [2].

Inputs:

• A: the operator for which we need to estimate the dominant generalized eigenpairs.
• B: the right-hand side operator.
• Binv: the inverse of the right-hand side operator.
• Omega: a random gassian matrix with \(m \geq k\) columns.
• k: the number of eigenpairs to extract.
• s: the number of power iterations for selecting the subspace.

Outputs:

• d: the estimate of the \(k\) dominant eigenvalues of \(A\).
• U: the estimate of the \(k\) dominant eigenvectors of \(A,\, U^T B U = I_k\).

hippylib.algorithms.steepestDescent module

class hippylib.algorithms.steepestDescent.SteepestDescent(model, parameters=<hippylib.utils.parameterList.ParameterList object>)

Prior-preconditioned Steepest Descent to solve constrained optimization problems in the reduced parameter space. Globalization is performed using the Armijo sufficient reduction condition (backtracking). The stopping criterion is based on a control on the norm of the gradient.

The user must provide a model that describes the forward problem, cost functionals, and all the derivatives for the gradient.

More specifically the model object should implement following methods:

• generate_vector() -> generate the object containing state, parameter, adjoint.
• cost(x) -> evaluate the cost functional, report regularization part and misfit separately.
• solveFwd(out, x,tol) -> solve the possibly non linear forward problem up to tolerance tol.
• solveAdj(out, x,tol) -> solve the linear adjoint problem.
• evalGradientParameter(x, out) -> evaluate the gradient of the parameter and compute its norm.
• Rsolver() –> A solver for the regularization term.

Type help(Model) for additional information.

Initialize the Steepest Descent solver. Type SteepestDescent_ParameterList().showMe() for list of default parameters and their descriptions.

solve(x)

Solve the constrained optimization problem with initial guess x = [u,a,p]. Return the solution [u,a,p].

Note

x will be overwritten.

termination_reasons = ['Maximum number of Iteration reached', 'Norm of the gradient less than tolerance', 'Maximum number of backtracking reached']

hippylib.algorithms.steepestDescent.SteepestDescent_ParameterList()

hippylib.algorithms.traceEstimator module

class hippylib.algorithms.traceEstimator.TraceEstimator(A, solve_mode=False, accurancy=0.1, init_vector=None, random_engine=<function rademacher_engine>, mpi_comm=<sphinx.ext.autodoc.importer._MockObject object>)

An unbiased stochastic estimator for the trace of \(A,\, d = \sum_{j=1}^k (v_j, A v_j)\), where

• \(v_j\) are i.i.d. Rademacher or Gaussian random vectors.
• \((\cdot,\cdot)\) represents the inner product.

The number of samples \(k\) is estimated at run time based on the variance of the estimator.

Reference: Haim Avron and Sivan Toledo, Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix, Journal of the ACM (JACM), 58 (2011), p. 17.

Constructor:

• A: an operator
• solve_mode: if True we estimate the trace of A\(^{-1}\), otherwise of A.

• code:accurancy: we stop when the standard deviation of the estimator is less then

  accurancy`*tr(:code:`A).

• init_vector: use a custom function to initialize a vector compatible with the

  range/domain of A.

• random_engine: which type of i.i.d. random variables to use (Rademacher or Gaussian).

hippylib.algorithms.traceEstimator.gaussian_engine(v)

Generate a vector of \(n\) i.i.d. standard normal variables.

hippylib.algorithms.traceEstimator.rademacher_engine(v)

Generate a vector of \(n\) i.i.d. Rademacher variables.

Module contents

hippylib.mcmc package

Submodules

hippylib.mcmc.chain module

class hippylib.mcmc.chain.MCMC(kernel)

Bases: object

consume_random()

run(m0, qoi=None, tracer=None)

class hippylib.mcmc.chain.NullQoi

Bases: object

eval(x)

class hippylib.mcmc.chain.SampleStruct(kernel)

assign(other)

hippylib.mcmc.diagnostics module

hippylib.mcmc.diagnostics._acorr(mean_free_samples, lag, norm=1)

hippylib.mcmc.diagnostics._acorr_vs_lag(samples, max_lag)

hippylib.mcmc.diagnostics.integratedAutocorrelationTime(samples, max_lag=None)

hippylib.mcmc.kernels module

class hippylib.mcmc.kernels.ISKernel(model, nu)

consume_random()

delta(sample)

derivativeInfo()

init_sample(current)

name()

proposal(current)

sample(current, proposed)

class hippylib.mcmc.kernels.MALAKernel(model)

acceptance_ratio(origin, destination)

consume_random()

derivativeInfo()

init_sample(s)

name()

proposal(current)

sample(current, proposed)

class hippylib.mcmc.kernels.gpCNKernel(model, nu)

Reference:

F. J. PINSKI, G. SIMPOSN, A. STUART, H. WEBER, Algorithms for Kullback-Leibler Approximation of Probability Measures in Infinite Dimensions, http://arxiv.org/pdf/1408.1920v1.pdf, Alg. 5.2

consume_random()

delta(sample)

derivativeInfo()

init_sample(current)

name()

proposal(current)

sample(current, proposed)

class hippylib.mcmc.kernels.pCNKernel(model)

consume_random()

derivativeInfo()

init_sample(current)

name()

proposal(current)

sample(current, proposed)

hippylib.mcmc.tracers module

class hippylib.mcmc.tracers.FullTracer(n, Vh, par_fid=None, state_fid=None)

Bases: object

append(current, q)

class hippylib.mcmc.tracers.NullTracer

Bases: object

append(current, q)

class hippylib.mcmc.tracers.QoiTracer(n)

Bases: object

append(current, q)

Module contents

hippylib.utils package

Submodules

hippylib.utils.checkDolfinVersion module

hippylib.utils.checkDolfinVersion.checkdlversion()

Check if FEniCS version is supported. Currently hIPPYlib requires FEniCS version 1.6.0 and newer.

hippylib.utils.checkDolfinVersion.dlversion()

hippylib.utils.parameterList module

class hippylib.utils.parameterList.ParameterList(data)

Bases: object

A small abstract class for storing parameters and their description. This class will raise an exception if the key one tries to access is not present.

data is a dictionary where each value is the pair (value, description)

showMe(indent='')

hippylib.utils.vector2function module

hippylib.utils.vector2function.vector2Function(x, Vh, **kwargs)

Wrap a finite element vector x into a finite element function in the space Vh. kwargs is optional keywords arguments to be passed to the construction of a dolfin Function.

hippylib.utils.random module

class hippylib.utils.random.Random(myid=0, nproc=1, blocksize=1000000, seed=1)

Bases: sphinx.ext.autodoc.importer._MockObject

This class handles parallel generation of random numbers in hippylib.

Create a parallel random number number generator.

INPUTS:

• myid: id of the calling process.
• nproc: number of processor in the communicator.
• blocksize: number of consecutive random number to be generated before jumping headed in the stream.
• seed: random seed to initialize the random engine.

normal(sigma, out=None)

Sample from normal distribution with given variance.

normal_perturb(sigma, out)

Add a normal perturbation to a Vector/MultiVector.

rademacher(out=None)

Sample from Rademacher distribution.

uniform(a, b, out=None)

Sample from uniform distribution.

hippylib.utils.random.parRandom

This class handles parallel generation of random numbers in hippylib.

hippylib.utils.nb module

hippylib.utils.nb._mesh2triang(mesh)

hippylib.utils.nb._mplot_cellfunction(cellfn)

hippylib.utils.nb._mplot_function(f, vmin, vmax, logscale)

hippylib.utils.nb.animate(Vh, state, same_colorbar=True, colorbar=True, subplot_loc=None, mytitle=None, show_axis='off', logscale=False)

Show animation for a :code:TimeDependentVector

hippylib.utils.nb.coarsen_v(fun, nx=16, ny=16)

hippylib.utils.nb.multi1_plot(objs, titles, same_colorbar=True, show_axis='off', logscale=False)

Plot a list of generic dolfin object in a single row

hippylib.utils.nb.plot(obj, colorbar=True, subplot_loc=None, mytitle=None, show_axis='off', vmin=None, vmax=None, logscale=False)

Plot a generic dolfin object (if supported)

hippylib.utils.nb.plot_eigenvalues(d, mytitle=None, subplot_loc=None)

Plot eigenvalues

hippylib.utils.nb.plot_eigenvectors(Vh, U, mytitle, which=[0, 1, 2, 5, 10, 15])

Plot specified vectors in a :code:MultiVector

hippylib.utils.nb.plot_pts(points, values, colorbar=True, subplot_loc=None, mytitle=None, show_axis='on', vmin=None, vmax=None, xlim=(0, 1), ylim=(0, 1))

Plot a cloud of points

hippylib.utils.nb.show_solution(Vh, ic, state, same_colorbar=True, colorbar=True, mytitle=None, show_axis='off', logscale=False, times=[0, 0.4, 1.0, 2.0, 3.0, 4.0])

Plot a :code:TimeDependentVector at specified time steps

Module contents

Module contents

hIPPYlib implements state-of-the-art scalable algorithms for PDE-based deterministic and Bayesian inverse problems. It builds on FEniCS (http://fenicsproject.org/) (a parallel finite element element library) for the discretization of the PDE and on PETSc (http://www.mcs.anl.gov/petsc/) for scalable and efficient linear algebra operations and solvers.

For building instructions, see the file INSTALL. Copyright information and licensing restrictions can be found in the file COPYRIGHT.

The best starting point for new users interested in hIPPYlib’s features are the interactive tutorials in the notebooks folder.

Conceptually, hIPPYlib can be viewed as a toolbox that provides the building blocks for experimenting new ideas and developing scalable algorithms for PDE-based deterministic and Bayesian inverse problems.