Modeling layer

The forward problem, priors, misfits, observation operators and the model.

hippymfem.modeling.variables

Indices of the three variables of a PDE-constrained inverse problem.

hippymfem.modeling.variables.STATE = 0

the state (forward solution)

hippymfem.modeling.variables.PARAMETER = 1

the parameter (inversion variable)

hippymfem.modeling.variables.ADJOINT = 2

the adjoint variable

hippymfem.modeling.variables.NVAR = 3

number of variables

hippymfem.modeling.variables.name(i)

Human-readable name of a variable index.

class hippymfem.modeling.variables.Variables(state=None, parameter=None, adjoint=None, *extra)

Bases: list

The [state, parameter, adjoint] triple of a PDE-constrained problem.

A list, so everything written for hIPPYlib’s convention keeps working: x[STATE], x[PARAMETER] = m, len(x), extra trailing entries. The three slots are also attributes, x.state, x.parameter, x.adjoint, readable and assignable, so application code need not import the indices. generate_vector() on the problems, the model and the QoI map returns one.

property state

the state, x[STATE]

property parameter

the parameter, x[PARAMETER]

property adjoint

the adjoint, x[ADJOINT]

hippymfem.modeling.PDEProblem

Abstract interface for the forward problem.

The signatures match hIPPYlib’s PDEProblem so that the algorithms, the model, and user driver scripts port across with no changes.

class hippymfem.modeling.PDEProblem.PDEProblem

Bases: object

Forward PDE, its adjoint, and the derivative blocks the Hessian needs.

Besides the interface, this carries what every concrete problem needs the same way: the four linear solvers as assignable attributes with a common default, vector generation from the three shapes, and the release of a solver’s operator before it is rebuilt.

solver = None

linear solvers; assignable, as in hIPPYlib. solver is the forward solve, solver_adj the adjoint’s own (it shares the forward preconditioner when the Jacobian is symmetric), solver_fwd_inc/solver_adj_inc the incremental systems’. None means “build the default on first use”.

solver_adj = None

solver_fwd_inc = None

solver_adj_inc = None

generate_state()

A vector in the shape of the state.

generate_parameter()

A vector in the shape of the parameter.

generate_adjoint()

A vector in the shape of the adjoint.

generate_vector(component='ALL')

A vector, or the Variables triple, of the problem.

spd_jacobian = False

Declare dR/du symmetric positive definite. The default solvers are then CG with BoomerAMG instead of GMRES: the same iterations at a lower cost per solve, markedly so with hypre on a device. Off by default because the symmetry probe cannot tell definite from indefinite, and CG on an indefinite Jacobian fails rather than converging slowly.

SOLVER_ATTRIBUTES = ('solver', 'solver_adj', 'solver_fwd_inc', 'solver_adj_inc')

forward, adjoint, and the two incremental ones

Type:

the solver attributes a problem holds

set_solvers(solver, *args, attributes=None, **kwargs)

Install one kind of linear solver for every solve this problem makes.

solver is either a factory, such as auto_solver() or KrylovSolver, called once per attribute with *args and the keywords the factory takes (each solve then gets a solver and hierarchy of its own), or a solver instance, which serves the forward solve and is cloned for the others. Keywords that name solver parameters (rel_tolerance, max_iter, amg_relax_type, …) are set on every solver that has them instead, so rel_tolerance=1e-13 on a direct solver is ignored. attributes restricts the install to some of SOLVER_ATTRIBUTES (the default is all four). For example:

pde.set_solvers(hm.auto_solver, Vu, comm, max_direct=0,
                rel_tolerance=1e-13, max_iter=3000)

Returns self.

init_parameter(m)

Resize m to the shape of the parameter.

solveFwd(state, x)

Solve the forward problem: given \(m\), find \(u\) with \(\delta_p F(u,m,p;\hat p) = 0\) for all \(\hat p\).

solveAdj(adj, x, adj_rhs)

Solve the adjoint problem: given \(u, m\), find \(p\) with \(\delta_u F(u,m,p;\hat u) = \mathrm{adj\_rhs}\).

evalGradientParameter(x, out)

out = \(\delta_m F(u,m,p;\hat m)\).

setLinearizationPoint(x, gauss_newton_approx)

Fix the point at which the incremental solves and second derivatives are evaluated.

solveIncremental(out, rhs, is_adj)

Solve the incremental forward (is_adj=False) or adjoint system.

apply_ij(i, j, dir, out)

out = \(\delta_{ij} F(u,m,p;\hat i,\tilde j)\) with \(\tilde j\) = dir.

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

out = the third derivative block contracted with jdir and kdir.

apply_third_dir(i, x, dirs, weights, out)

out = \(\sum_m w_m \sum_{j,k} \delta_{ijk} F[t_{m,j}, t_{m,k}]\), the weighted second directional derivatives of the slot-i gradient along directions t_m = (u, m, p) (None for an absent component).

This default sums apply_ijk() over the ordered pairs of nonzero components; a problem that can differentiate a whole direction at once overrides it.

hippymfem.modeling.PDEVariationalProblem

Stationary forward problem from a pointwise residual density.

The user writes the weak residual of the forward PDE as a density at one quadrature point; every operator the inverse problem needs follows by automatic differentiation (see hippymfem.fem.kernel).

The essential-dof treatment of each block (A with a unit diagonal, W_uu with a zero one, the rectangular blocks on the test side, the transposes carrying it across) is the table in hippymfem.fem.bcs; with it baked into the matrices, apply_ij is a plain matrix product with no masking.

hippymfem.modeling.PDEVariationalProblem.DIRECT_TRANSPOSES = 'auto'

Take the linearization point’s blocks from one differentiation pass instead of one per block. The blocks are bit-identical either way (for an element batch that is not split) and the shared pass is faster; HIPPYMFEM_SHARE_HESSIAN=0 turns it off. Whether the two blocks that every CG iteration applies transposed, C^T and W_um^T, are assembled as matrices of their own from the transposed element matrices of the same differentiation pass, rather than formed by MFEM at the first transposed product (HypreParMatrix::MultTranspose). With hypre on a device MFEM forms each such transpose by sorting, which asked the driver for 4.9 GiB a block at 120^3 quadratic hexahedra a GCD (2.4 GiB of them kept), and cost 1 to 2.4 % of two Newton-CG steps; on a host MFEM multiplies with the transpose without forming it, and a second matrix would only take memory. auto (default): with hypre on a device, 1 always, 0 never (HIPPYMFEM_DIRECT_TRANSPOSES). A problem with interior-facet terms keeps MFEM’s.

hippymfem.modeling.PDEVariationalProblem.set_share_hessian_pass(flag=True)

Turn the shared Hessian pass on or off; returns the old value.

class hippymfem.modeling.PDEVariationalProblem.PDEVariationalProblem(Vh, varf_handler, bc=None, bc0=None, is_fwd_linear=False, quadrature_degree=None, aux_spaces=(), bdr_varf=None, bdr_attributes='all', bdr_quadrature_degree=None, symmetric_jacobian='auto', release_linearization_on_move=False, transpose_free_adjoint=False, spd_jacobian=False, facet_varf=None, facet_quadrature_degree=None)

Bases: PDEProblem, KeepAlive

Forward problem defined by a pointwise residual density.

Parameters:

• Vh (sequence of 3 spaces) – [STATE, PARAMETER, ADJOINT]; each a FunctionSpace or a raw ParFiniteElementSpace.

• varf_handler (callable) – varf_handler(u, m, p, x) -> scalar, JAX-traceable, with Field arguments. This is the weak residual density: the forward problem is dR/dp = 0.

• bc (DirichletBC, list, BCSet, or None) – Essential conditions with data, for the forward solve.

• bc0 (same) – Homogeneous counterpart, for the adjoint and incremental solves.

• is_fwd_linear (bool) – Declares the residual linear in u. One Newton step is then exact; the residual is still checked, so a mislabeled nonlinear problem is reported rather than silently mis-solved.

• quadrature_degree (int, optional) – Integration rule degree. Defaults to 2*max(order) + 2. Every derivative block uses this one rule, so they stay mutually consistent.

• aux_spaces (sequence of spaces, optional) – Extra fields appended to varf_handler’s signature; set their values through aux_values.

• bdr_varf (callable, optional) – bdr_varf(u, m, p, x, n, *aux) -> scalar, a boundary residual density, with n the outward unit normal. This is the ds of the UFL form: Neumann and Robin conditions, boundary sources, penalty terms. It is differentiated exactly like the domain density, so every block picks up its boundary part; see hippymfem.fem.boundary. The fields carry their full gradient there, so dot(u.grad, n) works.

• bdr_attributes (sequence of int or "all", optional) – Which boundary attributes bdr_varf applies to (MFEM’s 1-based numbering).

• bdr_quadrature_degree (int, optional) – Rule degree on the boundary; defaults to quadrature_degree.

• facet_varf (callable, optional) – facet_varf(u, m, p, x, n, h, *aux) -> scalar, an interior facet residual density, the dS of a discontinuous Galerkin form. Each field arrives as a FacetField carrying its trace from both sides, so the form is written in jump and avg; n is the normal out of the first element and h the pair of face measures. It is differentiated like the domain density, so the residual and every block pick up their facet part; see hippymfem.fem.facets.

• facet_quadrature_degree (int, optional) – Rule degree on interior faces; defaults to quadrature_degree.

• symmetric_jacobian (bool or "auto", optional) – Whether dR/du is symmetric; "auto" (the default) probes the assembled matrix. See symmetric_jacobian.

• spd_jacobian (bool, optional) – Declare dR/du symmetric positive definite, so the default solvers use CG with BoomerAMG instead of GMRES (see spd_jacobian). Implies symmetric_jacobian=True.

• release_linearization_on_move (bool) – Drop the linearization point at a forward solve for a different parameter; see release_linearization_on_move.

• transpose_free_adjoint (bool) – Solve the adjoint through A’s MultTranspose with A’s own AMG hierarchy, without forming A^T; see transpose_free_adjoint.

symmetric_jacobian

Whether dR/du is symmetric. If it is, the adjoint and adjoint-incremental solves reuse the Jacobian and its preconditioner, saving a transposed copy, a second AMG hierarchy and one BoomerAMG setup per Newton step.

"auto", the default, decides from the assembled matrix rather than from a promise about the PDE: every newly assembled Jacobian is probed with _probe_symmetry() (four matrix-vector products, far cheaper than the AMG setup they can save), so a residual that is self-adjoint at one parameter and not at another is handled correctly. True and False declare it outright and skip the probe. With single-precision element matrices (HIPPYMFEM_PRECISION=mixed) the verdict on the first Jacobian is kept, since a symmetrized matrix passes any probe: such a residual has to say False there (a refined solve that stalls is then done again in double precision, with a warning). With single-precision solves alone the kept verdict is checked on every new Jacobian.

release_linearization_on_move

At a forward solve for a parameter other than the linearization point’s, release that point (C, the W blocks, A/At and the incremental solvers’ hierarchies) before the new Jacobian is assembled, so a line search never holds the old point’s operators beside the trial point’s. Opt-in, because it is a claim about the caller: a line-search Newton-CG has its direction before it moves and rebuilds the point after, but a trust-region step, a finite-difference check or anything that applies the Hessian after a trial solve needs the old point to survive it. See the GPU guide.

transpose_free_adjoint

Solve with A^T through A’s MultTranspose and A’s own AMG hierarchy as the preconditioner, instead of forming the transpose and a second hierarchy. For a residual that is not linear in the state a Newton step then holds A and one hierarchy instead of A, A^T and two, which can decide whether a large problem fits in device memory. GMRES accepts A’s V-cycle as a preconditioner for A^T; it is as good as the asymmetry is small (for example the k'(u) term of a temperature-dependent conductivity). Krylov adjoint solvers only; a direct solver keeps the explicit transpose. A Jacobian that is not symmetric then stays in double precision also where a single-precision hypre is loaded: there is no transpose wrapper of a matrix of that library, and it would need the transpose and a second hierarchy that this setting is asked to avoid.

aux_values

true-dof vectors for the auxiliary fields, in declaration order

Ct

C^T and W_um^T as matrices of their own (DIRECT_TRANSPOSES), or None

newton_parameters

Newton controls for a nonlinear forward solve. It ends when the residual is below max(rel_tolerance * ||r_0||, abs_tolerance), or when the Newton correction is below step_tolerance of the state (solveFwd(); with 0 the residual alone decides).

generate_state()

A vector in the shape of the state.

generate_parameter()

A vector in the shape of the parameter.

generate_adjoint()

A vector in the shape of the adjoint.

init_parameter(m)

Resize m to the shape of the parameter.

single_solves = True

Whether this problem’s solves may run in a single-precision hypre when one is named (HIPPYMFEM_HYPRE_SINGLE, hippymfem.algorithms.singlesolve). Change it with set_single_solves(), which drops what was built for the other precision.

set_single_solves(on=True)

Run the solves of this problem in the single-precision hypre (True, when one is named) or in MFEM’s double-precision one (False) from now on; returns the old choice.

The Jacobian, the solvers built on it and the linearization point are dropped when the choice changes: the next forward solve assembles the Jacobian into the library chosen, and setPointForHessianEvaluations must be called again before the next Hessian action. The stages of a Laplace approximation keep the single-precision solves well: at 64³ (doublePassG, k = 50) the eigenvalues came out within 1e-5 of double precision, the pointwise posterior variance within 6e-6, in half the time of the eigensolver; switch them off where eigenvalues are wanted to more digits than that.

SYMMETRY_PROBE_TOL = 1e-13

relative tolerance of the symmetry probe, and how many vector pairs it uses

SYMMETRY_PROBE_VECTORS = 2

SYMMETRY_PROBE_TOL_SINGLE = 1e-08

the matrix of a symmetric form then has a skew part of about 1e-7 of its entries, which the probe sees at 1e-7 / sqrt(n); a form that is not symmetric shows up at about 1 / sqrt(n), far above this for any size that fits a GPU

Type:

the same when the element matrices are computed in single precision

SYMMETRY_PROBE_TOL_SINGLE_PRODUCTS = 1e-05

the same, times sqrt(n), for a matrix of the single-precision hypre, whose products are rounded to single precision: the probe sees a symmetric form there at about 2e-7 / sqrt(n), and a skew part of a fraction s of the matrix at about s / sqrt(n) (an advection term next to diffusion: 2e-3 on 21^3 dofs, falling with the mesh size), so a skew part above about 1e-5 of the matrix is caught

REFINE_MAX_STEPS = 5

Iterative refinement of the forward and the adjoint solve when the Jacobian is assembled in single precision and the residuals in double ("mixed"): the most corrections after the first solve, and the relative accuracy that one solve with the single-precision matrix is asked for. Each solve reduces the double-precision residual by about that factor, so two solves at 1e-6 reach what one solve at 1e-12 reaches with a double-precision matrix, in about as many Krylov iterations.

REFINE_INNER_TOL = 1e-06

REFINE_PREDICT = True

The last pass of a refinement is not followed by an evaluation of the residual when the passes before it predict, with this margin, that it reaches the goal: a pass leaves of its right-hand side the tolerance it is asked for times a ratio that the measured passes give. REFINE_PREDICT = False measures every pass.

REFINE_MARGIN = 3.0

SINGLE_REFINE_GOAL = 0.0

With the Jacobian in a single-precision hypre, a relative residual at which the refinement of the forward and the adjoint solve may stop although the solver’s own tolerance asks for less. 0, the default, honours the solver’s tolerance (three passes for 1e-12). 1e-9 stops after two passes and one evaluation of the residual: the gradient is then exact to about nine digits, which Newton-CG with a tolerance of 1e-6 does not notice (the same steps, the same cost functional to nine digits), while an optimizer run to a tighter tolerance, such as BFGS to 1e-8, ends in a line search that finds no decrease.

REFINE_STALL = (1000.0, 1e-08)

A refinement that ends above this many times its goal and above this relative residual (or a hundred times the floor of single-precision element vectors) has not converged: single precision cannot carry the solve at this point (a parameter that makes the Jacobian too ill-conditioned for it, as a line search may propose), or the matrix is not the operator of its residual (a Jacobian taken for symmetric where it is not, symmetric_jacobian). The system is then solved once more with the Jacobian assembled in double precision (_solve_in_double()), with a RuntimeWarning once per process; so is one whose single-precision solver gives up in a pass. A refinement that converges ends at its goal (about 1e-11 for a tolerance of 1e-12) in two or three passes.

SINGLE_GIVE_UP = 'auto'

What a forward solve does where single precision cannot carry it (a refinement that stalls, or whose solver gives up): "double" solves it again with the Jacobian in double precision, as above; "raise" raises the RuntimeError, which a line search takes for a failed trial and answers by backtracking. The second is for a run whose card has no room for a double-precision Jacobian and its hierarchy beside the single-precision ones, which the first holds at once (with 1.7 million quadratic hexahedra a rank on a 64 GiB card, 36 GiB beside 12: such a trial ended the run, out of device memory); a forward solve outside a line search then fails where it would have been solved. An adjoint solve is solved in double precision either way: it is at a point the line search accepted. "auto" (the default since November 2026) does the first where a double-precision Jacobian and its hierarchy fit beside what the device holds, on every rank (_double_fits()), and the second where they do not; with hypre in host memory, or where the runtime does not say, it does the first. (checked at import, as the other settings read from the environment: a misspelt value was found only when a forward solve gave up, as a ValueError that a line search does not take for a failed trial)

GIVE_UP_BYTES_PER_NNZ = 50.0

Bytes of the double-precision Jacobian and its hierarchy that "auto" counts for each nonzero of the single-precision Jacobian: about 43 with quadratic hexahedra on Frontier (36 GiB for some 0.9 billion nonzeros, 120^3 elements a GCD), with a margin.

invalidate_jacobian()

Drop the cached Jacobian; call this if a coefficient changed in place.

release_linearization_point()

Drop every operator of the linearization point; the next use must set one.

One gc.collect at the end frees the dropped hypre matrices before the next point’s blocks exist, which lowers the memory peak.

solveFwd(state, x)

Newton solve of dR/dp = 0 for the state, given x[PARAMETER].

A linear residual converges in one step; the declaration is_fwd_linear only stops the loop early, it does not change the step.

A nonlinear solve ends when its residual is below the tolerance, and also when the Newton correction is below newton_parameters["step_tolerance"] of the state. The residual has a floor of round-off that grows with the mesh, and a tolerance may lie under it: 1e-14 against 1.2e-14 at 135 million unknowns, where a solve that a line search had warm-started ran to max_iter and raised with a state that had all its digits. The correction measures the error of the state it corrects, whatever the size of the residual. A solve whose own line search finds no step that lowers the residual, with a correction that is not negligible, raises at once and does not repeat that search max_iter times.

solveAdj(adj, x, adj_rhs)

Solve A^T p = adj_rhs with homogeneous essential conditions.

evalGradientParameter(x, out)

out = dR/dm at x; equals C^T p but needs no matrix.

setLinearizationPoint(x, gauss_newton_approx=False, matrix_free=False)

Assemble the blocks used by the incremental solves and the Hessian.

With matrix_free=True only the Jacobian is assembled (the incremental solves need it), and apply_ij() computes the products with the other blocks from the element kernels at x (apply_ij_at()); a collection of variables, matrix_free=(PARAMETER,) say, does that for the blocks that involve one of them and assembles the others. No block matrices and no scatter patterns for them, for one kernel pass per product: the trade for a problem that is short of memory, or for a point whose blocks are applied once or twice (the CG iterations of a Newton step apply them many times, and on a problem that fits they are cheaper assembled).

solveIncremental(out, rhs, is_adj, rel_tolerance=None)

Solve the incremental forward or adjoint system. rel_tolerance replaces the solver’s own relative tolerance for this solve where it is the looser.

block_is_zero(i, j)

Whether apply_ij() for (i, j) would return zero whatever the direction, at the current linearization point: a block the point does not hold (Gauss-Newton, a residual linear in the state) or, without assembled blocks, one that is known to vanish.

apply_ij(i, j, dir, out)

Apply the (i, j) second-derivative block to dir.

apply_ij_at(i, j, x, dir, out)

The (i, j) second-derivative block at x applied to dir, from the element kernels and without assembling the block: what apply_ij() gives after setLinearizationPoint() at x.

The essential dofs are treated as the assembled blocks treat them: a direction in the state or the adjoint loses its essential entries, the product in the state or the adjoint has them zeroed, and the Jacobian’s pair puts the identity back on its essential rows. One forward tangent through the element gradient per product, where assembling takes a tangent per element dof and a scatter: the cheaper route for a block applied once or twice at a point.

Like the residual and the gradients, a product is assembled without communication, each rank taking a face it shares with another from its own side; so an interior-facet density must be unchanged when the two sides trade places (as a DG form is, its jumps and normal changing sign together). The assembled blocks do not need that for their transposes; the products do.

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

Third-derivative block contracted with jdir and kdir.

apply_third_dir(i, x, dirs, weights, out)

Weighted second directional derivatives of the slot-i gradient, \(\sum_m w_m\, D^2(\partial_i R)[t_m, t_m]\).

dirs holds one direction t_m per weight w_m, each a sequence (u, m, p) of vectors (or longer, when Vh lists further variables), None where the direction has no component. It equals the sum of apply_ijk() over every ordered pair of nonzero components, \(\sum_{j,k} R_{ijk}[t_j, t_k]\), but the quadrature kernel differentiates the whole direction at once: one pass instead of up to nine, which is what a second-order adjoint over many modes spends its time in.

functional(x)

The integral of the residual density itself, R(u, m, p).

hippymfem.modeling.PDEVariationalProblem.require_finite_arrays(arrays, comm, what)

Raise if any element array has a non-finite entry, on every rank.

Checking the inputs is not enough: a parameter of 800 is finite while exp(m) is not, so the overflow happens inside the density and only the element arrays show it. Unchecked, the matrix reaches hypre, whose AMG setup either hangs or aborts the whole job (Error during setup! Error code: 12), neither of which a line search can recover from. The cost is one allreduce per assembly.

hippymfem.modeling.PDEVariationalProblem.require_finite(v, what)

Raise if v has any non-finite entry, on every rank.

A matrix assembled from non-finite values has inf or nan entries, on which hypre’s AMG setup hangs rather than fails. Checking first turns that into a RuntimeError, which the line searches handle by backtracking; the check costs one allreduce.

class hippymfem.modeling.PDEVariationalProblem.TransposeOf(*args, **kwargs)

Bases: TransposeOperator

A^T as an MFEM operator over A.MultTranspose, no transpose formed.

transposed_of is A itself: MFEM holds a raw pointer to it, so the wrapper keeps it alive, and _set_operator_once() reads it to hand over A’s preconditioner instead of building one on the wrapper.

hippymfem.modeling.TimeDependentPDEVariationalProblem

Time-dependent forward problems from a one-step residual density.

The user writes the residual of a single time step,

\[r_n(u_n, u_{n-1}, m, p_n; x, t_n, \Delta t) = 0,\]

as a pointwise density. Implicit Euler for a parabolic problem is

\[r_n = \frac{u_n - u_{n-1}}{\Delta t}\,p_n + a(u_n, m; p_n) - f(t_n)\,p_n .\]

Because the previous level \(u_{n-1}\) is a differentiation slot of the kernel and not a frozen coefficient, every operator the space-time inverse problem needs is an exact derivative of that one density:

block	what it is
\(A_n = \partial_{u_n}\partial_{p}r_n\)	the step Jacobian
\(B_n = \partial_{u_{n-1}}\partial_{p}r_n\)	coupling to the level before
\(C_n = \partial_{m}\partial_{p}r_n\)	the control block
\(W^{(n)}_{ab}\)	second derivatives in \(a,b \in \{u_n, u_{n-1}, m\}\)

The forward problem marches \(A_n u_n = \ldots\) forward; the adjoint marches \(A_n^{\!\top} p_n = \text{rhs}_n - B_{n+1}^{\!\top} p_{n+1}\) backward, the \(B^{\!\top}\) term being exactly where reverse-in-time coupling comes from. The second-order blocks are accumulated step by step, each step writing into the two time levels it touches.

Operators are built once per trajectory, and steps that agree share them. A step Jacobian that is linear in \(u_n\) with coefficients independent of \(t\) (implicit Euler for a linear PDE, the common case) is the same matrix at every step. Nothing here assumes that: every step is assembled and compared with the step before it through a random matrix-vector product, and a run of equal steps keeps one \(A\), one \(A^{\!\top}\), one \(B\) and one solver with its multigrid hierarchy. When all steps agree, one of each serves the whole march. Agreement at the first steps proves nothing about the later ones (a state at rest until a source switches on, a coefficient that changes in time), so every step is compared. Either way no operator is assembled or a hierarchy set up inside a solve loop, where a step-wise implementation would set BoomerAMG up once per step inside every Newton iteration of the forward solve and every CG iteration of the incremental solves. The blocks of a linearization point come from one differentiation pass per step rather than one per block.

The trajectory lives in a TimeDependentVector; the time and step size reach the kernel as traced scalars, so the compiled kernels are reused across all steps.

hippymfem.modeling.TimeDependentPDEVariationalProblem.NEW = 0

u_n

hippymfem.modeling.TimeDependentPDEVariationalProblem.OLD = 1

u_{n-1}

hippymfem.modeling.TimeDependentPDEVariationalProblem.PARAM = 2

m

hippymfem.modeling.TimeDependentPDEVariationalProblem.ADJ = 3

p_n

class hippymfem.modeling.TimeDependentPDEVariationalProblem.TimeDependentPDEVariationalProblem(Vh, varf_handler, bc, bc0, u0, t_init, t_final, dt, is_fwd_linear=False, quadrature_degree=None, spd_jacobian=False, bdr_varf=None, bdr_attributes='all', bdr_quadrature_degree=None, init_varf=None, init_linear=True)

Bases: PDEProblem, KeepAlive

One-step time-dependent forward problem.

Parameters:

• Vh (sequence of 3 spaces) – [STATE, PARAMETER, ADJOINT], the spatial spaces for one time level.

• varf_handler (callable) – varf(u, u_old, m, p, x, t, dt) -> scalar, JAX-traceable, with Field arguments for the four fields.

• bdr_varf (callable, optional) – bdr_varf(u, u_old, m, p, x, n, t, dt) -> scalar, a boundary density of the same step residual with the outward unit normal n, integrated over the boundary elements of bdr_attributes and added to the domain term: a Robin condition, a prescribed flux, a boundary source. Its derivative blocks are added to the step blocks before the essential rows are eliminated, as in PDEVariationalProblem.

• bdr_attributes ("all" or sequence of int) – Which boundary attributes bdr_varf applies to (MFEM’s 1-based attributes); all of them by default.

• bdr_quadrature_degree (int, optional) – Quadrature degree on the boundary; quadrature_degree by default.

• bc (boundary conditions on the state, with data and homogeneous.)

• bc0 (boundary conditions on the state, with data and homogeneous.)

• u0 (ParVector) – Initial condition at t_init, fixed; None with init_varf.

• init_varf (callable, optional) – init_varf(u0, m, p, x) -> scalar: the initial density, the integrand of a weak equation for the initial state, r_0(u_0, m; p) = 0, so that the initial state may depend on the parameter: (u0.val - g(m.val, x)) * p.val is the projection of a pointwise map g (the parameter as the initial state where it is m.val). Every operator of that equation (its Jacobian, its coupling to the parameter and its second derivatives) comes from automatic differentiation like those of the steps; the adjoint, the gradient and the Hessian then reach back to the initial level. It is a PDEVariationalProblem of its own, init_problem, with the state’s boundary conditions (its solvers can be set there).

• init_linear (bool) – The initial density is linear in u0 (default), as a projection is.

• t_init (float)

• t_final (float)

• dt (float)

• is_fwd_linear (bool) – The step residual is linear in u_n; one Newton step is then exact.

• spd_jacobian (bool) – The step Jacobian is symmetric positive definite; the default solvers are then CG with BoomerAMG (see spd_jacobian).

INVARIANCE_PROBE_TOL = 1e-13

relative tolerance of the “same operator” probe (see _same_operator)

slots

(u_n, u_{n-1}, m, p_n)

Type:

slot spaces

newton_parameters

Newton controls of a nonlinear step. A step ends when its residual is below max(rel_tolerance * ||r_0||, abs_tolerance), with r_0 the residual at the previous level, or when the Newton correction is below step_tolerance of the state (0: the residual alone decides). Near a steady state r_0 is the small change of one step, and the round-off of a residual in large units can lie above that tolerance; the correction measures the error of the state whatever the size of the residual.

init_problem

the problem of the initial level, where an initial density is given

generate_state()

A vector in the shape of the state.

generate_adjoint()

A vector in the shape of the adjoint.

generate_parameter()

A vector in the shape of the parameter.

generate_static_state()

generate_static_adjoint()

init_parameter(m)

Resize m to the shape of the parameter.

solveFwd(out, x)

March the state forward, storing every level in out.

For a residual linear in u_n the step Jacobian of every step is assembled and compared with the one the solver holds; a run of equal steps keeps one matrix and one hierarchy. A Jacobian that was the same at every step of a march is kept across forward solves at the same parameter and then serves every step without assembly.

solveAdj(out, x, adj_rhs)

March the adjoint backwards.

\(A_n^{\!\top}p_n = \text{rhs}_n - B_{n+1}^{\!\top} p_{n+1}\), with the step operators of the trajectory x built once (see _trajectory) and the adjoint solver holding one transposed Jacobian for the whole march when the steps agree.

evalGradientParameter(x, out)

out = \(\sum_n \partial_m r_n\).

setLinearizationPoint(x, gauss_newton_approx=False)

Build (or complete) the step operators at x and point the incremental solvers at them, once.

With Gauss-Newton the W blocks are not needed and not assembled; a trajectory that already has them keeps them.

solveIncremental(out, rhs, is_adj)

Solve the space-time incremental system.

Forward (is_adj=False): march \(A_n \hat u_n = \text{rhs}_n - B_n \hat u_{n-1}\) from \(\hat u_0 = 0\) for a fixed initial state, or from the solution of \(A_0 \hat u_0 = \text{rhs}_0\) with an initial density.

Adjoint (is_adj=True): march \(A_n^{\!\top}\hat p_n = \text{rhs}_n - B_{n+1}^{\!\top}\hat p_{n+1}\) backward, and with an initial density end with \(A_0^{\!\top}\hat p_0 = \text{rhs}_0 - B_1^{\!\top}\hat p_1\). The solvers were pointed at their operators by setLinearizationPoint(); nothing is set up here.

applyC(dm, out)

out_n = C_n dm at every level.

applyCt(dp, out)

out = \(\sum_n C_n^{\!\top} dp_n\).

applyWuu(du, out)

State-state block of the Lagrangian, including the time coupling.

applyWum(dm, out)

State rows of the state-parameter block.

applyWmu(du, out)

Parameter row of the parameter-state block (the transpose of Wum).

applyWmm(dm, out)

apply_ij(i, j, dir, out)

out = \(\delta_{ij} F(u,m,p;\hat i,\tilde j)\) with \(\tilde j\) = dir.

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

out = the third derivative block contracted with jdir and kdir.

release_linearization_point()

Drop the trajectory operators and the solvers’ hold on them.

exportState(u, basename, field='u')

Write a trajectory as a ParaView time series.

class hippymfem.modeling.TimeDependentPDEVariationalProblem.ImplicitEulerTimeDependentPDEVariationalProblem(Vh, varf_handler, bc, bc0, u0, t_init, t_final, dt, is_fwd_linear=False, quadrature_degree=None, spd_jacobian=False, bdr_varf=None, bdr_attributes='all', bdr_quadrature_degree=None, init_varf=None, init_linear=True)

Bases: TimeDependentPDEVariationalProblem

Alias kept for source compatibility with hIPPYlib.

The base class already covers any one-step scheme, implicit Euler included; which scheme it is lives in the residual density the user writes.

hippymfem.modeling.prior

Gaussian prior models.

Matern-class Gaussian priors whose precision is a differential operator. Writing \(\mathcal{A}\) for the square-root precision operator and \(M\) for the mass matrix, the bi-Laplacian prior has precision

\[R = \mathcal{A} M^{-1} \mathcal{A},\]

so a sample is \(s = \mathcal{A}^{-1} \sqrt{M}\,\xi\) with \(\xi \sim N(0, I)\), because then \(\mathrm{Cov}(s) = \mathcal{A}^{-1} M \mathcal{A}^{-1} = R^{-1}\).

The square root of the mass matrix, without a square root. hIPPYlib(x) gets \(\sqrt{M}\) from a quadrature element space: with a diagonal quadrature mass matrix \(M_q\) it forms \(\sqrt{M} = G M_q^{-1/2}\). The same operator has an explicit entry-wise form, which is what is used here:

\[\sqrt{M}_{i,(e,q)} = N_i(x_q^e)\,\sqrt{w_q\,|\det J^e_q|},\]

and indeed \(\sqrt{M}\sqrt{M}^{\!\top} = \sum_{e,q} w_q |\det J| N_i N_j = M\) exactly for that quadrature rule. The Laplacian prior needs the same thing for \(R = \gamma L + \delta M\) directly, with the gradient components carrying \(\sqrt{\gamma}\) and the value component \(\sqrt{\delta}\).

So no square root is ever computed, and the white-noise vector lives on quadrature points, which are element-local: drawing a prior sample needs no communication. Keying the noise on each element’s centroid, not on its partition-dependent number, makes it bit-identical under repartitioning, so results can be compared across rank counts.

hippymfem.modeling.prior.prior_h()

class hippymfem.modeling.prior.QuadratureSqrtPrecision(space, quadrature_degree, value_coeff=1.0, grad_coeff=None)

Bases: Operator

Map element-local white noise at quadrature points to the FE dual space.

Builds

\[S_{i,(e,q,c)} = c_0\,N_i(x_q)\,\sqrt{w_q|\det J|} \quad\text{and}\quad c_1\,\partial_d N_i(x_q)\,\sqrt{w_q|\det J|}\]

so that \(S S^{\!\top} = c_0^2 M + c_1^2 L\). With (value_coeff, grad_coeff) = (1, None) this is \(\sqrt{M}\); with (sqrt(delta), sqrt(gamma)) it is the square root of \(\delta M + \gamma L\).

Applied matrix-free. The noise space is purely local (no ghost dofs), and its entries are keyed on element centroids so that a sample does not depend on the partition.

noise_vector()

A zero white-noise vector (element-local layout); no communication.

sample_noise(sigma=1.0, out=None, rng=None)

Fill a noise vector with N(0, sigma^2), partition-independently.

Entries are keyed on a hash of the element centroid (see _element_keys()), so the same sample comes out on any number of ranks. Nothing here communicates.

init_vector(x, dim)

mult(noise, out)

out = S noise (quadrature space -> FE dual space).

multTranspose(v, noise_out)

noise_out = S^T v (FE space -> quadrature space).

class hippymfem.modeling.prior.SqrtPrecisionPDE_Prior(Vh, domain_integrators, bdr_integrators=(), mean=None, rel_tol=1e-12, max_iter=1000, solver_type='krylov', quadrature_degree=None, systems_dim=None, lumped_mass=False)

Bases: _Prior

Prior with precision \(R = \mathcal{A} M^{-1} \mathcal{A}\).

Parameters:

• Vh (FunctionSpace)

• domain_integrators (sequence of mfem.BilinearFormIntegrator) – The square-root precision form \(\mathcal{A}\).

• bdr_integrators (sequence, optional) – Boundary contributions to \(\mathcal{A}\) (the Robin term).

• mean (ParVector, optional)

• solver_type ({"krylov", "lu"}) – "lu" is serial only (see hippymfem.algorithms.linSolvers).

• lumped_mass (bool) – The mass matrix \(M\) of \(R = \mathcal{A} M^{-1} \mathcal{A}\) by the rule at the vertices of an element, which makes it diagonal: the sum of each row of the consistent one where the elements are affine, the lumped mass of the SPDE approach to Matern fields (Lindgren, Rue and Lindstrom, 2011). A product with \(R\), which the cost, the gradient and every Hessian action take, then solves with a diagonal matrix instead of running some seventy CG iterations on the consistent one, and the square root that samples are drawn with is by the same rule, so the samples have the covariance \(R^{-1}\) of this precision exactly. It is another discretization of the same prior, not the same matrix: off by default. For a parameter space of first order on a mesh of one kind of element.

init_vector(x, dim)

noise_vector()

A white-noise vector of the right shape for sample().

sample_noise(sigma=1.0, out=None, rng=None)

class hippymfem.modeling.prior.LaplacianPrior(Vh, gamma, delta, mean=None, rel_tol=1e-12, max_iter=1000, solver_type='krylov', quadrature_degree=None)

Bases: _Prior

Prior with precision \(R = \gamma L + \delta M\) assembled directly.

Its square root is available in closed form on quadrature points, so samples need one \(R\) solve rather than two \(\mathcal{A}\) solves.

init_vector(x, dim)

noise_vector()

A white-noise vector of the right shape for sample().

sample_noise(sigma=1.0, out=None, rng=None)

hippymfem.modeling.prior.BiLaplacianPrior(Vh, gamma, delta, Theta=None, mean=None, rel_tol=1e-12, max_iter=1000, robin_bc=False, solver_type='krylov', quadrature_degree=None, lumped_mass=False)

Matern-class prior with \(\mathcal{A} = \gamma\,\mathrm{div}(\Theta\nabla) + \delta\).

Theta may be None (isotropic), a scalar, a numpy array (a constant tensor), or an mfem.MatrixCoefficient. robin_bc=True adds \(\sqrt{\gamma\delta}/1.42\) on the boundary, which reduces the boundary-layer artefact in the pointwise variance. lumped_mass=True takes the diagonal mass matrix of the rule at the vertices for the \(M\) of \(R = \mathcal{A} M^{-1} \mathcal{A}\) (SqrtPrecisionPDE_Prior), with which a product with \(R\) needs no iteration.

hippymfem.modeling.prior.VectorBiLaplacianPrior(Vh, gamma, delta, mean=None, rel_tol=1e-12, max_iter=1000, robin_bc=False, solver_type='krylov', quadrature_degree=None)

Componentwise Matern prior on a vector space.

gamma and delta are per-component sequences whose entries must agree: MFEM’s vector integrators apply one scalar coefficient to every component, so differing values are rejected rather than silently averaged.

hippymfem.modeling.prior.MollifiedBiLaplacianPrior(Vh, gamma, delta, locations, m_true, Theta=None, pen=10.0, order=2, rel_tol=1e-12, max_iter=1000, solver_type='krylov', quadrature_degree=None)

Matern prior with a mollifier pinning the field near given locations.

Adds \(\mathrm{pen}\,\delta\, \mathrm{mollifier}\, m\, \hat m\) to the square-root precision form and sets the mean so that the prior mean matches m_true where the mollifier is active.

hippymfem.modeling.prior.BiLaplacianComputeCoefficients(sigma2, rho, ndim)

(gamma, delta) giving marginal variance sigma2 and correlation length rho for the Matern-class bi-Laplacian prior.

class hippymfem.modeling.prior.GaussianRealPrior(Vh=None, covariance=None, mean=None, comm=None)

Bases: _Prior

Gaussian prior on a finite-dimensional (Real-type) space.

The covariance is a dense n x n matrix; used for hyper-parameters and for testing, where an exact reference is wanted.

init_vector(x, dim)

noise_vector()

A white-noise vector of the right shape for sample().

sample_noise(sigma=1.0, out=None, rng=None)

trace(method='Exact', **kw)

Trace of \(R^{-1} M\), i.e. the integrated prior variance.

ensemble ("Randomized" only): a communicator whose ranks each hold a complete copy of the problem; the solves are divided among them (MatMvMult()).

pointwise_variance(method='Exact', **kw)

Diagonal of \(R^{-1}\), the pointwise prior variance.

"Exact" applies \(R^{-1}\) to every unit vector (small problems only). "Randomized" is hIPPYlib’s: the diagonal of the rank-r truncation of \(R^{-1}\) from a randomized eigensolver. It is biased low, even for large r, because a prior covariance has a slowly decaying spectrum. "MonteCarlo" averages the squares of n prior samples (n defaults to r): unbiased at one solve per sample, so it is the method to use when samples are cheap, as they are with hypre on a device.

ensemble ("MonteCarlo" and "Randomized"): a communicator whose ranks each hold a complete copy of the problem; the samples, or the solves, are divided among them, and every rank returns the field that one rank alone would compute.

hippymfem.modeling.misfit

Misfit (data likelihood) functionals.

A misfit exposes cost, grad, setLinearizationPoint and apply_ij with hIPPYlib’s signatures, so the model and the algorithms do not care which one is in use.

class hippymfem.modeling.misfit.Misfit

Bases: object

Abstract misfit term of the cost functional.

cost(x)

grad(i, x, out)

setLinearizationPoint(x, gauss_newton_approx=False)

apply_ij(i, j, dir, out)

block_is_zero(i, j)

Whether the second-derivative block (i, j) vanishes identically, so that a caller can skip applying and adding it. False unless a class knows better.

class hippymfem.modeling.misfit.DiscreteStateObservation(B, data=None, noise_variance=None)

Bases: Misfit, KeepAlive

Gaussian misfit for a linear observation operator: \(\tfrac{1}{2\sigma^2}\|Bu - d\|^2\).

cost(x)

grad(i, x, out)

setLinearizationPoint(x, gauss_newton_approx=False)

block_is_zero(i, j)

Whether the second-derivative block (i, j) vanishes identically, so that a caller can skip applying and adding it. False unless a class knows better.

apply_ij(i, j, dir, out)

hippymfem.modeling.misfit.PointwiseStateObservation(Vh, obs_points, data=None, noise_variance=None)

Misfit for pointwise observations of the state at obs_points.

class hippymfem.modeling.misfit.MultDiscreteStateObservation(B, data=None, Mpar=1.0)

Bases: Misfit, KeepAlive

Multiplicative Gamma(M, M) noise model.

The negative log-likelihood is \(M \sum_t \left(\log (Bu)_t + d_t / (Bu)_t\right)\), so the state must stay strictly positive; a nonpositive observed value is reported rather than producing a silent NaN.

cost(x)

grad(i, x, out)

setLinearizationPoint(x, gauss_newton_approx=False)

block_is_zero(i, j)

Whether the second-derivative block (i, j) vanishes identically, so that a caller can skip applying and adding it. False unless a class knows better.

apply_ij(i, j, dir, out)

hippymfem.modeling.misfit.MultPointwiseStateObservation(Vh, obs_points, Mpar, data=None)

Multiplicative-noise misfit for pointwise observations.

class hippymfem.modeling.misfit.ContinuousStateObservation(Vh, attributes=None, bcs=None, data=None, noise_variance=None, boundary=False, coefficient=None)

Bases: Misfit, KeepAlive

Misfit in a weighted \(L^2\) norm over the domain or a subdomain: \(\tfrac{1}{2\sigma^2}(u-d)^{\!\top} W (u-d)\).

Parameters:

• Vh (FunctionSpace) – State space.

• attributes (sequence of int, optional) – Mesh element attributes defining the observed subdomain; the default is the whole domain.

• bcs (optional) – Essential conditions whose rows and columns are removed from W.

• boundary (bool) – Observe on the boundary (ds) rather than the interior (dx).

createVecLeft()

cost(x)

grad(i, x, out)

setLinearizationPoint(x, gauss_newton_approx=False)

block_is_zero(i, j)

Whether the second-derivative block (i, j) vanishes identically, so that a caller can skip applying and adding it. False unless a class knows better.

apply_ij(i, j, dir, out)

class hippymfem.modeling.misfit.MultiStateMisfit(misfits=None)

Bases: Misfit

Sum of misfits over several independent states (experiments).

append(misfit)

cost(x)

grad(i, x, out)

setLinearizationPoint(x, gauss_newton_approx=False)

apply_ij(i, j, dir, out)

class hippymfem.modeling.misfit.MisfitTD(misfits, sim_times)

Bases: Misfit

Time-dependent misfit: a list of stationary misfits, one per observation time.

cost(x)

grad(i, x, out)

setLinearizationPoint(x, gauss_newton_approx=False)

apply_ij(i, j, dir, out)

hippymfem.modeling.pointwiseObservation

Pointwise observation operators.

Mesh::FindPoints locates each target in the local partition; the shape functions of the containing element, evaluated at the reference point it returns, are the row of the observation matrix.

Ownership. A target on a shared face is found by more than one rank and a target outside the domain by none. The first case is resolved by giving the target to the lowest-numbered rank that found it, so no observation is counted twice; the second is by default an error, reported with the offending coordinates rather than silently dropped (a missing observation changes the inverse problem).

Ordering. Observation vectors are distributed by owning rank, so their local order is not the order of the targets array. That is invisible as long as data comes from the operator itself; PointwiseObservation.gather() and PointwiseObservation.scatter() convert to and from the original target order for I/O.

class hippymfem.modeling.pointwiseObservation.PointwiseObservation(Vh, targets, components=None, component=None, warn_missing=None, missing='raise')

Bases: Operator

Evaluation of a finite element field at a set of points.

Parameters:

• Vh (FunctionSpace)

• targets ((ntargets, sdim) array)

• components ((ntargets, vdim) array, optional) – Weights across the solution components, giving a line-of-sight observation \(\sum_c w_c u_c(x_t)\). Defaults to observing component 0 for scalar spaces, and is required for vector spaces unless component is given.

• component (int, optional) – Observe a single component of a vector space.

init_vector(x, dim)

mult(u, obs)

obs = B u: evaluate the field at the owned targets.

multTranspose(obs, u)

u = B^T obs: spread observation weights back to the dofs.

transpmult(obs, u)

u = B^T obs: spread observation weights back to the dofs.

getSize()

getComm()

createVecLeft()

createVecRight()

gather(obs)

Values at every target, in the order of targets, on every rank.

scatter(values, obs=None)

Fill an observation vector from an array in targets order.

noise(sigma=1.0, rng=None, out=None)

Observation noise keyed on the global target index.

Drawing in the partition-dependent local order would make synthetic data, and with it the MAP point, depend on the rank count. Targets are few, so every rank draws the whole global vector and keeps its own entries.

perturb(obs, sigma, rng=None)

Add partition-independent noise of standard deviation sigma.

observe(u, rel_noise=None, sigma=None, rng=None)

Synthetic data from a state: B u plus noise; returns (data, sigma).

The noise standard deviation is sigma, or rel_noise times the largest observed value; it is returned because the misfit needs it as the noise variance. rel_noise=0 or neither given returns clean data with sigma = 0.

hippymfem.modeling.pointwiseObservation.assemblePointwiseObservation(Vh, targets, **kw)

Observation operator evaluating the state at targets.

hippymfem.modeling.pointwiseObservation.assemblePointwiseLOSObservation(Vh, targets, directions)

Line-of-sight observation: the component of a vector field along a direction.

hippymfem.modeling.pointwiseObservation.exportPointwiseObservation(B, data, filename, comm=None)

Write targets and observed values as CSV on rank 0.

hippymfem.modeling.model

The inverse problem: forward model + prior + misfit.

The cost functional is

\[J(m) = \underbrace{\Phi(u(m))}_{\text{misfit}} + \underbrace{\tfrac12\|m-\bar m\|_R^2}_{\text{regularization}},\]

and Model assembles the adjoint-based gradient and the action of the reduced Hessian from the blocks the PDE problem and misfit provide. Signatures follow hIPPYlib exactly.

class hippymfem.modeling.model.ReducedMap(problem, functional)

Bases: object

A PDE problem and a functional of (u, m) reduced over the state.

Both the inverse problem (Model, whose functional is the data misfit) and the parameter-to-QoI map (Parameter2QoiMap, whose functional is the QoI) reduce a functional \(f(u(m), m)\) through the same adjoint machinery: the adjoint source is \(-\partial_u f\), the reduced gradient \(\partial_m f + C^{\!\top}p\), and the reduced Hessian action combines the PDE’s blocks with the functional’s. This class holds that machinery once; the subclasses add what differs (the prior and the cost for the inverse problem, the value for the QoI map).

The functional supplies grad(i, x, out), setLinearizationPoint and apply_ij(i, j, d, out); the problem is a PDEProblem.

generate_vector(component='ALL')

A vector (or the triple) in the shape of the problem’s variables.

init_parameter(m)

solveFwd(out, x)

solveAdj(out, x)

Solve the adjoint problem with right-hand side \(-\partial_u f\).

setLinearizationPoint(x, gauss_newton_approx=False, matrix_free=False)

Fix the linearization point of the reduced Hessian.

matrix_free is handed to the problem (see setLinearizationPoint()): its second-derivative products are then computed from the element kernels rather than from assembled blocks.

solveFwdIncremental(sol, rhs, rel_tolerance=None)

The incremental forward solve; rel_tolerance loosens the solver’s own tolerance for this solve (a problem whose solveIncremental takes it).

solveAdjIncremental(sol, rhs, rel_tolerance=None)

applyC(dm, out)

applyCt(dp, out)

applyWuu(du, out)

W_uu from the functional plus, unless Gauss-Newton, from the PDE.

applyWum(dm, out)

applyWmu(du, out)

applyWmm(dm, out)

apply_ij(i, j, d, out)

Dispatch to the corresponding applyXY.

class hippymfem.modeling.model.Model(problem, prior, misfit, gradient_norm='M')

Bases: ReducedMap

Full description of a PDE-constrained Bayesian inverse problem.

Parameters:

• problem (PDEProblem)

• prior (a prior from hippymfem.modeling.prior)

• misfit (a misfit from hippymfem.modeling.misfit)

gradient_norm

"M" (the default) gives the discrete L2 norm sqrt(g^T M^{-1} g) that hIPPYlib reports, so iteration histories, Newton-CG’s stopping test and its Eisenstat-Walker tolerances match hIPPYlib’s; "R" gives the prior-preconditioned norm sqrt(g^T R^{-1} g). Both are mesh-independent.

Type:

Riesz map used to measure the gradient

cost(x)

[total, regularization, misfit] at x = [u, m, p].

evalGradientParameter(x, mg, misfit_only=False)

Gradient of the cost in mg; returns its norm.

The gradient is a dual-space object, so its norm needs a Riesz map, set by gradient_norm (by default the discrete \(L^2\) norm \(\sqrt{g^{\!\top}M^{-1}g}\) that hIPPYlib reports).

setPointForHessianEvaluations(x, gauss_newton_approx=False, matrix_free=False)

Fix the linearization point of the reduced Hessian (matrix_free: see ReducedMap.setLinearizationPoint()).

applyR(dm, out)

Rsolver(tolerance=None, cycles=None)

The solver of the prior’s precision. With tolerance, as the preconditioner of a Krylov iteration: its solves stopped at that relative tolerance; with cycles, each of them a fixed polynomial of that many multigrid cycles (prior.getHessianPreconditioner).

hippymfem.modeling.reducedHessian

The reduced Hessian, applied matrix-free.

For a direction \(\tilde m\), the Gauss-Newton Hessian is

\[H_{GN}\tilde m = C^{\!\top} A^{-\top} W_{uu} A^{-1} C \tilde m + R \tilde m,\]

and the full Hessian adds the \(W_{um}\), \(W_{mu}\) and \(W_{mm}\) terms. The Gauss-Newton product needs no explicit sign: the incremental forward solve returns \(\hat u = A^{-1}C\tilde m = -\,\mathrm{d}u/\mathrm{d}m[\tilde m]\), and that sign cancels against the one in the adjoint source. Each application costs one incremental forward and one incremental adjoint solve, which is why those solves dominate Newton-CG.

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

Bases: Operator

Matrix-free reduced Hessian at the point set on the model.

Parameters:

• model (Model) – setPointForHessianEvaluations must have been called.

• misfit_only (bool) – Drop the R term, giving the Hessian of the data misfit alone (the operator whose dominant eigenvectors define the Laplace approximation).

accuracy

relative tolerance of the two incremental solves of an action, where looser than the solvers’ own (set_accuracy()); None: their own

init_vector(x, dim)

set_accuracy(tolerance)

Ask the next actions for incremental solves stopped at the relative tolerance tolerance, if that is looser than the solvers’ own (None: their own again). A Krylov iteration can use products that become less accurate as it converges (CGSolverSteihaug, relax_operator). The prior’s part of an action keeps its accuracy: a relative error of its mass solve is the same relative error of R x, with nothing to damp it, and at the accuracies that the incremental solves tolerate it changed the path of the model problem.

mult(x, y)

multTranspose(x, y)

inner(x, y)

y^T A x.

GNHessian(x, y)

TrueHessian(x, y)

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

Bases: Operator

Reduced Hessian by central finite differences of the gradient.

Slow, and only as accurate as the step size, but it needs nothing but solveFwd/solveAdj/evalGradientParameter, which makes it the right reference for checking ReducedHessian.

init_vector(x, dim)

mult(x, y)

multTranspose(x, y)

inner(x, y)

y^T A x.

hippymfem.modeling.posterior

Low-rank Gaussian (Laplace) approximation of the posterior.

If \((d_i, u_i)\) are the dominant generalized eigenpairs of

\[H_{\rm misfit}\, u_i = d_i\, R\, u_i, \qquad U^{\!\top} R U = I,\]

then the Laplace approximation at the MAP point has

• precision apply: \((R + R U D U^{\!\top} R)\,x\);

• covariance apply: \(\left(R^{-1} - U (I + D^{-1})^{-1} U^{\!\top}\right) x\);

• sampling: \(y = (I - U S U^{\!\top} R) x\) with \(S = I - (I+D)^{-1/2}\) and \(x \sim N(0, R^{-1})\).

Once the eigenpairs are in hand, all three need only prior operations and k inner products, with no further PDE solves.

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

Bases: Operator

Low-rank posterior precision and its inverse.

init_vector(x, dim)

inner(x, y)

y^T A x.

mult(x, y)

\(y = (R + R U D U^{\!\top} R) x\).

multTranspose(x, y)

\(y = (R + R U D U^{\!\top} R) x\).

solve(sol, rhs)

\(\mathrm{sol} = (R^{-1} - U (I + D^{-1})^{-1} U^{\!\top})\,\mathrm{rhs}\).

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

Bases: KeepAlive

Turn a prior sample into a posterior sample: \(y = (I - U S U^{\!\top} R)x\).

init_vector(x, dim)

sample(noise, s)

noise is a prior sample (zero mean); s becomes a posterior sample.

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

Bases: KeepAlive

Laplace approximation of the posterior at the MAP point.

Parameters:

• prior (a prior object)

• d (the dominant generalized eigenpairs of the data-misfit Hessian)

• U (the dominant generalized eigenpairs of the data-misfit Hessian)

• mean (ParVector, optional) – The MAP point.

init_vector(x, dim)

cost(m)

sample(*args, noise=None, s_prior=None, s_post=None, add_mean=True, rng=None)

Draw a posterior sample; returns s_post.

hIPPYlib’s two positional forms work: sample(noise, s_prior, s_post) with white noise, filling the prior and the posterior sample, and sample(s_prior, s_post) from an existing zero-mean prior sample. The keyword form makes the roles explicit, and every vector may be left out: post.sample() draws the noise (from rng or parRandom) and returns a new posterior sample, post.sample(s_prior=s) transforms a given prior sample.

trace(**kwargs)

(posterior, prior, correction) traces of the covariance.

trace_update()

Trace of the data-informed correction, \(\mathrm{tr}(U(I+D^{-1})^{-1}U^{\!\top}M)\).

pointwise_variance(**kwargs)

(posterior, prior, correction) pointwise variance fields.

klDistanceFromPrior(sub_comp=False)

Kullback-Leibler divergence of the Laplace posterior from the prior.

hippymfem.modeling.modelVerify

Finite-difference verification of the gradient and the Hessian.

Two checks, both from hIPPYlib:

• the cost should satisfy \(J(m+h\,\tilde m) = J(m) + h\,g^{\!\top}\tilde m + O(h^2)\), so the error in the first-order model falls like \(h\);

• the gradient should satisfy \(g(m+h\tilde m) = g(m) + h\,H\tilde m + O(h^2)\), likewise.

Seeing those slopes shows that the adjoint gradient and the second-order blocks are consistent with the cost.

hippymfem.modeling.modelVerify.modelVerify(model, m0, is_quadratic=False, misfit_only=False, verbose=True, eps=None, plotting=False, filename=None)

Run the finite-difference checks at m0.

Parameters:

• model (Model)

• m0 (ParVector) – The parameter at which the derivatives are checked.

• is_quadratic (bool) – The cost is quadratic in m, so the Hessian error should be round-off.

• misfit_only (bool) – Check the misfit part of the cost only.

• verbose (bool) – Print the table of errors and observed slopes on rank 0.

• eps (array_like, optional) – Finite-difference steps; halvings from 1 down to 1e-10 by default.

• plotting (bool) – Draw both error curves against a first-order reference. The figure is left open, so a notebook shows it inline.

• filename (str, optional) – Also save the figure there, and close it.

Returns:

• eps (ndarray) – Step sizes used.

• err_grad (ndarray) – \(|(J(m+h\tilde m) - J(m))/h - g^{\!\top}\tilde m|\).

• err_H (ndarray) – \(\|(g(m+h\tilde m) - g(m))/h - H \tilde m\|_\infty\).

hippymfem.modeling.modelVerify.fd_slopes(eps, err)

Observed convergence rates log2(err[i]/err[i+1]) between steps.

hippymfem.modeling.modelVerify.best_slope(eps, err, lo=1e-07, hi=0.01)

Median observed slope over the step sizes where round-off has not taken over.

NaN where an error at one of those step sizes is not a finite number. A slope beside an error of zero is left out, which is round-off; an error that is not a number is a cost or a gradient that is not one, and the slope of the other steps must not stand for it.

hippymfem.modeling.timeDependentVector

A vector with one spatial field per discrete time level.

Time-dependent inverse problems treat the whole trajectory as the state, so the state vector is a list of spatial vectors indexed by time. Lookup is by time value rather than index, with a tolerance, because the misfit and the time stepper refer to times, not to positions in an array.

class hippymfem.modeling.timeDependentVector.TimeDependentVector(times, tol=1e-10, comm=None)

Bases: KeepAlive

A trajectory: one ParVector per entry of times.

Parameters:

• times (array_like) – Discrete times, ascending.

• tol (float) – Tolerance for matching a time to a level.

• comm (mpi4py communicator)

initialize(space_or_vector)

Allocate a zero field at every time level.

copy(other=None)

A deep copy; with other given, copy that trajectory into self.

index(t)

Index of the level nearest t; raises if none is within tol.

view(t)

The spatial vector at time t (a reference, not a copy).

store(u, t)

Copy the spatial vector u into the level at time t.

retrieve(u, t)

Copy the level at time t into u.

zero()

set(alpha)

scale(alpha)

axpy(a, other)

assign(other)

duplicate()

inner(other)

Sum over time levels of the spatial inner products.

dot(other)

Sum over time levels of the spatial inner products.

element_wise_inner(other)

Per-time-level inner products, as an array.

norm(time_norm='linf', space_norm='l2')

Norm over space at each level, combined over time.

time_norm is "linf", "l2" or "l1".

matmul(mat, out)

Apply a spatial operator or matrix at every level.