Priors
The prior is a Gaussian measure on the parameter field whose precision is a differential operator, so sampling from it and applying \(R^{-1}\) are linear solves rather than dense factorizations, which is what makes the whole approach work at mesh-independent cost.
prior = hm.BiLaplacianPrior(Vm, gamma=0.1, delta=0.5, robin_bc=True)
\(R = A M^{-1} A\) with \(A = \gamma(-\nabla\cdot\Theta\nabla) + \delta I\),
giving a Matern-class covariance. robin_bc=True adds the boundary term that
removes the spurious variance inflation a homogeneous Neumann condition produces
near the boundary.
parameter |
effect |
|---|---|
|
pointwise variance scales as \(1/(4\pi\gamma\delta)\) in 2D; the correlation length as \(\sqrt{\gamma/\delta}\) |
|
an anisotropic tensor for the differential
operator, as a |
|
the prior mean, a
|
|
the boundary term above |
|
|
|
the mass matrix \(M\) by the rule at the vertices: diagonal (below) |
A lumped mass. Every product with \(R\), which the cost, the gradient and each
Hessian action take, solves a system with the mass matrix; with the consistent one that
is some seventy CG iterations, 50 to 80 ms on a GPU at two million parameter unknowns,
as much as a solve with \(R\) itself. lumped_mass=True integrates the mass
matrix by the rule at the vertices of an element, which makes it diagonal (on affine
elements the sums of the rows of the consistent one): the lumped mass of the SPDE
approach to Matern fields (Lindgren, Rue and Lindstrom, 2011). The square root for
sampling is by the same rule, so samples have the covariance \(R^{-1}\) of that
precision exactly. It is another discretization of the same prior, with another
matrix \(R\), which is why it is not the default; it is for a parameter space of
first order on a mesh of one kind of element. On the MI250X cards of Frontier, with the
Newton benchmark at 1.8 million parameter unknowns a GPU, the operations of two
Newton-CG steps took 18.3 s instead of 19.4 s on one GPU and 23.9 s instead of 24.9 s on
64.
Use BiLaplacianComputeCoefficients() to go the other
way, from a desired pointwise standard deviation and correlation length to
gamma and delta. Do this rather than guessing: a prior whose pointwise
standard deviation is 28 for a field of amplitude 1 gives a Hessian with condition
number \(5\times10^8\), and Newton-CG will not converge on it.
Sampling
noise = prior.noise_vector()
prior.sample_noise(1.0, noise) # white noise, partition independent
m = Vm.vector()
prior.sample(noise, m) # m ~ N(mean, R^{-1})
Sampling needs \(\sqrt{M}\), the square root of the mass matrix, which is built in closed form at quadrature points rather than by a factorization: the element mass matrix of a single quadrature point is rank one, so its square root is explicit, and assembling those gives a matrix with \(\sqrt{M}\sqrt{M}^{\top} = M\) exactly. The test suite checks that identity to round-off rather than to a tolerance.
The white noise behind a sample is the same field on any number of ranks, so prior samples agree across rank counts to the tolerance of the prior’s solve; see Parallel correctness.
Other priors
\(R = A\), a rougher field |
|
Matern, the usual choice |
|
a vector-valued parameter |
|
a finite-dimensional Gaussian |
VectorBiLaplacianPrior applies one scalar gamma/delta pair to all
components, because MFEM’s vector integrators take one scalar coefficient;
component-dependent values raise rather than silently averaging.