Algorithms

Optimizers, linear solvers, randomized linear algebra.

hippymfem.algorithms.NewtonCG

Inexact Newton-CG in the reduced parameter space.

At each outer iteration the Newton system \(H \hat m = -g\) is solved inexactly by preconditioned CG (preconditioner \(R\)), truncated by an Eisenstat-Walker tolerance \(\min(\tau_0, \sqrt{\|g\|/\|g_0\|})\). The first GN_iter iterations use the Gauss-Newton Hessian, which is positive definite far from the minimum; the full Hessian takes over once the iterate is close enough for its indefiniteness not to matter.

Globalization is either an Armijo line search or a trust region; both are ported from hIPPYlib, with the same parameters, the same termination codes and the same printed table.

One default departs from hIPPYlib: cg_reorthogonalize. The CG iteration of a Newton step keeps its residuals orthogonal explicitly instead of by its recurrence (cgsolverSteihaug). On the model problem of the benchmarks with 2.1 million state dofs this took 131 CG iterations for twelve Newton steps where the recurrence took 193 to 210, a number that changed from run to run and from one GPU to another, and the 131 did not change when the incremental solves of the Hessian action were stopped at a relative residual of 1e-6 instead of 1e-12; with the recurrence, incremental solves stopped at 1e-8 already cost a third more CG iterations. The gradient stays exact (its forward and adjoint solves keep their tolerance), so the MAP point is the same; only the Newton directions come from a Hessian of lower accuracy. cg_reorthogonalize = False gives hIPPYlib’s iteration, for a comparison step by step.

Three parameters make the solves inside that CG cheaper. cg_preconditioner_tolerance (1e-6 by default) stops the prior’s solves where they precondition it; 0 leaves them at their own tolerance, as hIPPYlib does. cg_preconditioner_cycles (off by default) replaces each of those solves by a fixed polynomial of that many multigrid cycles, which is a linear map and needs no tolerance. cg_hessian_relaxation (off by default) lets the incremental solves of a Hessian action lose accuracy as the CG converges. They are described, with what they save and why they cannot go further, in the guide (docs/source/guide/optimization.rst).

exception hippymfem.algorithms.NewtonCG.ModelConvergenceError

Bases: RuntimeError

Raised by a forward solve that fails, so the line search can back off.

hippymfem.algorithms.NewtonCG.LS_ParameterList()

hippymfem.algorithms.NewtonCG.TR_ParameterList()

hippymfem.algorithms.NewtonCG.ReducedSpaceNewtonCG_ParameterList()

class hippymfem.algorithms.NewtonCG.ReducedSpaceNewtonCG(model, parameters=None, callback=None)

Bases: object

Inexact Newton-CG for the reduced (parameter-space) problem.

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', 'Forward solve failed during backtracking']

initial_grad_norm

gradient norm at the starting point, so callers can judge the reduction achieved independently of whether a requested tolerance was met: with inexact forward solves the attainable floor is set by the linear solver, not by the optimizer

grad_reduction

||g||/||g_0|| when a line search exhausted its backtracks

property comm

solve(x)

Minimize the cost from the initial guess x = [u, m, p].

hippymfem.algorithms.bfgs

(Limited-memory) BFGS for the reduced problem.

An alternative to Newton-CG when second-order information is unavailable or too expensive: it needs only gradients, but typically many more iterations. The inverse-Hessian approximation is applied by the two-loop recursion, with Powell damping so that the curvature condition \(s^{\!\top}y > 0\) holds even when the line search returns a poor step.

H0inv is the initial inverse-Hessian guess. Using the prior solver \(R^{-1}\) there, rather than a scaled identity, makes the iteration mesh-independent, for the same reason it preconditions Newton-CG.

hippymfem.algorithms.bfgs.BFGSoperator_ParameterList()

hippymfem.algorithms.bfgs.BFGS_ParameterList()

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

Bases: object

Default H0inv: multiplication by a scalar \(d_0\).

The scalar is refreshed from the oldest stored secant pair, which is the standard scaling that keeps the first steps from being wildly wrong.

init_vector(x, dim)

solve(x, b)

class hippymfem.algorithms.bfgs.BFGS_operator(parameters=None)

Bases: object

The BFGS inverse-Hessian approximation, applied by the two-loop recursion.

set_H0inv(H0inv)

solve(x, b)

x = H_k b, the current approximation to \(H^{-1}b\).

update(s, y)

Add the secant pair (s, y), damping it if curvature is too small.

The damping needs \(H y\), written into help. hIPPYlib’s version (and this one in 0.1.0) computed it with solve() working in help as well, so H0inv.solve got its input as its output: the default rescaled identity zeroed it, \(y^{\!\top} H y\) came out 0, and a pair that needed damping was damped to s = 0 and raised.

class hippymfem.algorithms.bfgs.BFGS(model, parameters=None)

Bases: object

BFGS (or L-BFGS) for the reduced inverse problem.

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']

solve(x, H0inv=None, bounds_xPARAM=None)

Minimize the cost from the initial guess x.

H0inv defaults to a rescaled identity. bounds_xPARAM, if given, is a (lower, upper) pair of ParVectors, clipped after each step.

hippymfem.algorithms.steepestDescent

Steepest descent in the prior-preconditioned metric.

The descent direction is \(-R^{-1}g\), not \(-g\): gradients live in the dual space, and using the Euclidean gradient on a mesh-dependent discretization gives steps whose size changes with refinement.

hippymfem.algorithms.steepestDescent.SteepestDescent_ParameterList()

class hippymfem.algorithms.steepestDescent.SteepestDescent(model, parameters=None, callback=None)

Bases: object

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

solve(x)

hippymfem.algorithms.cgsolverSteihaug

Preconditioned CG with the Steihaug trust-region stopping rule.

Ported from hIPPYlib, including the termination codes, so that iteration counts and stopping behaviour are directly comparable between the two libraries, with one departure: with a trust region, a direction of nonpositive curvature is followed to the boundary, as Steihaug’s method prescribes. hIPPYlib takes the full first direction wherever that lands (outside the region when the radius is small) and, at a later iteration, stops at the current iterate. Without a trust region the two agree: the first direction in full, or the current iterate.

Termination reasons:

0. the iteration limit was reached (no convergence);

1. the preconditioned residual met the tolerance (convergence);

2. a direction of nonpositive curvature was found (the operator is not SPD);

3. the trust-region boundary was reached.

The stopping test is on the \(B^{-1}\)-norm of the residual, relative to its initial value or against an absolute floor.

A second departure is optional (reorthogonalize). The recurrence of CG keeps a residual orthogonal to the earlier ones only for one symmetric operator in exact arithmetic. A prior-preconditioned Hessian has a few large, well separated eigenvalues, for which that orthogonality is lost to rounding within a few tens of iterations, and a Hessian action computed with inexact solves is not one symmetric operator at all: in both cases the iteration is delayed, and by an amount that changes with the last digit of the operands. With reorthogonalize every new residual is made orthogonal, in the inner product of the preconditioner, to all the earlier ones explicitly, which costs k inner products and updates at iteration k and two stored vectors per iteration. The iteration then takes the number of steps of exact arithmetic, the same in every run, and keeps taking it when the operator carries an error far above rounding. ReducedSpaceNewtonCG turns it on (cg_reorthogonalize); here it is off by default, since a solve of a thousand iterations would store two thousand vectors.

hippymfem.algorithms.cgsolverSteihaug.CGSolverSteihaug_ParameterList()

class hippymfem.algorithms.cgsolverSteihaug.CGSolverSteihaug(parameters=None, comm=None)

Bases: object

Solve \(Ax=b\) with preconditioner \(B\) and Steihaug truncation.

A must provide mult(x, y) and init_vector(x, dim); B_solver must provide solve(x, b).

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

set_operator(A)

set_preconditioner(B_solver)

set_TR(radius, B_op)

Restrict steps to the ball x^T B_op x <= radius^2.

update_x_without_TR(x, alpha, d)

update_x_with_TR(x, alpha, d)

Step, and if that leaves the trust region, stop on its boundary.

solve(x, b)

Solve; returns the iteration count. With relax_operator the accuracy asked of the operator is put back to its own when the solve ends.

hippymfem.algorithms.linSolvers

Linear solvers.

hIPPYlib leans on PETSc here, and in particular on PETScLUSolver for the forward, adjoint and incremental solves. PyMFEM’s standard parallel build links no direct solver (no MUMPS, SuperLU_dist, STRUMPACK or PETSc), so the exact-solve capability is supplied from the Python side in directSolvers:

hIPPYlib	hIPPyMFEM
PETScKrylovSolver	KrylovSolver (MFEM Krylov + hypre PC), or PETScKrylovSolver when petsc4py is importable
PETScLUSolver	LUSolver, exact on any rank count, by replicating the matrix and factorizing with scipy’s SuperLU; or PETScLUSolver, which uses a genuinely distributed factorization when the local PETSc has one
amg_method()	hypre BoomerAMG

In parallel the factorization is replicated on every rank, so LUSolver refuses a matrix larger than max_global_size rather than silently exhausting memory.

hippymfem.algorithms.linSolvers.KrylovSolver_ParameterList()

hippymfem.algorithms.linSolvers.HYPRE_PCG = True

Whether a CG solve with a hypre preconditioner runs in hypre’s own PCG (MFEM’s HyprePCG) instead of MFEM’s CGSolver (HIPPYMFEM_HYPRE_PCG, on by default). The two make the same iteration: the same iterates, counts and stopping test from a zero initial guess. hypre’s PCG tells BoomerAMG that the initial guess of a V-cycle is zero, which saves the first relaxation its matrix-vector product on the finest level; MFEM zeroes the vector without saying so. With l1-Jacobi relaxation on a GPU that is one product of four per iteration: a solve of 2.1 million dofs to 1e-12 took 0.103 s and 0.089 s on an H100, 0.274 s and 0.227 s on an L40S (24 iterations each).

hippymfem.algorithms.linSolvers.set_hypre_pcg(on=True)

Choose hypre’s own PCG (True) or MFEM’s CG (False) for the CG solves with a hypre preconditioner that are set up from now on (see HYPRE_PCG); returns the old choice.

class hippymfem.algorithms.linSolvers.KrylovSolver(comm, method='cg', precond='amg', systems_dim=None, elasticity=False, fes=None, systems_bynodes=None)

Bases: _SolverBase

MFEM Krylov solver with a hypre preconditioner.

Parameters:

• comm (mpi4py communicator)

• method ({"cg", "gmres", "fgmres", "bicgstab", "minres"})

• precond ({"amg", "ilu", "jacobi", "l1jacobi", "gs", "parasails", "euclid", "none"})

• systems_dim (int, optional) – Pass the number of components for BoomerAMG’s systems-of-PDEs mode; use this for vector-valued (elasticity-like) operators, where scalar AMG converges poorly.

• elasticity (bool) – Turn on BoomerAMG’s elasticity interpolation (needs systems_dim). With hypre on a GPU, where MFEM has none, the systems options alone, with a warning.

• systems_bynodes (bool, optional) – How the unknowns of the components are ordered, which BoomerAMG has to be told to pair the rows of one component: component after component (mfem.Ordering.byNODES, the ordering of this library’s spaces unless another is asked for, and the default here), or interleaved (byVDIM: pass False). Taken from fes where that is given. With systems_dim and neither of the two, the default is taken and a warning says so once: on an interleaved space BoomerAMG then pairs rows of different components, and a solve takes 15 to 50 times the iterations.

• fes (FunctionSpace or mfem.ParFiniteElementSpace, optional) – The space of the unknowns, for their ordering and for the elasticity options.

ERROR_MODE

How hypre setup failures are handled. MFEM’s default aborts the whole job via MFEM_VERIFY, which is wrong inside an optimizer: a line search routinely proposes a parameter whose operator defeats AMG’s coarsening (exp(m) spanning hundreds of decades), and that step should be rejected, not the run killed. Warning lets the solve fail to converge and surface as a RuntimeError, which the line search already catches.

set_operator(A, pc=None, donor=None)

Point the solver at A; pc shares an existing preconditioner of A.

The forward solver and the forward-incremental solver hold the same Jacobian; the PDE passes the forward solver’s preconditioner here when the operators are the same object, so the hierarchy is built once and held once. donor is the solver pc comes from: when that one runs hypre’s PCG on this very matrix (HYPRE_PCG), its PCG object is shared too, since each such object sets its preconditioner up at its first solve.

A SingleParMatrix is solved in the single-precision hypre it lives in (CG with BoomerAMG only); pc is then the engine of another solver of the same matrix.

SetOperator(A, pc=None, donor=None)

Point the solver at A; pc shares an existing preconditioner of A.

The forward solver and the forward-incremental solver hold the same Jacobian; the PDE passes the forward solver’s preconditioner here when the operators are the same object, so the hierarchy is built once and held once. donor is the solver pc comes from: when that one runs hypre’s PCG on this very matrix (HYPRE_PCG), its PCG object is shared too, since each such object sets its preconditioner up at its first solve.

A SingleParMatrix is solved in the single-precision hypre it lives in (CG with BoomerAMG only); pc is then the engine of another solver of the same matrix.

property tolerance_floor

The relative residual below which a solve with the present operator does not go, whatever tolerance is asked: that of single precision for a matrix of the single-precision hypre, zero otherwise.

can_solve_single()

Whether this solver could run in the single-precision hypre: CG or GMRES with a plain BoomerAMG, the combinations set up there (GMRES for a Jacobian that is not symmetric, and for its transpose).

solve(x, b)

x = A^{-1} b; returns the iteration count.

class hippymfem.algorithms.linSolvers.ChebyshevSolver(krylov, degree=4, lanczos=12, safety=(0.85, 1.05))

Bases: KeepAlive

A fixed polynomial in place of a Krylov solve: \(x = q(CA)\,C\,b\).

krylov is a KrylovSolver whose operator \(A\) is set, with its preconditioner \(C\) (a multigrid cycle). degree steps of the Chebyshev iteration on the interval of the eigenvalues of \(CA\) apply \(C\) that many times (a BoomerAMG cycle from the iterate before it, any other preconditioner to a residual updated with a product of \(A\)), with coefficients that do not depend on the right-hand side. The result is one linear map, symmetric and positive definite where \(C\) is, and its error in the energy norm is at most

\[1 / T_{\mathrm{degree}}\Bigl(\frac{\lambda_{\max} + \lambda_{\min}} {\lambda_{\max} - \lambda_{\min}}\Bigr)\]

of the solution (error_bound): 1.6 % with four cycles of a multigrid with which CG gains a factor of 0.3 an iteration, where CG needs twelve iterations for six digits.

That is what a preconditioner has to be. A Krylov solve stopped at a tolerance is a different map for every right-hand side, so it has to be stopped far below the tolerance of the iteration it preconditions (ReducedSpaceNewtonCG, cg_preconditioner_tolerance); a fixed polynomial needs no tolerance, and the iteration it preconditions sees an operator whose spectrum is that of the exact preconditioner moved by at most twice the bound. It does not replace the solve where the solution itself is wanted (a sample of the prior, a variance).

The interval is estimated once, at the first solve, from the Lanczos matrix of lanczos CG iterations on a vector of random signs; the Ritz values lie inside the spectrum, so the interval used is wider (safety). Outside the interval the polynomial’s error grows slowly and stays below one up to twice its centre, so an estimate that is somewhat off costs accuracy, not definiteness.

interval

the interval of the polynomial, and the Ritz values it came from

ritz

the interval of the polynomial, and the Ritz values it came from

static fits(krylov)

Whether krylov can be replaced: CG with a preconditioner, in double.

init_vector(v, dim)

property error_bound

The bound on the relative error of a solve in the energy norm, for eigenvalues inside interval (None before the first solve).

setup()

Estimate the interval (and let MFEM set the preconditioner up for a direct application, which it does at the first one whatever hypre’s PCG did before).

solve(x, b)

x = q(C A) C b; returns the number of cycles.

A step of the iteration needs C r for the residual r of the current iterate y. With BoomerAMG a cycle applied from y (MFEM’s iterative_mode) returns y + C (b - A y), the residual formed inside it by the relaxation that a cycle begins with: the step is the difference, and no product of the operator is made here. (Applied from zero, a cycle begins with that same product, of a zero vector, since MFEM does not mark the vector as zero for hypre: one product with its exchange more in every cycle.) Any other preconditioner is applied to a residual that is updated with a product, since not every one of them starts from the vector it is given.

class hippymfem.algorithms.linSolvers.CycleOf(cycle, krylov)

Bases: object

A Krylov solver’s operator with another preconditioner than its own, as ChebyshevSolver reads a solver: cycle (Mult(b, x) on MFEM’s vectors, from zero) in place of the preconditioner of krylov.

method = 'cg'

precond_type = 'cycle'

the polynomial then applies the cycle to a residual that it updates)

Type:

(not “amg”

property A

init_vector(v, dim)

class hippymfem.algorithms.linSolvers.LUSolver(comm=None, method='superlu', max_global_size=None, replicate=False)

Bases: _SolverBase

Sparse direct solve via scipy’s SuperLU, on any number of ranks.

Stands in for hIPPYlib’s PETScLUSolver. On one rank this factorizes the matrix in place. In parallel it gathers the matrix onto rank 0 and factorizes it there (on every rank with replicate=True), which is exact and independent of the rank count; see ReplicatedLUSolver, which this delegates to, for the memory trade and the size guard.

Parameters:

• comm (mpi4py communicator)

• method (str) – Kept for source compatibility; "superlu" is the only factorization.

• max_global_size (int, optional) – Passed through to the replicated solver.

• replicate (bool) – Passed through: the factorization on every rank rather than on rank 0.

set_operator(A)

SetOperator(A)

solve(x, b)

solveTranspose(x, b)

x = A^{-T} b.

mult(x, y)

y = A^{-1} x, so the solver can be used as an operator.

class hippymfem.algorithms.linSolvers.TransposeSolver(solver_factory)

Bases: KeepAlive

Solve with A^T by forming the transpose once and wrapping a solver.

release()

As KrylovSolver.release(): drop the transpose and the inner solver.

set_operator(A)

SetOperator(A)

property parameters

solve(x, b)

init_vector(v, dim)

class hippymfem.algorithms.linSolvers.LumpedMassSolver(M, comm=None)

Bases: KeepAlive

Inverse of the row-sum-lumped diagonal of a mass matrix.

A cheap, robust stand-in for an exact mass solve when only a spectrally equivalent operator is needed.

set_operator(M)

solve(x, b)

mult(x, y)

init_vector(v, dim)

hippymfem.algorithms.linSolvers.auto_solver(space_or_size, comm=None, max_direct=None, method='cg', precond='amg', **kw)

A direct solve where one is affordable, Krylov where it is not.

An exact solve makes a Hessian spectrum a property of the discretization rather than of a Krylov tolerance, and it is rank-count independent. It is the wrong default at scale, because the replicated factorization (the only one available without petsc4py) costs O(nnz(L+U)) on every rank. The choice therefore follows the problem size, not the rank count: a one-rank run on a fine mesh is exactly where a direct solve fails.

Parameters:

• space_or_size (FunctionSpace or int) – The space to be solved on, or its global number of unknowns.

• max_direct (int, optional) – Largest problem to factorize. Defaults to the replicated solver’s own refusal threshold, so this never selects a solver that would then decline.

• preferred (Where petsc4py is available its distributed factorization is)

• since

• fill. (it is exact without replicating the)

hippymfem.algorithms.linSolvers.make_solver(comm, kind='krylov', method='cg', precond='amg', **kw)

Factory used by the priors and PDE problems.

kind is one of

"krylov" / "iterative"

MFEM Krylov with a hypre preconditioner.

"lu" / "direct"

Exact scipy SuperLU, replicated in parallel (LUSolver).

"petsc_lu"

PETSc’s best available distributed factorization.

"petsc"

PETSc Krylov.

"auto"

auto_solver(); needs space_or_size.

The PETSc kinds need petsc4py to have been imported before PyMFEM; see directSolvers.

hippymfem.algorithms.directSolvers

Parallel direct (exact) linear solvers.

hIPPYlib uses PETScLUSolver, a distributed MUMPS factorization, whenever an exact forward, adjoint or incremental solve is wanted: an exact Hessian action lets a randomized eigensolver or a Newton-CG tolerance study measure the discretization rather than the Krylov tolerance. PyMFEM’s standard parallel build links none of MUMPS, SuperLU_dist, STRUMPACK or PETSc, so this module supplies the capability from the Python side, by two exact routes:

ReplicatedLUSolver

Gathers the distributed matrix onto rank 0 and factorizes it there with scipy’s SuperLU; a solve gathers the right-hand side, solves on rank 0 and scatters the result. Needs nothing beyond scipy, works on any rank count, and matches a distributed direct solve to round-off. It costs one serial factorization of the whole matrix in time and O(nnz(L+U)) on rank 0 in memory, the standard “redundant coarse solve” trade and the right one up to a few hundred thousand unknowns. replicate=True factorizes on every rank instead, which the name comes from: every rank then solves for itself and no rank waits on rank 0, for O(nnz(L+U)) on each. Beyond that size use a Krylov solver.

PETScLUSolver

A distributed factorization when petsc4py is importable, picking the best package the local PETSc was built with (MUMPS, SuperLU_dist, STRUMPACK, PaStiX) and falling back to PETSc’s own redundant + lu, which is exact as well. PETScKrylovSolver exposes PETSc’s Krylov methods and preconditioners for the same matrices.

Import order matters for the PETSc route. petsc4py and PyMFEM are linked against MPI libraries that must be resolved in that order; importing petsc4py after mfem raises ImportError: libmpi_mpifh.so.40: undefined symbol: mpi_conversion_fn_null_. Set HIPPYMFEM_PETSC=1 before importing hippymfem (which then imports petsc4py first), or import petsc4py yourself at the top of the script. petsc_available() reports what happened instead of guessing.

class hippymfem.algorithms.directSolvers.ReplicatedLUSolver(comm=None, max_global_size=None, permc_spec=None, replicate=False)

Bases: _Base

Exact sparse LU of the gathered matrix, on rank 0 or on every rank.

A drop-in stand-in for hIPPYlib’s PETScLUSolver that needs no external package. By default the matrix is gathered onto rank 0 and factorized there; solve gathers the right-hand side, solves on rank 0 and scatters the result, so the other ranks hold no factorization and a node’s memory is charged once, not once per rank. With replicate=True every rank holds the factorization and solves for itself from an all-gathered right-hand side. Either way the result is independent of the rank count and the same on every rank.

Parameters:

• comm (mpi4py communicator)

• max_global_size (int) – Refuse matrices larger than this. The default keeps a typo from turning into an out-of-memory kill; raise it deliberately if the fill fits.

• replicate (bool) – Factorize on every rank rather than on rank 0 alone.

set_operator(A)

SetOperator(A)

solve(x, b)

x = A^{-1} b.

solveTranspose(x, b)

x = A^{-T} b.

mult(x, y)

y = A^{-1} x, so the solver can be used as an operator.

class hippymfem.algorithms.directSolvers.PETScLUSolver(comm=None, package='auto')

Bases: _PETScBase

Distributed exact LU through PETSc.

package="auto" tries the parallel factorizations the local PETSc was built with and, if none is available, falls back to redundant with an inner lu: every rank factorizes the gathered matrix, which is still an exact solve. package_used records what was chosen, so a script can report it rather than assume it.

set_operator(A)

SetOperator(A)

solveTranspose(x, b)

x = A^{-T} b.

class hippymfem.algorithms.directSolvers.PETScKrylovSolver(comm=None, method='cg', precond='gamg', **options)

Bases: _PETScBase

PETSc Krylov solver, the counterpart of hIPPYlib’s PETScKrylovSolver.

method and precond are PETSc names ("cg", "gmres"; "gamg", "hypre", "ilu", "jacobi", "none"). Extra PETSc options can be passed as keywords, e.g. pc_hypre_type="boomeramg".

set_operator(A)

SetOperator(A)

hippymfem.algorithms.directSolvers.petsc_available()

(True, "") when the PETSc route can be used, else (False, reason).

hippymfem.algorithms.directSolvers.preload_petsc()

Import petsc4py and return the PETSc module, or None on failure.

hippymfem.algorithms.directSolvers.gather_matrix(A, comm=None, format='csc', root=None)

Assemble a distributed HypreParMatrix into a scipy matrix on every rank, or on root alone (None on the other ranks) when a root is given.

Blocks are ordered by their true row offset rather than by rank, so this does not assume that hypre’s partition runs in rank order.

hippymfem.algorithms.multivector

A collection of distributed vectors with BLAS-3 reductions.

hIPPYlib stores a MultiVector as a Python list of vectors, so V^T W costs nvec^2 separate global reductions. Here the columns are rows of a single (nvec, n_local) C-contiguous array, so V^T W is one local GEMM plus one Allreduce. Each column is still exposed as a ParVector aliasing its row, so code written against the hIPPYlib interface works unchanged.

class hippymfem.algorithms.multivector.MultiVector(v=None, nvec=None)

Bases: object

nvec distributed vectors sharing one layout.

Parameters:

• v (ParVector or MultiVector, optional) – Layout template, or a MultiVector to deep-copy.

• nvec (int, optional) – Number of columns; required when v is a ParVector.

property data

The (nvec, n_local) backing array, current on the host.

classmethod from_array(comm, data)

Adopt a (nvec, n_local) array (no copy).

nvec()

setSizeFromVector(v, nvec)

Reset to nvec zero columns with v’s layout.

copy()

dot(v)

self^T v.

Returns a 1-D array of length nvec for a vector argument and a flattened nvec x v.nvec array for a MultiVector argument, which is hIPPYlib’s convention; dot_mv() reshapes it.

dot_v(v)

dot_mv(mv)

norm(norm_type='l2')

Per-column norms as a 1-D array.

reduce(v, alpha)

v += sum_i alpha[i] * self[i].

axpy(a, y)

self[i] += a[i]*y[i] for a MultiVector y, or self[i] += a*y.

scale(k, a=None)

Scale every column by k[i], or column k by a.

zero()

swap(other)

orthogonalize()

Euclidean QR in place; self becomes Q, returns R.

Borthogonalize(B)

B-orthogonal QR in place.

Returns (Bq, r) where Bq holds B @ Q and r is the triangular factor, with Q^T B Q = I. self becomes Q.

hippymfem.algorithms.multivector.MatMvMult(A, x, y, ensemble=None)

y[i] = A x[i] for every column.

ensemble is a communicator whose ranks each hold a complete copy of the problem (operators and vectors on MPI.COMM_SELF, or on a communicator of their own). The columns are then divided among its ranks, rank r applies A to the columns r, r + size, ..., and the results are exchanged so that every rank ends with all of y. The applications are independent, so this scales with the number of ranks, which a domain decomposition of a small problem over the same GPUs does not (see ensemble_columns()).

hippymfem.algorithms.multivector.ensemble_columns(ensemble, nvec, rank=None)

The columns of an nvec-column MultiVector that a rank of ensemble computes.

hippymfem.algorithms.multivector.ensemble_exchange(ensemble, y)

Give every rank of ensemble all columns of y, each taken from the rank that computed it (ensemble_columns()). The copies of y must have the same local size on every rank, which they have when each rank holds the whole problem.

hippymfem.algorithms.multivector.MatMvTranspmult(A, x, y)

y[i] = A^T x[i] for every column.

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

Y = X A for a small dense A (X.nvec x Y.nvec).

hippymfem.algorithms.randomizedEigensolver

Randomized eigensolvers for the (generalized) Hermitian eigenvalue problem.

The Laplace approximation needs the dominant eigenpairs of the prior-preconditioned data misfit Hessian, i.e. \(H_{\rm misfit} u = \lambda R u\). The dominant spectrum decays fast for ill-posed inverse problems, which is exactly the regime where randomized methods are appropriate.

References

1. N. Halko, P.G. Martinsson, J.A. Tropp, Finding structure with randomness, SIAM Review 53(2), 2011.

2. A.K. Saibaba, J. Lee, P.K. Kitanidis, Randomized algorithms for generalized Hermitian eigenvalue problems, NLAA 23(2), 2016.

hippymfem.algorithms.randomizedEigensolver.singlePass(A, Omega, k, s=1, check=False, ensemble=None)

Dominant eigenpairs of A from one pass over the operator.

Returns (d, U) with \(U^{\!\top}U = I_k\). ensemble: see doublePassG().

hippymfem.algorithms.randomizedEigensolver.doublePass(A, Omega, k, s=1, check=False, ensemble=None)

Dominant eigenpairs of A with s power iterations and two passes.

More accurate than singlePass() at the cost of k extra operator applications; this is the default for the Laplace approximation. ensemble: see doublePassG().

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

Single-pass solver for \(A u = \lambda B u\), with \(U^{\!\top}BU=I\).

ensemble: see doublePassG().

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

Double-pass solver for \(A u = \lambda B u\), with \(U^{\!\top}BU=I\).

A is the data-misfit Hessian, B the prior precision R and Binv a solver for it.

ensemble is a communicator whose ranks each hold a complete copy of the problem, with the same Omega; the applications of A and of Binv to the columns are then divided among its ranks (MatMvMult()), and every rank returns the same eigenpairs. The B-orthogonalization, which is sequential in the columns, is carried out by every rank.

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

Orthogonality and residual diagnostics for a standard eigenproblem.

Returns (err_AU_UD, err_UtU_I, err_UtAU_D); all three should be small.

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

Diagnostics for a generalized eigenproblem: residual and B-orthogonality.

hippymfem.algorithms.randomizedSVD

Randomized SVD.

Needed for non-symmetric reduced operators, above all the parameter-to-observable map \(B A^{-1} C\) of a linear inverse problem, whose singular values say how much of the parameter the data can see.

hippymfem.algorithms.randomizedSVD.accuracyEnhancedSVD(A, Omega, k, s=1, check=False)

Randomized SVD of a (possibly non-symmetric) operator.

Returns (U, sigma, V) with U in the range, V in the domain, and \(A \approx U \,\mathrm{diag}(\sigma)\, V^{\!\top}\).

s power iterations are applied with re-orthogonalization at every step, which is what keeps the small singular values from being swamped.

hippymfem.algorithms.randomizedSVD.singlePassSVD(A, Omega, k, check=False)

Randomized SVD: accuracyEnhancedSVD() with s=1, its default.

hippymfem.algorithms.randomizedSVD.check_SVD(A, U, sigma, V)

Diagnostics: (||A v_i - s_i u_i||, ||U^T U - I||, ||V^T V - I||).

hippymfem.algorithms.lowRankOperator

Low-rank operators built from a spectral decomposition.

Represents \(U D U^{\!\top}\) (optionally in a \(B\) inner product, where the columns of U are \(B\)-orthonormal). Used by the Laplace approximation, where d and U come from the randomized generalized eigensolver.

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

Bases: Operator

\(A = U\,\mathrm{diag}(d)\,U^{\!\top} B\).

Parameters:

• d (ndarray)

• U (MultiVector) – Columns \(B\)-orthonormal: \(U^{\!\top}BU = I\).

• my_init_vector (callable, optional) – Supplies the layout, when U’s own is not what the caller wants.

init_vector(x, dim)

mult(x, y)

y = U diag(d) U^T x.

multTranspose(x, y)

y = U diag(d) U^T x.

inner(x, y)

y^T A x.

solve(sol, rhs)

sol = U diag(1/d) U^T rhs (the pseudo-inverse on the range of U).

get_diagonal(diag)

Diagonal of \(U \mathrm{diag}(d) U^{\!\top}\).

trace(W=None)

Trace of \(U D U^{\!\top} W\); W = None means the identity.

hippymfem.algorithms.traceEstimator

Hutchinson trace estimation with a running error estimate.

Used for the integrated prior variance when the exact diagonal would cost one solve per degree of freedom.

class hippymfem.algorithms.traceEstimator.TraceEstimator(A, solve_mode=False, accurancy=0.1, init_vector=None, random_engine=None, accuracy=None)

Bases: object

Estimate \(\mathrm{tr}(A)\) by averaging \(z^{\!\top}Az\) over Rademacher vectors z.

Parameters:

• A (operator)

• solve_mode (bool) – Use A.solve instead of A.mult (so the trace of the inverse).

• accuracy (float) – Target relative standard error; iteration stops once reached. (accurancy, hIPPYlib’s spelling, is accepted too.)

• init_vector (callable, optional)

property accurancy

hippymfem.algorithms.cgsampler

Sampling from a Gaussian using conjugate gradients.

Draws \(x \sim N(0, A^{-1})\) by accumulating the Lanczos/CG search directions, which avoids ever forming a square root of \(A^{-1}\): if \(d_i\) are the CG directions for a random right-hand side, then \(\sum_i (z_i / \sqrt{d_i^{\!\top} A d_i})\, d_i\) has covariance \(A^{-1}\) restricted to the Krylov subspace explored.

Reference: Parker & Fox, Sampling Gaussian distributions in Krylov spaces with conjugate gradients, SIAM J. Sci. Comput. 34(3), 2012.

hippymfem.algorithms.cgsampler.CGSampler_ParameterList()

class hippymfem.algorithms.cgsampler.CGSampler(parameters=None)

Bases: object

Sample \(N(0, A^{-1})\) with conjugate gradients.

set_operator(A)

sample(noise, s)

s receives a sample; noise supplies the standard normals.

noise must have at least as many entries as CG takes iterations; an exhausted supply raises rather than reusing numbers. A numpy array must be the same on every rank. A ParVector contributes its global entries, in global order, so every rank uses the same variate at each iteration.