Inversion of convection–diffusion equation with discrete sources

Abstract

We present a convection–diffusion inverse problem that aims to identify an unknown number of sources and their locations. We model the sources using a binary function, and we show that the inverse problem can be formulated as a large-scale mixed-integer nonlinear optimization problem. We show empirically that current state-of-the-art mixed-integer solvers cannot solve this problem and that applying simple rounding heuristics to solutions of the relaxed problem can fail to identify the correct number and location of the sources. We develop two new rounding heuristics that exploit the value and a physical interpretation of the continuous relaxation solution, and we apply a steepest-descent improvement heuristic to obtain satisfactory solutions to both two- and three-dimensional inverse problems. We also provide the code used in our numerical experiments in open-source format.

Appendices

Appendix 1: Finite-element discretization of the source inversion problem

Here, we briefly review how we discretize the variational problem (1) on a computational mesh. We employ the first-discretize-then-optimize approach, which is a common strategy in PDE-constrained optimization; see; for example, (Gunzburger 2003, Sec. 2.9).

We begin by discretizing the PDE-constraint in (1) following the usual procedure of finite-element methods. To obtain a weak form of the constraint, we multiply both sides of the PDE constraint in (1) with a test function $\phi \in H^1(\varOmega ,\mathbb R)$, the Sobolev space of all functions whose first derivative is square integrable over $\varOmega$, with $\phi (x) = 0$ for $x\in \varGamma _D$, $\frac{\partial \phi (x)}{\partial n} = 0$ for $x\in \varGamma _N$, and then apply Green’s identity and integration by parts:

$$\begin{aligned} \int _{\varOmega } w \phi dx&= \int _{\varOmega } (- c \varDelta u + v^\top \nabla u) \phi dx \end{aligned}$$

(17)

$$\begin{aligned}&= \int _{\varOmega } c (\nabla u)^\top \nabla \phi + (v^\top \nabla u) \phi dx + \int _{\partial \varOmega } \phi \left( - c \frac{\partial u}{\partial n} \right) ds \end{aligned}$$

(18)

$$\begin{aligned}&= \int _{\varOmega } c (\nabla u)^\top \nabla \phi + (v^\top \nabla u) \phi dx, \end{aligned}$$

(19)

where the last identity is obtained by using the Neumann boundary conditions on $\varGamma _N$ and the fact that $\phi (x)=0$ on $\varGamma _D$. The weak problem then is to find a $u \in H^1(\varOmega ,\mathbb R)$ that satisfies (19) for all test functions $\phi \in H^1(\varOmega ,\mathbb R)$.

Next, we discretize the state u and the control w in the weak form of the PDE constraint (19) using compactly supported Ansatz functions on a computational mesh. For ease of presentation we assume that our computational domain $\varOmega = [0,1]^d$ is divided into $N^d$ quadrilateral finite elements of edge length $L=1/N$. We assume a lexicographical ordering of the elements $\varOmega _1, \ldots , \varOmega _{N^d}$ in the mesh, which in our case consists of pixels and voxels for $d=2,3$, respectively. We note that extensions to more general domains and anisotropic or unstructured meshes are straightforward. For example, the implementation used in our numerical experiments uses rectangular meshes whose elements have different edge lengths along the coordinate direction. We use the standard bi-/trilinear Ansatz functions $\phi _1, \phi _2, \ldots ,\phi _{(N+1)^d}$ for $d=2$ and $d=3$, respectively, as test functions and to discretize the state variable u:

$$\begin{aligned} u(x) = \sum _{i=1}^{(N+1)^d} \mathbf {u}_i \phi _i(x). \end{aligned}$$

(20)

Because the Ansatz functions form a Lagrange basis (i.e., for the jth node of the mesh we have $\phi _i(x_j) = \delta _{ij}$), the elements in $\mathbf {u}$ correspond to the value of u at the nodes of the mesh. We represent the control variable w as a linear combination of piecewise constant Ansatz functions $\psi _1, \psi _2, \ldots , \psi _{N^d}$:

$$\begin{aligned} w(x) = \sum _{i=1}^{N^d} \mathbf {w}_i \psi _i(x). \end{aligned}$$

(21)

Similar to above, the Ansatz functions form a Lagrange basis (i.e., $\psi _i(x)=1$ if $x\in \varOmega _i$ and $\psi (x)=0$ else); hence the entries in $\mathbf {w}$ correspond to the values of w in the pixels/voxels of our mesh.

The finite-dimensional approximations of u and w in (20) and (21) allow us to write the weak form of the PDE-constraint (19) in terms of the coefficient vectors $\mathbf {u}$ and $\mathbf {w}$ as

$$\begin{aligned} \mathbf {S}\mathbf {u}= \mathbf {M}\mathbf {w}- \mathbf {r}, \end{aligned}$$

(22)

where $\mathbf {S}\in \mathbb R^{(N+1)^d \times (N+1)^d}$ is the (nonsingular) stiffness and $\mathbf {M}\in \mathbb R^{(N+1)^d \times N^d}$ is the (full-rank) mass matrix, respectively. We compute the entries of both matrices, $\mathbf {M}$ and $\mathbf {S}$, approximately by applying a 5th-order Gaussian quadrature rule to (19):

$$\begin{aligned} \mathbf {S}_{ij}&\approx \int _{\varOmega } c (\nabla \phi _j)^\top \nabla \phi _i - \phi _j v^\top \nabla \phi _i dx \end{aligned}$$

(23)

$$\begin{aligned} \mathbf {M}_{ij}&\approx \int _{\varOmega } \psi _j \phi _i dx.\end{aligned}$$

(24)

$$\begin{aligned} \mathbf {r}_j&\approx \int _{\varGamma _D} \phi _j g (v^\top n) ds. \end{aligned}$$

(25)

Because of the compact support of the basis functions, the integrands vanish on all but a few elements, which leads to both matrices being sparse.

We now derive our discretization of the total variation regularizer (3). In this work, we use a first-order finite-difference discretization of the regularizer. This choice is motivated by its simplicity and computational efficiency. We refer to the recent work in  Herrmann et al. (2018) for a more extensive discussion and finite-element discretizations of the total variation. Because we assume a quadrilateral mesh, the discrete gradient operators for $d=2$ and $d=3$ are

$$\begin{aligned} \mathbf {G}= \left( \begin{array}{r} \mathbf {I}_N \otimes \mathbf {D}\\ \mathbf {D}\otimes \mathbf {I}_N \\ \end{array}\right) \quad \text { and } \quad \mathbf {G}= \left( \begin{array}{r} \mathbf {I}_N \otimes \mathbf {I}_N \otimes \mathbf {D}\\ \mathbf {I}_N \otimes \mathbf {D}\otimes \mathbf {I}_N \\ \mathbf {D}\otimes \mathbf {I}_N \otimes \mathbf {I}_N \\ \end{array}\right) . \end{aligned}$$

(26)

Here, $\mathbf {I}_N\in \mathbb R^{N\times N}$ is the identity matrix, $\otimes$ denotes the Kronecker product, and the one-dimensional difference operator is

$$\begin{aligned} \mathbf {D}= \frac{1}{L} \left( \begin{array}{rrrrrr} -1 &{}\quad 0 &{}\quad 0 &{}\quad \cdots &{}\quad \cdots &{}\quad 0 \\ -1 &{}\quad 1 &{}\quad 0 &{}\quad \cdots &{}\quad \cdots &{}\quad 0 \\ 0 &{}\quad -1 &{}\quad 1 &{}\quad 0 &{}\quad \cdots &{}\quad 0 \\ \vdots &{}\quad &{}\quad \ddots &{}\quad \ddots &{}\quad &{}\quad \vdots \\ \vdots &{}\quad &{}\quad &{}\quad \ddots &{}\quad \ddots &{}\quad 0 \\ 0 &{}\quad \cdots &{}\quad \cdots &{}\quad 0&{}\quad -1 &{}\quad 1 \\ 0 &{}\quad \cdots &{}\quad \cdots &{}\quad 0&{}\quad 0 &{}\quad 1 \\ \end{array} \right) \in \mathbb R^{N+1 \times N}. \end{aligned}$$

We note that this discretization involves a grid change (from the N cell-centers to the $N+1$ nodes). To accommodate for the grid change, the relaxed form includes the average matrix $\mathbf {A}$

$$\begin{aligned} \mathbf {A}= \left( \mathbf {I}_N \otimes \mathbf {B}\,|\,\mathbf {B}\otimes \mathbf {I}_N \right) \end{aligned}$$

for $d=2$ and

$$\begin{aligned} \mathbf {A}= \left( \mathbf {I}_N \otimes \mathbf {I}_N \otimes \mathbf {B}\,|\,\mathbf {I}_N \otimes \mathbf {B}\otimes \mathbf {I}_N\,|\,\mathbf {B}\otimes \mathbf {I}_N \otimes \mathbf {I}_N \right) , \end{aligned}$$

for $d=3$, respectively, where we use the one-dimensional average matrix

$$\begin{aligned} \mathbf {B}= \frac{1}{2} \left( \begin{array}{rrrrrr} 2 &{}\quad 0 &{}\quad &{}\quad \cdots &{}\quad \cdots &{}\quad 0 \\ 1 &{}\quad 1 &{}\quad 0 &{}\quad \cdots &{}\quad \cdots &{}\quad 0 \\ 0 &{}\quad 1 &{}\quad 1 &{}\quad 0 &{}\quad \cdots &{}\quad 0 \\ \vdots &{}\quad &{}\quad \ddots &{}\quad \ddots &{}\quad &{}\quad \vdots \\ \vdots &{}\quad &{}\quad &{}\quad \ddots &{}\quad \ddots &{}\quad 0 \\ 0 &{}\quad \cdots &{}\quad \cdots &{}\quad 0&{}\quad 1 &{}\quad 1 \\ 0 &{}\quad \cdots &{}\quad \cdots &{}\quad &{}\quad 0 &{}\quad 1 \\ \end{array} \right) \in \mathbb R^{N \times N+1}. \end{aligned}$$

This leads to the discrete regularization function

$$\begin{aligned} R_{\mathrm {iso}}(\mathbf {w}) = L^d \mathbf {e}^\top \sqrt{ \mathbf {A}(\mathbf {G}\mathbf {w})^2 + \kappa }, \end{aligned}$$

(27)

where $\mathbf {e}$ is a vector of all ones, the factor $L^d$ represents the volume of the cells, and $\kappa >0$ is a conditioning parameter (in our experiments, we use $\kappa =10^{-3}$). Note that the square and the square root are applied component-wise.

Appendix 2: Finite-difference discretization of the source inversion problem

We discretize both the source, w, and the state, u, in the cell-centered points of our $N_x \times N_y$ computational mesh:

$$\begin{aligned} \mathbf {W}_{kl} \simeq w\left( k L_x - L_x/2, l L_y - L_y/2\right) , \quad \mathbf {U}_{ij} \simeq u\left( iL_x - L_x/2, jL_y - L_y/2\right) , \end{aligned}$$

where $L_x=2/N_x$ and $L_y=1/N_y$ are the discretization steps in the x and y direction, respectively, and $k=1,\ldots ,N_x, l=1,\ldots ,N_y$, $i=0,\ldots ,N_x+1, j=0,\ldots ,N_y$. Here, the variables $\mathbf {U}_{i,0}, \mathbf {U}_{N_x+1,j}, \mathbf {U}_{i,N_y+1}$ approximate the PDE solution at ghost points placed along $\varGamma _N$, and the variables $\mathbf {U}_{0,j}$ are the ghost points near $\varGamma _D$.

Let us further define $\mathbf {V}\in \mathbb R^{m}$ to obtain the bilinear interpolation from the cell-centered points of the mesh closest to the receiver locations $r^1, r^2, \ldots , r^m$. Mathematically, for each receiver location $r^k=(r^k_x, r^k _y), k=1,\ldots , m$, we define variable $V_k$ as

$$\begin{aligned} V_{k} = \dfrac{1}{L_xL_y} \begin{pmatrix} iL_x+\dfrac{L_x}{2}-r^k _x \\ r^k _x-iL_x+\dfrac{L_x}{2} \end{pmatrix}^T \begin{pmatrix} \mathbf {U}_{i,j} &{} \mathbf {U}_{i,j+1} \\ \mathbf {U}_{i+1,j} &{} \mathbf {U}_{i+1,j+1} \end{pmatrix} \begin{pmatrix} jL_y+\dfrac{L_y}{2}-r^k _y\\ r^k _y-jL_y+\dfrac{L_y}{2} \end{pmatrix}, \end{aligned}$$

(28)

where, $i = \lfloor \dfrac{r^k_x}{L_x} + 0.5 \rfloor$ and $j = \lfloor \dfrac{r^k_y}{L_y} + 0.5 \rfloor$, leading to the following mixed-integer quadratic program:

$$\begin{aligned}&\min _{\mathbf {U},\mathbf {W}}\dfrac{1}{2\sigma }\Big (\sum _{k=1}^{m} \Big (\mathbf {V}_{k}-b_k\Big )^2\Big ) \nonumber \\&\quad+ \alpha {L_x L_y} {\sum _{k=2}^{N_x}\sum _{l=2}^{N_y} \sqrt{{\left( \frac{\mathbf {W}_{kl}-\mathbf {W}_{(k-1)l}}{L_x}\right) }^2 + {\left( \frac{\mathbf {W}_{kl}-\mathbf {W}_{k(l-1)}}{L_y}\right) }^2 + \kappa }}, \nonumber \\&\text{ s.t. }&c \frac{4 \mathbf {U}_{ij} - \mathbf {U}_{(i-1)j} - \mathbf {U}_{(i+1)j} - \mathbf {U}_{i(j-1)}-\mathbf {U}_{i(j+1)}}{{L_x}^2 {L_y}^2} \nonumber \\&\quad+ \frac{\mathbf {U}_{ij} - \mathbf {U}_{(i-1)j} }{L_x} =\mathbf {W}_{ij}, \; i=1,\ldots ,N_x, \quad j=1,\ldots ,N_y \nonumber \\&\quad\mathbf {U}_{N_x+1,j}= \mathbf {U}_{N_x,j}, \; \mathbf {U}_{i,0}= \mathbf {U}_{i,1}, \; \mathbf {U}_{i,N_y+1}= \mathbf {U}_{i,N_y}, \mathbf {U}_{0,j}=-\mathbf {U}_{1,j}, \nonumber \\&\quad i=0,\ldots ,N_x,\quad j=0,\ldots ,N_y \nonumber \\&\quad \mathbf {W}\in \{0,1\}^{N_x\times N_y}, \mathbf {U}\in \mathbb R^{(N_x+2)\times (N_y+2)} . \end{aligned}$$

(29)

Here, the fifth row explicitly encodes the Neumann and Dirichlet boundary conditions. As before, we can again eliminate the state variables $\mathbf {U}$ using the discretized PDE and boundary conditions and the state variables $\mathbf {V}$ using (28), resulting in a problem that has similar structure to (6). As before, the objective function is second-order cone representable.

Appendix 3: Selection of regularization parameter for the relaxed problem

To find an effective regularization parameter, we consider the continuous relaxation (7) and follow the L-curve procedure; see (Hansen 1998) for details. In the inversions we use coarser meshes with $256\times 128$ and $96\times 48\times 48$ cells for the 2D instance and 3D instance, respectively. We consider the datasets generated in the preceding section and perturb the generated data with 10% iid Gaussian white noise.

To compute the L-curve, we solve 30 instances of the continuous relaxation for different values of $\alpha$ that are logarithmically spaced between 1 and $10^{-6}$. To accelerate the computation, we initialize the optimization with the solution from the previous $\alpha$. For each value of $\alpha$, we use up to 20 Gauss-Newton iterations and approximately compute the search direction using 5 iterations of projected preconditioned CG that use the Hessian of the regularization function as a preconditioner. For each experiment, we store the reconstructed source, the predicted data, the value of the misfit function, and the values of the regularization function (without the factor $\alpha$). The L-curve shown in Fig. 9 show the value of the regularizer and the value of the misfit of these optimal solutions. As is common, the axes are scaled logarithmically; and to provide additional insight, we have added visualizations of the reconstructed sources for the largest and smallest value of $\alpha$ (resulting in overly smoothed and very noisy reconstructions, respectively) as well as solutions that provide a good trade-off. Using this process we select the regularization parameters $\alpha =8.531 \times 10^{-3}$ for the two-dimensional instance and $\alpha =5.298\times 10^{-3}$ for the three-dimensional instance, respectively. Computing the L-curves took about 4 and 48 minutes and involved 32,580 and 30,892 PDE solves in 2D and 3D, respectively. The large number of PDE solves underscores the importance of computing a factorization (or a good preconditioner in large-scale problems) apriori.

Fig. 9

L-curve plots for the relaxed optimization problem (7) for the two-dimensional instance (left) and three-dimensional instance (right). In both cases, we solve the relaxed problem for 30 $\alpha$ values logarithmically spaced between 1 and $10^{-6}$. We plot the value of the regularizer and misfit at the computed solution in a loglog plot. To highlight the impact of $\alpha$ on the smoothness of the reconstructed images, we provide snapshots of the reconstructed source at the extremal values and one value that provides a good trade-off (values are marked with a circle)