Numerical Gradient Computation for Simultaneous Detection of Geometry and Spatial Random Fields in a Statistical Framework

Abstract

The target of this chapter is the evaluation of gradients in inverse problems where spatial field parameters and geometry parameters are treated separately. Such an approach can be beneficial especially when the geometry needs to be detected accurately using L2-norm-based regularization. Emphasis is laid upon the computation of the gradients directly from the governing equations. Working in a statistical framework, the Karhunen-Loève (K-L) expansion is used for discretization of the spatial random field and inversion is done using the gradient-based Hamiltonian Monte Carlo (HMC) algorithm. The HMC gradients involve sensitivities w.r.t the random spatial field and geometry parameters. Building on a method developed by the authors, a procedure is developed which considers the gradients of the associated integral eigenvalue problem (IEVP) as well as the interaction between the gradients w.r.t random spatial field parameters and the gradients w.r.t the geometry parameters. The same mesh and linear shape functions are used in the finite element method employed to solve the forward problem, the artificial elastic deformation problem and the IEVP. Analysis of the rate of convergence using seven different meshes of increasing density indicates a linear rate of convergence of the gradients of the log posterior.

1. Introduction

Accurate computation of gradients, w.r.t parameters of interest, is a key aspect of deterministic algorithms like Gauss-Newton, Levenberg–Marquardt, Occam’s inversion [1] as well as statistical algorithms like Hamiltonian Monte Carlo (HMC) [2]. Common nonintrusive methods like finite differences compute the gradient by taking differences between the response at the current model and at a perturbed model, such methods suffer from certain drawbacks. Two types of errors stand out in particular: numerical error involved with truncation of Taylor’s series and round-off error involved with finite precision arithmetic of computers [3]. For more robust analysis, this chapter focuses on methods where the gradient is computed directly from the analytical/numerical model by enlisting the sensitivity equations.

The type of solutions that can be obtained from inverse problems is guided by the regularization term. When the accurate detection of the geometry of embedded objects (or detection of discontinuities) is of interest, $L_1$-norm-based difference priors have proven to be useful [4, 5]. However, special techniques need to be used to accommodate the nondifferentiability of the $L_1$-norm [1]. Hence, to make use of a large library of adjoint-based inversion solvers (that use gradients w.r.t the parameters), $L_2$-norm-based priors which readily allow for differentiation are more popular in practice. However, $L_2$-norm regularization or Gaussian priors in a stochastic sense only admit smooth solutions. Hence, in order to still be able to capture discontinuities, this paper explicitly parameterizes the shape of the boundary (or the geometry of the domain). This approach thereby considers two sets of parameters, one related to the spatial field and the second with the geometry parameters. This approach is of course only applicable when the unknown geometry can be parameterized explicitly.

Gradients have to be computed w.r.t both spatial as well geometry parameters. While sensitivity analysis of spatial parameters is usually straightforward, computation of the gradients w.r.t geometry parameters [6, 7] needs to be done more carefully and consider aspects like mesh distortion. This chapter develops on the method presented in [8] and briefly details the simultaneous spatial field and geometry update within the HMC statistical framework. Similar to [8], the Karhunen-Loève (K-L) expansion is used for discretization of the random field, with the difference that the complete theoretical basis is considered. The complete integral eigenvalue problem or IEVP is solved and the associated gradients are computed. The interaction between the geometry and spatial parameters due to the domain of definition of the IEVP is also detailed. It should be noted that the procedure detailed above is applicable to both the direct differentiation and adjoint methods [9] of sensitivity analysis.

The entire numerical study is done with a focus on aspects related to gradients and not with the aim to solve the inverse problem. Nevertheless, a forward model is still required for the computation of gradients. The chapter begins with a description of the forward model, observation equation and the discretization of the spatial random field in Section 2. This is followed by a brief description in Section 3 of the HMC-based methods developed in [8, 10] for simultaneous spatial field and geometry detection. Section 3.3 introduces the new gradients obtained when the complete IEVP is considered. It is shown how the gradients of the eigenvectors involve the computation of a Moore-Penrose pseudoinverse. Finally, the gradient computation procedure is validated and a convergence study is done in Section 4.

2. Inversion preliminaries

2.1 Governing and observation equations

Consider a linear steady seepage flow problem defined on a domain $\mathbf{z}\mathbf{\in }\mathrm{\Omega }\subset {\mathbb{R}}^d$, $d\in \left\{2,3\right\}$, where $\mathbf{k}\left(\mathbf{z}\right)$ is a symmetric spatially varying hydraulic conductivity matrix, $h\left(\mathbf{z}\right)$ is the hydraulic head, and $Q\left(\mathbf{z}\right)$ is a source term as shown in Eq. (1) below

$$\nabla .\mathit{-}\mathbf{k}\left(\mathbf{z}\right)\nabla h\left(\mathbf{z}\right)\mathbf{=}Q\left(\mathbf{z}\right).$$E1

Standard Dirichlet: $h=\overline{h}$ on ${\mathrm{\Gamma }}_D$, and Neumann: $f_n=\mathit{f}\mathit{.}\mathit{n}\mathit{=}\mathit{-}\mathbf{k}\nabla h.\mathit{n}$ on ${\mathrm{\Gamma }}_N$ boundary conditions are applied. Following a spatial discretization of the weak form of the PDE using the finite element method, the governing equation can be written as:

$$\mathbf{K}\left(\mathbf{\theta }\right)\mathbf{h}=\mathbf{q},$$E2

where $\mathbf{K}$ is the global hydraulic conductivity matrix and $\mathbf{h}$ and $\mathbf{q}$ are the nodal hydraulic head and flux vectors respectively. The parameter vector $\mathbf{\theta }\mathbf{=}\left({}^1\mathbf{\theta },{}^2\mathbf{\theta }\right)\in {\mathbb{R}}^K$ in Eq. (2), includes all the unknowns related to the inverse problem. In this study, the unknowns are divided into two sets [8]: ${}^1\mathbf{\theta }\in {\mathbb{R}}^{K_1}$, related to the spatial discretization of a random field and ${}^2\mathbf{\theta }\in {\mathbb{R}}^{K_2}$, and related to the definition of the geometry of the domain $\mathrm{\Omega }$. Also, consider an observation model relating the state vector $\mathbf{m}=\left(\mathbf{h},\mathbf{q}\right)$ to the discrete observations $\mathbf{y}$, through a map $\mathbf{H}$ that is independent of $\mathbf{\theta }$ i.e.

$$\mathbf{y}=\mathbf{Hm}\left(\mathbf{\theta }\right)+\mathbf{r}.$$E3

The error in Eq. (3) is modeled as Gaussian $\mathbf{r}\sim \mathbb{N}\left(\mathbf{0},\mathbf{R}\right)$ with a known covariance matrix $\mathbf{R}$. Parameter estimation is done in a probabilistic sense using Bayesian inference. Starting with a Gaussian prior distribution $p\left(\mathbf{\theta }\right)=\mathbb{N}\left(\mathbf{\theta }|\mathbf{0},{\mathbf{\Sigma }}_{\mathbf{\theta }}\right)$ and a likelihood distribution $p\left(\mathbf{y}|\mathbf{\theta }\right)=\mathbb{N}\left(\mathbf{y}|\mathbf{Hm}\left(\mathbf{\theta }\right),\mathbf{R}\right)$, the posterior distribution is written as:

$$p\left(\mathbf{\theta }|\mathbf{y}\right)\propto p\left(\mathbf{y}|\mathbf{\theta }\right)p\left(\mathbf{\theta }\right).$$E4

Except for linear Gaussian observation models, the posterior cannot be computed analytically and is usually evaluated using MCMC sampling algorithms.

2.2 Karhunen-Loève (K-L) expansion

The parameter vector ${}^1\mathbf{\theta }$ defined in Section 2.1 is associated with a continuous hydraulic conductivity spatial random field $k\left(\mathbf{z},\omega \right)$, where $\mathbf{z}$ is defined on the domain $\mathrm{\Omega }$ and $\omega$ belongs to the space of random events $\mathrm{\Theta }$. Let the expected value of the random field be denoted as $\overline{k}\mathit{:}\mathrm{\Omega }\to \mathbb{R}$ and the autocovariance function $C:\mathrm{\Omega }\times \mathrm{\Omega }\to \mathbb{R}$ be defined as $C\left(\mathbf{z},{\mathbf{z}}^{\prime }\right)\mathit{=}\sigma \left(\mathbf{z}\right)\sigma \left({\mathbf{z}}^{\prime }\right)\rho \left(\mathbf{z},{\mathbf{z}}^{\prime }\right)$. Here $\sigma \mathit{:}\mathrm{\Omega }\to \mathbb{R}$ is the standard deviation function and $\rho :\mathrm{\Omega }\times \mathrm{\Omega }\to \left[-1,1\right]$ is the autocorrelation coefficient function. The study in this chapter is confined to Gaussian random fields that can be defined completely by their mean and autocovariance functions.

The Karhunen-Loève (K-L) expansion method is a series expansion method for the discretization of random fields which is based on the spectral decomposition of the autocovariance function. It can be shown that a random field can be written as an infinite sum [11]:

$$k\left(\mathbf{z},\omega \right)=\overline{k}\left(\mathbf{z}\right)+\sum _{k=1}^{\infty }\sqrt{{\lambda }_k}{{}^1\theta }_k{\varphi }_k\left(\mathbf{z}\right),$$E5

where ${{}^1\theta }_k\mathit{:}\mathrm{\Theta }\to \mathbb{R}$ are standard uncorrelated random variables, ${\lambda }_k$ are the eigenvalues (always non-negative) and ${\varphi }_k\left(\mathbf{z}\right)$ are the eigenfunctions of the linear operator related to the covariance kernel $C$. They can be obtained by solving the homogeneous Fredholm integral eigenvalue problem (IEVP) on the domain $\mathrm{\Omega }$:

$${\int }_{\mathrm{\Omega }}C\left(\mathbf{z},{\mathbf{z}}^{\prime }\right){\varphi }_k\left({\mathbf{z}}^{\prime }\right)\mathrm{d}{\mathbf{z}}^{\prime }={\lambda }_k{\varphi }_k\left(\mathbf{z}\right).$$E6

The autocovariance function is symmetric, bounded, and positive semi-definite and has the spectral decomposition $C\left(\mathbf{z},{\mathbf{z}}^{\prime }\right)\mathit{=}{\sum }_{k=1}^{\infty }{\lambda }_k{\varphi }_k\left(\mathbf{z}\right){\varphi }_k\left({\mathbf{z}}^{\prime }\right)$. The eigenfunctions are orthogonal, and in a normalized form satisfy the condition ${\int }_{\mathrm{\Omega }}{\varphi }_k\left(\mathbf{z}\right){\varphi }_{\mathit{l}}\left(\mathbf{z}\right)\mathrm{d}\mathbf{z}={\delta }_{\mathit{kl}}$, where ${\delta }_{\mathit{kl}}$ is the Kronecker delta. In the case of Gaussian random fields, the random variables ${{}^1\theta }_k$ are also independent and follow the standard normal distribution.

In practice, the eigenvalues decay exponentially fast for smooth functions and algebraically fast for non-smooth autocovariance kernels and the K-L expansion is usually truncated after $K_1$ terms. If the eigenvalues are arranged in descending order such that ${\lambda }_1>{\lambda }_2>\dots >{\lambda }_{K_1}$, then accompanied by the associated eigenfunctions, the truncated K-L expansion approximation of the random field can be written as

$$\widehat{k}\left(\mathbf{z},\omega \right)=\overline{k}\left(\mathbf{z}\right)+\sum _{k=1}^{K_1}\sqrt{{\lambda }_k}{{}^1\theta }_k{\varphi }_k\left(\mathbf{z}\right),$$E7

The truncated K-L expansion approximation is optimal in the sense that, for a fixed number of terms $K_1$, the mean square error over the domain is minimized [12]. A global error measure related to random field discretization is called the mean error variance ${\overline{\epsilon }}_{\sigma }$ and is defined as [13]:

$${\overline{\epsilon }}_{\sigma }\left(\mathbf{z}\right)=\frac{1}{\left|\mathrm{\Omega }\right|}{\int }_{\mathrm{\Omega }}\frac{\mathrm{Var}\left[k\left(\mathbf{z},\omega \right)-\widehat{k}\left(\mathbf{z},\omega \right)\right]}{\mathrm{Var}\left[k\left(\mathbf{z},\omega \right)\right]}\mathrm{d}\mathbf{z}.$$E8

It can be shown that the variance of the truncated K-L expansion $\widehat{k}\left(\mathbf{z},\omega \right)$ is:

$$\mathrm{Var}\left[\widehat{k}\left(\mathbf{z},\omega \right)\right]=\sum _{k=1}^{K_1}{\lambda }_k{\varphi }_k^2\left(\mathbf{z}\right).$$E9

Using the property $\mathbb{E}\left[{{}^1\theta }_k{{}^1\theta }_l\right]={\delta }_{\mathit{kl}}$, the mean error variance can be calculated as [14]:

$${\overline{\epsilon }}_{\sigma ,\mathrm{KL}}=1-\frac{1}{\left|\mathrm{\Omega }\right|{\sigma }^2}\sum _{k=1}^{K_1}{\lambda }_k.$$E10

The derivation for Eq. (10) assumes the random field to be homogeneous, i.e., $\sigma \left(\mathbf{z}\right)\mathbf{=}\sigma$. We only consider the case where the prior random field is homogeneous and Gaussian. This assumption is for numerical convenience and does not limit the posterior, which can be non-Gaussian and non-homogeneous [15].

2.3 Galerkin finite element method to solve the integral eigenvalue problem

The Galerkin Finite Element Method (FEM) is used to solve the IEVP on $\mathrm{\Omega }$. The eigenfunctions are approximated with the help of the shape functions $N_j\mathbf{:}\mathrm{\Omega }\to \mathbb{R}$ of the FE mesh, and is represented as:

$${\varphi }_k\left(\mathbf{z}\right)\mathbf{\approx }\sum _{j=1}^n{d}_{\mathit{kj}}{N}_j\left(\mathbf{z}\right),$$E11

where the coefficients $d_{\mathit{kj}}\in \mathbb{R}$ are unknown and $n$ is the number of nodes in the FE mesh. Substitution of Eq. (11) into Eq. (6), yields the residual:

$$r\left(\mathbf{z}\right)\mathbf{=}\sum _{j=1}^n{d}_{\mathit{kj}}\left({\int }_{\mathrm{\Omega }}C\left(\mathbf{z},{\mathbf{z}}^{\prime }\right)N_j\left({\mathbf{z}}^{\prime }\right)\mathrm{d}{\mathbf{z}}^{\prime }-{\lambda }_j{N}_j\left(\mathbf{z}\right)\right),$$E12

In the Galerkin method, the unknown coefficients are determined by making the residual $r\left(\mathbf{z}\right)$ orthogonal to the space spanned by the shape functions i.e.

$${\int }_{\mathrm{\Omega }}r\left(\mathbf{z}\right)N_i\left({\mathbf{z}}^{\prime }\right)\mathrm{d}{\mathbf{z}}^{\prime }=0\forall j=1,\dots ,n.$$E13

This results in a generalized eigenvalue problem

$$\mathbf{B}{\mathbf{d}}_k\mathit{=}{\lambda }_k\mathbf{M}{\mathbf{d}}_k,$$E14

where

$$B_{\mathit{ij}}={\int }_{\mathrm{\Omega }}{N}_i\left(\mathbf{z}\right){\int }_{\mathrm{\Omega }}C\left(\mathbf{z},{\mathbf{z}}^{\prime }\right)N_j\left({\mathbf{z}}^{\prime }\right)\mathrm{d}{\mathbf{z}}^{\prime }\mathrm{d}\mathbf{z}\text{and}$$

$$M_{\mathit{ij}}={\int }_{\mathrm{\Omega }}{N}_i\left(\mathbf{z}\right)N_j\left(\mathbf{z}\right)\mathrm{d}\mathbf{z}.$$E15

Both $\mathbf{B}$ and $\mathbf{M}$ are $n\times n$ matrices that involve integrals over the domain $\mathrm{\Omega }$. Hence the actual geometry of the domain has to be considered for integration. The maximum number of available eigenpairs is $n$, but in practice, the K-L expansion can usually be truncated at $K_1$ terms such that $K_1\ll n$. As such, for computational efficiency, it is sufficient to compute the first $K_1$ eigenpairs only, which can be done through the Lanczos algorithm.

3. Simultaneous geometry and spatial field detection

3.1 Hamiltonian Monte Carlo

Consider a parameter space $\mathbf{\theta }\in {\mathbb{R}}^K$ augmented with equidimensional momentum variables $\mathbf{p}\in {\mathbb{R}}^K$ and a joint probability distribution with density $p\left(\mathbf{\theta },\mathbf{p}\right)$ defined over this augmented space. If the underlying distribution over the momentum variables is chosen to be a Gaussian: $p\left(\mathbf{p}\right)\mathit{\equiv }\mathbb{N}\left(\mathbf{p}|\mathbf{0},\mathcal{M}\right)$, where $\mathcal{M}$ is user-specified and independent of $\mathbf{\theta }$, then a joint probability distribution can be defined as $p\left(\mathbf{\theta },\mathbf{p}\right)=p\left(\mathbf{p}\right)p\left(\mathbf{\theta }|\mathbf{y}\right)$. The Hamiltonian $H:{\mathbb{R}}^K\times {\mathbb{R}}^K\to \mathbb{R}$ is then defined as:

$$H\left(\mathbf{\theta },\mathbf{p}\right)=-logp\left(\mathbf{p}\right)-logp\left(\mathbf{\theta }|\mathbf{y}\right)$$E16

The second term on the right of Eq. (17) is generally called $\phi :{\mathbb{R}}^K\to \mathbb{R}$ and is given as $\phi \left(\mathbf{\theta }\right)=-logp\left(\mathbf{\theta }|\mathbf{y}\right)$.

The introduction of the momentum variables allows for the generation of trajectories through conservative Hamiltonian dynamics [16], which are given as:

$$\left(\begin{array}{c}\frac{d\mathbf{\theta }}{\mathit{dt}}\\ \frac{d\mathbf{p}}{\mathit{dt}}\end{array}\right)=\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]\left(\begin{array}{c}\frac{\partial H}{\partial \mathbf{\theta }}\\ \frac{\partial H}{\partial \mathbf{p}}\end{array}\right)=\left(\begin{array}{c}-\frac{\partial \phi \left(\mathbf{\theta }\right)}{\partial \mathbf{\theta }}\\ {\mathcal{M}}^{-1}\mathbf{p}\end{array}\right).$$E17

These dynamics are exactly reversible (provided the gradient $\frac{\partial \phi \left(\mathbf{\theta }\right)}{\partial \mathbf{\theta }}$ is one-to-one) and preserve volume as Eq. (18) is just a rotation transformation in $\mathbf{\theta }\mathbf{-}\mathbf{p}$ space. Except for simple problems Eq. (18) cannot be solved analytically and is usually solved using the leapfrog method, which is a second-order accurate numerical integrator given as:

$$\mathbf{p}\left(t+\frac{ϵ}{2}\right)\mathbf{=}\mathbf{p}\left(t\right)\mathbf{-}\frac{ϵ}{2}\frac{\partial \phi \left(\mathbf{\theta }\left(t\right)\right)}{\partial \mathbf{\theta }},$$E18

$$\mathbf{\theta }\left(t+ϵ\right)\mathit{=}\mathbf{\theta }\left(t\right)\mathit{+}ϵ{\mathcal{M}}^{\mathit{-}1}\mathbf{p}\left(t+\frac{ϵ}{2}\right)\text{and}$$E19

$$\mathbf{p}\left(t+ϵ\right)\mathbf{=}\mathbf{p}\left(t+\frac{ϵ}{2}\right)\mathbf{-}\frac{ϵ}{2}\frac{\partial \phi \left(\mathbf{\theta }\left(t+\epsilon \right)\right)}{\partial \mathbf{\theta }}.$$E20

Starting from a point $\left({\mathbf{\theta }}^j,{\mathbf{p}}^j\right)$, these equations are applied repeatedly for $L$ steps, each with a step-size $ϵ$, to determine a transition to a new point $\left({\mathbf{\theta }}^{j+1},{\mathbf{p}}^{j+1}\right)$, which lies on the same Hamiltonian level-set as $\left({\mathbf{\theta }}^j,{\mathbf{p}}^j\right)$. The deterministic part of Hamiltonian Monte Carlo (HMC) [2] is defined by Eqs. (19)–(21). The stochastic part of HMC comes from resampling $\mathbf{p}\sim \mathbb{N}\left(\mathbf{0},\mathcal{M}\right)$. The statistical efficiency of Hamiltonian Monte Carlo stems from the fact that the gradient-guided transitions can propose new points that are “far-away” from the starting point, thereby enabling efficient sampling of the posterior. This is in contrast to the random nature of transitions, which suffer from the curse of dimensionality [17], in conventional MCMC algorithms. As shown in Eqs. (19)–(21), critical to the success of HMC, is the computation of the gradient $\frac{\partial \phi \left(\mathbf{\theta }\right)}{\partial \mathbf{\theta }}$. Special attention must be paid to maintaining the reversibility of the transitions to satisfy the detailed balance condition [18] for MCMC algorithms. This is detailed along with the gradient computation procedure in the following sections.

3.2 Parameter update using the mesh moving method

The leapfrog equations determine an update in $\mathbf{\theta }\mathbf{-}\mathbf{p}$ space. In particular, the update from $\mathbf{\theta }\left(t\right)\mathit{\to }\mathbf{\theta }\left(t+ϵ\right)$ in Eq. (20), defines not only a new realization of the random field, but also a new domain, i.e., $\mathrm{\Omega }\left({}^2\mathbf{\theta }\right)$. Without loss of generality, consider a domain in 2D discretized such that $\mathbf{Z}\left({}^2\mathbf{\theta }\right)\in {\mathbb{R}}^2$ is the nodal coordinate vector of all the nodes of the finite element mesh. Let ${\mathbf{Z}}_v\left({}^2\mathbf{\theta }\right)\in {\mathbb{R}}^2$ represent a subset of this vector that includes only the node coordinates at the piping zone boundary. Koch et al. [10] show that the computation of the gradient $\frac{\partial \phi \left(\mathbf{\theta }\right)}{\partial \mathbf{\theta }}$, by analytical methods [6], ultimately involves the computation of the gradient of the nodal coordinate vector $\frac{\partial \mathbf{Z}\left(\mathbf{\theta }\right)}{\partial \mathbf{\theta }}$.

Computation of the nodal coordinate vector gradient requires the definition of a differentiable map which is additionally reversible and one-to-one, to satisfy the detailed balance condition of MCMC. One such map proposed in [10], is to update from an arbitrarily fixed reference domain ${\mathrm{\Omega }}^{\mathit{ref}}\left({{}^2\mathbf{\theta }}^{\mathit{ref}}\right)$, defined by an arbitrary parameter ${{}^2\mathbf{\theta }}^{\mathit{ref}}$. Let ${\mathbf{Z}}^{\mathit{ref}}\left({{}^2\mathbf{\theta }}^{\mathit{ref}}\right)\in {\mathbb{R}}^{n\times 2}$ be the nodal coordinate vector on the discretized domain ${\mathrm{\Omega }}^{\mathit{ref}}$ and ${\mathbf{Z}}_v^{\mathit{ref}}\left({{}^2\mathbf{\theta }}^{\mathit{ref}}\right)\in {\mathbb{R}}^{n_v\times 2}$ represent a subset of this vector that includes only the coordinates of the $n_v$ nodes at the piping zone boundary. The nodal coordinates ${\mathbf{Z}}_v\left({}^2\mathbf{\theta }\right)$ and ${\mathbf{Z}}_v^{\mathit{ref}}\left({{}^2\mathbf{\theta }}^{\mathit{ref}}\right)$ can be determined explicitly if ${}^2\mathbf{\theta }$ and ${{}^2\mathbf{\theta }}^{\mathit{ref}}$ are known respectively. Following the update, ${}^2\mathbf{\theta }\left(t\right)\mathit{\to }{}^2\mathbf{\theta }\left(t+\epsilon \right)$, ${\mathbf{Z}}_v\left({}^2\mathbf{\theta }\left(t+\epsilon \right)\right)$ is available, and an artificial elastic deformation problem can be set up from the arbitrary known reference domain ${\mathrm{\Omega }}^{\mathit{ref}}\left({{}^2\mathbf{\theta }}^{\mathit{ref}}\right)$ to the current domain $\mathrm{\Omega }\left({}^2\mathbf{\theta }\left(t+\epsilon \right)\right)$. The prescribed displacements ${\mathbf{u}}_v^{\mathit{ref}}\in {\mathbb{R}}^{n_v\times 2}$ for the elastic deformation problem are given as:

$${\mathbf{u}}_v^{\mathit{ref}}={\mathbf{Z}}_v\left({}^2\mathrm{\theta }\left(t+\epsilon \right)\right)-{\mathbf{Z}}_v^{\mathit{ref}}\left({{}^2\mathrm{\theta }}^{\mathit{ref}}\right).$$E21

The entire mesh is moved and the new nodal coordinates can be determined as:

$$\mathbf{Z}\left({}^2\mathbf{\theta }\left(t+\epsilon \right)\right)\mathit{=}{\mathbf{Z}}^{\mathit{ref}}\left({{}^2\mathbf{\theta }}^{\mathit{ref}}\right)\mathit{+}{\mathbf{u}}^{\mathit{ref}},$$E22

where ${\mathbf{u}}^{\mathit{ref}}\in {\mathbb{R}}^{n\times 2}$ represents the displacement of all the nodes from the reference domain to the current domain.

The displacements in the elastic deformation step can cause distortions in the mesh, especially in regions where large deformation is expected, i.e., near the piping zone boundary. To maintain a good mesh quality for computation purposes, a mesh moving method [19] is used. The simple idea is to scale the elastic modulus $E_e^{\mathit{ref}}$ of each element (in the reference domain) with the determinant of the Jacobian $\left|{\mathbf{J}}_e^{\mathit{ref}}\right|$ in the elastic deformation step:

$${\stackrel{\sim }{E}}_e^{\mathit{ref}}\mathbf{=}{E}_e^{\mathit{ref}}{\left(1/\left|{\mathbf{J}}_e^{\mathit{ref}}\right|\right)}^{{\chi }^{\mathit{ref}}},$$E23

where ${\chi }^{\mathit{ref}}$ is an arbitrary positive scaling parameter. The Poisson’s ratio of the reference domain ${\nu }^{\mathit{ref}}$ is also chosen arbitrarily. The performance of the algorithm has been shown [20] to be invariant to the choice of these reference parameters. The net effect of such scaling is that small elements become rigid and larger elements become more flexible. Hence, if the reference domain mesh is constructed carefully such that small elements are placed in regions where large distortion is expected, i.e., near the piping zone, and larger elements are placed in regions of less expected distortion, the method is expected to yield a good mesh quality in the elastic deformation stage. The map in Eq. (23) is differentiable and helps determine $\frac{\partial \mathbf{Z}\left(\mathbf{\theta }\right)}{\partial \mathbf{\theta }}$.

3.3 Gradient computation

Analytical methods for the computation of the gradient are intrusive and involve gradients of the steady seepage flow forward solver. From the definition of $\phi \left(\mathbf{\theta }\right)=-logp\left(\mathbf{\theta }|\mathbf{y}\right)$, it is apparent that the computation of $\frac{\partial \phi \left(\mathbf{\theta }\right)}{\partial \mathbf{\theta }}$ requires the computation of $\frac{\partial \mathbf{m}}{\partial \mathbf{\theta }}=\left(\frac{\partial \mathbf{h}}{\partial \mathbf{\theta }},\frac{\partial \mathbf{q}}{\partial \mathbf{\theta }}\right)$. These terms can be computed by taking the derivative of Eq. (2) and is given as:

$$\frac{\partial \mathbf{K}}{\partial \mathbf{\theta }}\mathbf{h}+\mathbf{K}\frac{\partial \mathbf{h}}{\partial \mathbf{\theta }}=\frac{\partial \mathbf{q}}{\partial \mathbf{\theta }}.$$E24

Given standard boundary conditions, this equation can be solved to obtain $\frac{\partial \mathbf{h}}{\partial \mathbf{\theta }}$ and $\frac{\partial \mathbf{q}}{\partial \mathbf{\theta }}$, provided $\frac{\partial \mathbf{K}}{\partial \mathbf{\theta }}$ is known.

Following the standard procedure in FEM, the hydraulic conductivity matrix at the element level ${\mathrm{\Omega }}_e$ in the current domain is given as:

$${\mathbf{K}}_e=\underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{G}}^{\mathrm{T}}\widehat{k}\left(\mathbf{\theta },\mathbf{z}\right)\mathbf{G}\left|{\mathbf{J}}_e\right|\mathrm{d}\mathbf{\xi },$$E25

where $\widehat{k}$ represents the hydraulic conductivity spatial field obtained from the truncated K-L expansion, $\mathbf{G}$ contains the derivatives (w.r.t $\mathbf{z}$) of the shape functions $N_j$ described earlier in Section 2.3, $\left|{\mathbf{J}}_e\right|$ is the determinant of the Jacobian matrix associated with the isoparametric transformation $\left({\xi }_1,{\xi }_2\right)\to \left({\mathrm{z}}_1,{\mathrm{z}}_2\right)$ and ${\mathrm{\Omega }}_e$ is the region occupied by the parent element related to the isoparametric transformation. The gradient of ${\mathbf{K}}_e$ with respect to the spatial parameters ${}^1\mathbf{\theta }$ can be computed as:

$$\frac{\partial {\mathbf{K}}_e}{\partial {}^1\mathbf{\theta }}\mathit{=}\underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{G}}^{\mathrm{T}}\frac{\partial \widehat{k}}{\partial {}^1\mathbf{\theta }}\mathbf{G}\left|{\mathbf{J}}_e\right|\mathrm{d}\mathbf{\xi },$$E26

where the gradient of the hydraulic conductivity field w.r.t ${}^1\mathbf{\theta }$ is easily obtained by differentiating Eq. (7) and is given as:

$$\frac{\partial \widehat{k}}{\partial {}^1{\theta }_j}\mathbf{=}\sqrt{{\lambda }_j}{\varphi }_j,$$E27

The gradient of ${\mathbf{K}}_e$ w.r.t the geometry parameters ${}^2\mathbf{\theta }$ can be written as [6]:

$$\frac{\partial {\mathbf{K}}_e}{\partial {}^2\mathbf{\theta }}={\int }_{{\mathrm{\Omega }}_e}\left(\frac{\partial {\mathbf{G}}^{\mathrm{T}}}{\partial {}^2\mathbf{\theta }}\widehat{k}\mathbf{G}\left|{\mathbf{J}}_{\mathrm{e}}\right|+{\mathbf{G}}^{\mathrm{T}}\widehat{k}\frac{\partial \mathbf{G}}{\partial {}^2\mathbf{\theta }}\left|{\mathbf{J}}_{\mathrm{e}}\right|+{\mathbf{G}}^{\mathrm{T}}\widehat{k}\mathbf{G}\frac{\partial \left|{\mathbf{J}}_{\mathrm{e}}\right|}{\partial {}^2\mathbf{\theta }}\mathbf{+}{\mathbf{G}}^{\mathrm{T}}\frac{\partial \widehat{k}}{\partial {}^2\mathbf{\theta }}\mathbf{G}\left|{\mathbf{J}}_{\mathrm{e}}\right|\right)\mathrm{d}\mathbf{\xi }.$$E28

Formulas for the calculation of the gradients $\frac{\partial \mathbf{G}}{\partial {}^2\mathbf{\theta }}$ and $\frac{\partial \left|{\mathbf{J}}_{\mathrm{e}}\right|}{\partial {}^2\mathbf{\theta }}$ can readily be found in literature [6]. It is clear from Eq. (6) that the computation of the eigenvalues and eigenfunctions depends on the definition of the domain $\mathrm{\Omega }\left({}^2\mathbf{\theta }\right)$. Hence, the gradient $\frac{\partial \widehat{k}\left(\mathbf{\theta },\mathbf{z}\right)}{\partial {}^2\mathbf{\theta }}$ will involve the gradients of the eigenvalues and eigenvectors and is written as:

$$\frac{\partial \widehat{k}}{\partial {}^2\mathbf{\theta }}\mathbf{=}\sum _{j=1}^{K_1}\left(\frac{1}{2\sqrt{{\lambda }_j}}\frac{\partial {\lambda }_j}{\partial {}^2\mathbf{\theta }}{\varphi }_j\mathbf{+}\sqrt{{\lambda }_q}\frac{\partial {\varphi }_j}{\partial {}^2\mathbf{\theta }}\right){{}^1\theta }_j.$$E29

The gradient of the eigenvalues and eigenvectors of the generalized eigenvalue problem in Eq. (14) w.r.t. ${}^2\mathbf{\theta }$ are written as [21]:

$$\frac{\partial {\lambda }_j}{\partial {}^2\mathbf{\theta }}\mathit{=}{\mathbf{d}}_j^{\mathrm{T}}\left(\frac{\partial \mathbf{B}}{\partial {}^2\mathbf{\theta }}-{\lambda }_j\frac{\partial \mathbf{M}}{\partial {}^2\mathbf{\theta }}\right){\mathbf{d}}_j,$$E30

$$\frac{\partial {\mathbf{d}}_j}{\partial {}^2\mathbf{\theta }}\mathit{=}{\left({\lambda }_j\mathbf{M}\mathit{-}\mathbf{B}\right)}^{†}\left(\frac{\partial \mathbf{B}}{\partial {}^2\mathbf{\theta }}-{\lambda }_j\frac{\partial \mathbf{M}}{\partial {}^2\mathbf{\theta }}\right){\mathbf{d}}_j\mathit{-}\frac{1}{2}\left({\mathbf{d}}_j^{\mathrm{T}}\frac{\partial \mathbf{M}}{\partial {}^2\mathbf{\theta }}{\mathbf{d}}_j\right){\mathbf{d}}_j.$$E31

As ${\mathbf{d}}_j$ lies in the null space of ${\lambda }_j\mathbf{M}\mathit{-}\mathbf{B}$, the inverse ${\left({\lambda }_j\mathbf{M}\mathit{-}\mathbf{B}\right)}^{-1}$ cannot be computed and a generalized inverse called the Moore-Penrose pseudoinverse ${\left(\bullet \right)}^{†}$ is employed. In this study, the Moore-Penrose pseudoinverse is calculated by first carrying out an SVD of the matrix ${\lambda }_j\mathbf{M}\mathit{-}\mathbf{B}$ and then eliminating the smallest singular values below a tolerance level. This is then followed by taking a standard inverse of the SVD.

Considering the same mesh that is used for the solution of the steady seepage flow problem to be used for the discretization of the random field, the gradients of the matrices $\mathbf{B}$ and $\mathbf{M}$ at an elemental level in isoparametric space are given as:

$$\frac{\partial B_{\mathit{ij}}}{\partial {}^2\mathbf{\theta }}\mathit{=}\begin{array}{l}{\int }_{{\mathrm{\Omega }}_e}{N}_i{\int }_{{\mathrm{\Omega }}_e}C\left(\mathbf{z},{\mathbf{z}}^{\prime }\right)N_j\left(\frac{\partial \left|{\mathbf{J}}_{\mathrm{e}}\left({\mathbf{z}}^{\prime }\right)\right|}{\partial {}^2\mathbf{\theta }},\left|{\mathbf{J}}_{\mathrm{e}}\left({\mathbf{z}}^{\prime }\right)\right|+\left|{\mathbf{J}}_{\mathrm{e}}\left(\mathbf{z}\right)\right|\frac{\partial \left|{\mathbf{J}}_{\mathrm{e}}\left({\mathbf{z}}^{\prime }\right)\right|}{\partial {}^2\mathbf{\theta }}\right)\mathrm{d}{\mathbf{\xi }}^{\prime }\mathrm{d}\mathbf{\xi }\\ +{\int }_{{\mathrm{\Omega }}_e}{N}_i{\int }_{{\mathrm{\Omega }}_e}\left(\frac{\partial C\left(\mathbf{z},{\mathbf{z}}^{\prime }\right)}{\partial {}^2\mathbf{\theta }}\right)N_j\left|{\mathbf{J}}_{\mathrm{e}}\left(\mathbf{z}\right)\right|\left|{\mathbf{J}}_{\mathrm{e}}\left({\mathbf{z}}^{\prime }\right)\right|\mathrm{d}{\mathbf{\xi }}^{\prime }\mathrm{d}\mathbf{\xi },\end{array}$$E32

$$\frac{\partial M_{\mathit{ij}}}{\partial {}^2\mathbf{\theta }}\mathit{=}{\int }_{{\mathrm{\Omega }}_e}{N}_i{N}_j\frac{\partial \left|{\mathbf{J}}_{\mathrm{e}}\left(\mathbf{z}\right))\right|}{\partial {}^2\mathbf{\theta }}\mathrm{d}\mathbf{\xi }.$$E33

In this study, the squared exponential autocorrelation coefficient function has been used, i.e.,

$$\rho \left(\mathbf{z},{\mathbf{z}}^{\prime }\right)=exp\left(\frac{-{\left|\mathbf{z}-{\mathbf{z}}^{\prime }\right|}^{\mathbf{2}}}{l_c^2}\right),$$E34

where $l_c$ is the correlation length of the random field. Assuming $\sigma \left(\mathbf{z}\right)\mathbf{=}\sigma \left({\mathbf{z}}^{\prime }\right)=\sigma$, the gradient of the autocovariance function $C\left(\mathbf{z},{\mathbf{z}}^{\prime }\right)$ w.r.t ${}^2\mathbf{\theta }$ is given as:

$$\frac{\partial C\left(\mathbf{z},{\mathbf{z}}^{\prime }\right)}{\partial {}^2\mathbf{\theta }}=-\frac{2}{l_c^2}\left({\left(\frac{\partial \mathbf{z}}{\partial {}^2\mathbf{\theta }}-\frac{\partial \mathbf{z}}{\partial {}^2\mathbf{\theta }}\right)}^{\mathrm{T}}\left(\mathbf{z}\mathbf{-}{\mathbf{z}}^{\prime }\right)\right)C\left(\mathbf{z},{\mathbf{z}}^{\prime }\right).$$E35

This completes the definition of all the terms required to compute the HMC gradient

$$\frac{\partial \phi \left(\mathbf{\theta }\right)}{\partial \mathbf{\theta }}.$$

4. Numerical implementation and results

4.1 Observation data

A seepage zone containing a piping region of length $l$ and width $w$, shown in Figure 1, is chosen to generate synthetic observations for inversion. A steady seepage flow problem is solved on the domain defined by $l=0.15$ m and $w=0.05$ m. The left and right boundaries marked in blue are Dirichlet boundaries and the top and bottom boundaries are no-flow Neumann boundaries. Linear shape functions are used in the FE solution of the forward problem. Ten sets of observations (hydraulic head and total outward normal flux) are made by increasing the hydraulic head at the left boundary as shown by point A in Figure 2. The hydraulic head at the right boundary is fixed at 0. The corresponding hydraulic head recorded at observation points B, C, D, and E are shown in Figure 2. The total outward normal flux from the right boundary is summed and also shown on the right in Figure 2. Observation data is generated assuming a constant hydraulic conductivity field, i.e., $k^{\text{true}}\left(\mathbf{z}\right)=0.001$ m/s. The standard deviation of the observation noise for the hydraulic head data and outward normal flux data is taken to be $1$ and 5% respectively.

Figure 1.

Discretized seepage domain containing piping zone, used to obtain observation data, i.e., hydraulic head $h$ data at points B, C, D and E, and total outward normal flux data $q$ from the boundary marked in blue on the right. The hydraulic head is constant on the left boundary. The number of nodes in the discretization is 3943.

Figure 2.

Observation data generated from numerical seepage flow experiment. Data at each time $t$ is obtained by solving a new seepage flow experiment with an input hydraulic head at the left boundary represented by A.

4.2 Inversion setup

Considering a log-normal hydraulic conductivity field $k\left(\mathbf{z}\right)$, inversion is carried out on the Gaussian field $\stackrel{\sim }{k}\left(\mathbf{z}\right)=log\left[k\left(\mathbf{z}\right)-\breve{k}\right]$, where $\breve{k}$ is a lower bound taken as ${10}^{-5}$ m/s. The log-normal construction enables a convenient choice of the standard deviation of the random field which is chosen as $\sigma \left(\mathbf{z}\right)\mathit{=}1$. The correlation length of the random field is assumed to be known and is taken to be equal to the length of the largest dimension of the domain i.e. $l_c=0.4\mathrm{m}$.

A preliminary assessment of the eigenvalues obtained by solving the generalized eigenvalue problem in Eq. (14) reveals that the K-L expansion can be truncated at $K_1=12$ terms. The eigenvalues are found to decay by approximately $3\times {10}^4$ times. This enables the determination of the number of terms in the spatial parameter vector ${}^1\mathbf{\theta }\mathit{=}\left({\theta }_1,\dots ,{\theta }_{12}\right)$. As mentioned in Section 3.2, once the geometry parameterization ${}^2\mathbf{\theta }\mathit{=}\left({}^2{\theta }_1,{}^2{\theta }_2\right)\mathit{=}\left({\theta }_{13},{\theta }_{14}\right)$ is known (see Figure 3(a)), an explicit function ${\mathbf{Z}}_v\left({}^2\mathbf{\theta }\right)$ can be constructed as:

Figure 3.

(a) Parameterizations of the piping zone boundary. Discretization of the reference domain defined by $\left({}^2{\theta }_1^{\mathit{ref}},{}^2{\theta }_2^{\mathit{ref}}\right)\mathit{=}\left(0.15\mathrm{m},0.05\mathrm{m}\right)$ with (b) 307 (c) 1476 and (d) 3018 nodes. Red dots indicate the position of observation points.

$${\mathbf{z}}_{v_j}^1\mathbf{=}\left(\begin{array}{c}{L}_1\mathbf{-}{}^2{\theta }_1\mathbf{+}\left(j-1\right)\frac{{}^2{\theta }_1}{n_1-1}\\ L_2-{}^2{\theta }_2\end{array}\right),$$E36

$${\mathbf{z}}_{v_j}^2\mathbf{=}\left(\begin{array}{c}{L}_1\mathbf{-}{}^2{\theta }_1\\ L_2\mathbf{-}{}^2{\theta }_2\mathbf{+}\left(j-1\right)\frac{{}^2{\theta }_2}{n_2-1}\end{array}\right).$$E37

Eqs. (37) and (38) can be readily differentiated to obtain $\frac{\partial {\mathbf{Z}}_v}{\partial \mathbf{\theta }}$, which can be used to compute other gradients $\frac{\partial \mathbf{Z}}{\partial \mathbf{\theta }}$, $\frac{\partial \phi }{\partial \mathbf{\theta }}$ etc. This completes the definition of the parameter vector.

As all updates are designed to take place from an arbitrary reference domain as mentioned in Eq. (23), the study of the convergence properties requires a focus on the reference mesh. The reference domain is arbitrarily chosen as the true domain used to generate the observations. Seven different reference configuration meshes with 128, 306, 497, 1038, 1476, 1860, and 3018 nodes are considered for the study of the convergence properties. Three of these meshes are shown in Figure 3. The artificial elastic properties selected for the mesh moving method are $E^{\mathit{ref}}\mathit{=}25\mathrm{MPa}$, ${\nu }^{\mathit{ref}}\mathit{=}0.25$ and ${\chi }^{\mathit{ref}}=1$. To reduce mesh distortion during the mesh moving stage, all the reference meshes have a unique construction, i.e., smaller elements (which behave rigidly) are placed close to the piping zone boundary and larger elements (which are more flexible) are placed further away.

4.3 Inversion and convergence analysis

There are no analytical solutions to verify the correctness of the gradient computation procedure mentioned in Section 3. Hence, for verification purposes, a short inversion analysis is done using HMC, where $i=150$ samples are drawn from the posterior. The number of leapfrog steps is variable and drawn from a Gaussian $L\sim \mathbb{N}\left(5,2\right)$. The prior for the parameters are chosen as ${}^1\mathbf{\theta }\sim \mathrm{N}\left(\mathbf{0},{\mathbf{I}}_{12}\right)$ and ${}^2\mathbf{\theta }\sim \mathrm{N}\left(\mathbf{0},0.1{\mathbf{I}}_2\right)$. The results presented in Figure 4 correspond to inversion done considering the reference mesh shown in Figure 3(b). The potential energy $\phi$ decreases rapidly in the first 20 steps and thereafter HMC begins the exploration of the region of the high posterior probability. Prior experience leads the authors to consider the performance of HMC to be appropriate and as such, is an indirect indicator that the gradient computation procedure is correct.

Figure 4.

Results from HMC showing (a) Decrease in potential energy as HMC progresses towards the high posterior probability region and samples of (b) ${\theta }_{13}$ and (c) ${\theta }_{14}$. The dashed lines represent the true values of the parameters used to generate the observations.

In order to study the convergence properties, suitable error measures are required. The rate of convergence of the truncated KL expansion of the prior random field, for a fixed number of parameters $K_1=12$, is determined through the relative mean error variance [14]:

$${\epsilon }_{\mathrm{Var},\mathrm{rel}}=\frac{\left|{\overline{\epsilon }}_{\sigma ,\mathrm{KL}}-{\overline{\epsilon }}_{\sigma ,\text{analytic}}\right|}{{\overline{\epsilon }}_{\sigma ,\text{analytic}}}.$$E38

The analytical mean error variance ${\overline{\epsilon }}_{\sigma ,\text{analytic}}$ cannot be computed for squared exponential type autocorrelation coefficient functions (Eq. (35)) and is calculated numerically using the fine mesh containing 3943 nodes, as shown in Figure 1. The relative error is computed for the seven different meshes mentioned in Section 4.2 and plotted in a log–log plot in Figure 5. A closer look at Eq. (10) indicates that the mean error variance is dependent only on the cumulative sum of the $K_1$ eigenvalues. As the eigenvalues are arranged in descending order, the relative error of the largest eigenvalue $\left({\lambda }_1\right)$ is also shown in Figure 5. The corresponding relative error measure can be written in a generalized manner as:

Figure 5.

Convergence of relative errors for different variables for a fixed number of terms $K_1=12$ in the K-L expansion.

$${\epsilon }_{{\lambda }_1}=\frac{\left|{\lambda }_1-{{\lambda }_1}_{\text{analytic}}\right|}{{{\lambda }_1}_{\text{analytic}}}.$$E39

Other error measures can be computed in a similar manner just by replacing ${\lambda }_1$ with the relevant variable. Finally, the relative error of the gradient of the largest eigenvalue w.r.t the geometry parameters $\frac{\partial {\lambda }_1}{\partial {}^2\mathbf{\theta }}$ is also shown in Figure 5. To compute the error related to the gradient, one mesh moving step is carried out, where the mesh is moved from the reference domain $\left({}^2{\theta }_1^{\mathit{ref}},{}^2{\theta }_2^{\mathit{ref}}\right)\mathit{=}\left(0.15\mathrm{m},0.05\mathrm{m}\right)$ to an arbitrary domain defined by $\left({}^2{\theta }_1,{}^2{\theta }_2\right)\mathit{=}\left(0.16\mathrm{m},0.04\mathrm{m}\right)$. The same linear shape functions are used for discretization of the random field, the mesh moving method and the solution of the forward problem. The same realization of the parameter vector ${}^1\mathbf{\theta }$ is used to generate the random field in all the cases. A least squares fit of the relative error plots with a first-order polynomial reveals a linear convergence rate between −2.5 and − 5. Linear rates of convergence for ${\epsilon }_{\mathrm{Var},\mathrm{rel}}$, using linear FEM, are also observed by Betz et al. [14].

Finally, the rate of convergence of the gradient of the potential energy function w.r.t the spatial parameter related to the largest eigenvalue ${\theta }_1$, and the two geometry parameters ${\theta }_{13}$ and ${\theta }_{14}$ is plotted in Figure 6. The convergence rates of the different gradients are almost parallel to each other. Once again, the slope of each plot is determined by a least squares fit with a first-order polynomial. Each plot shows a linear convergence rate with a slope of approximately −2.7.

Figure 6.

Convergence of relative errors of potential energy gradients for a fixed number of terms $K_1=12$ in the K-L expansion.

5. Conclusions

This paper details the numerical procedure for the computation of gradients in probabilistic inverse problems involving the simultaneous estimation of spatial fields and geometry. The method is analytical (in the sense of [6]), intrusive and involves the computation of gradients of the forward problem. Emphasis is laid upon the calculation of gradient of eigenvalues and eigenvectors involved with the truncated K-L expansion method for discretization of the random field. The eigenvalues and eigenvectors are obtained by solving a generalized eigenvalue problem on a defined domain. This implies that as the geometry parameters are updated, the domain is updated and the generalized eigenvalues and eigenvectors change. Computation of the gradient of the eigenvectors w.r.t the geometry parameters involve the computation of a generalized inverse.

The gradients are validated through an inverse analysis using HMC. The potential energy decreases rapidly as the Markov chain related to the geometry parameters approaches the region of high posterior probability, indicating the correctness of the computed gradients. Overall the rate of convergence of various quantities is observed to be linear on a log–log plot. This means computation costs can grow exponentially to achieve better results. While the same mesh that is used to solve the forward problem can also be used for the Galerkin FE method-based discretization of the IEVP, the repeated need to solve the forward problem, the mesh moving method and the generalized eigenvalue problem at every step can be computationally prohibitive. Elimination of any one, two or even all three of these numerical problems can significantly improve the computational feasibility of simultaneous spatial field and geometry detection using gradient-based stochastic samplers like HMC.

Acknowledgments

This work was supported by JSPS KAKENHI Grant Number JP21H02304.

Conflict of interest

The authors declare no conflict of interest.