In this work we extend the reduced order model (ROM) approach which was used previously for inverse impedance, scattering and diffusion [borcea2011resistor, borcea2014model, druskin2016direct, druskin2018nonlinear, borcea2017untangling, borcea2019robust, BoDrMaMoZa, borcea2020reduced, DrMoZa, Borcea2021ReducedOM] to nonsymmetric data sets. Here, we consider the inverse diffusion problem given spectral data, although the technique to handle non-symmetric transfer functions can be applied to a much wider class of problems.
In the ROM framework for solving multidimensional inverse problems, the ROM is constructed from a symmetric data set, that is, for coinciding source and receiver pairs. That is, the ROM is chosen precisely to match this data set, see [druskin2016direct, druskin2018nonlinear, borcea2017untangling, borcea2019robust, BoDrMaMoZa, Borcea2021ReducedOM, borcea2021reduced]. Then, the ROM is transformed to a sparse form (tridiagonal for single input single output (SISO) problems, block tridiagonal for multiple input/output (MIMO) problems) by Lanczos orthogonalization. The data-driven ROM, in this orthogonalized form, can be viewed as a discrete network which has entries for which their dependence on the unknown pde coefficients is approximately linear [borcea2011resistor, druskin2016direct, borcea2014model, borcea2020reduced, Borcea2021ReducedOM]. That is, the main nonlinearity of the inverse problem was absorbed by the orthogonalization process. This process is related to the seminal works of Marchenko, Gelfand, Levitan and Krein on inverse spectral problems, and to the idea of spectrally matched second order staggered finite-difference grids first introduced in [druskin1999gaussian] and first used for direct inversion in [BoDr].
The data-driven ROM can be viewed as a projection operator, or Galerkin system [doi:10.1137/1.9781611974829.ch7]. A crucial property, which was first noticed in [druskin2016direct], is that the Galerkin basis corresponding to the unknown coefficient is very close to the one from the homogeneous problem. This led to the back projection algorithm [druskin2018nonlinear], which allowed for direct reconstructions in multiple dimensions. In [BoDrMaMoZa], it was found that thanks to this weak dependence of the basis functions on the unknown coefficients, the ROM can also be directly used to generate internal solutions from boundary data only.
In the work [DrMoZa], the data-generated internal solutions (corresponding to unknown coefficient ) introduced in [BoDrMaMoZa] were used in the Lippmann-Schwinger integral equation framework. If we consider the data for the unknown coefficient and background data corresponding to a known background, then the Lippmann-Schwinger integral equation with respect to the unknown can be written as
where is the continuous inner product on the PDE domain, and where and are the unknown and background internal solutions respectively. For the Lippmann-Schwinger Lanczos integral equation (LSL IE), we use the data generated internal solution in place of :
Since is precomputed directly from the data without knowing , (2) becomes linear with respect to .
The Lippmann-Schwinger-Lanczos approach resolves directly one of the main limitations of the above mentioned works; that it applies only to symmetric data sets. This is the subject of this manuscript. Indeed, while the ROM must still be constructed from a symmetric subset of the data, the Lippmann-Schwinger Lanczos equation (2) can be directly applied to any additional receivers, using the same internal solutions. Furthermore, additional sources (with receivers in the symmetric part) can also be added thanks to reciprocity. Adding this additional data increases the range of the Lippmann-Schwinger-Lanczos linear operator on the r.h.s. of (2), thus improving the quality of the solution of the inverse problem.
The structure of our extended multi-input/multi-output (MIMO) measurement array can be summarized with the help of (3). The column and row numbers correspond respectively to the indices of the receivers (outputs) and transmitters (inputs). The conventional data-driven ROM requires collocated receiver and transmitter arrays (the symmetric data set), so the measurements in this case are given by a square matrix, that is, the left upper block of (3). We assume reciprocity, implying symmetric transfer functions, thus it is sufficient just to measure their upper triangular parts. The extended receivers are shown in the upper right block. A common case is when all added receivers measure responses from one of the transmitters, resulting in adding a column. However, added receivers can potentially just measure responses from certain transmitters, resulting in an arbitrary sparsity pattern in the upper right block. Likewise, one can add transmitters by adding elements to the lower left block. By reciprocity, they can be symmetrically reflected into the right upper one.
This paper is organized as follows. In Section 2 we describe the entire process in detail for a one dimensional single input multiple output (SIMO) problem. We briefly describe the construction of the ROM from the single input/ single output (SISO) part of data, the Lanczos orthogonalization process and the generation of the internal solutions. We then show its use in the Lippmann-Schwinger equation for the full SIMO data. The generalization of this process to SIMO or nonsymmetric MIMO arrays in higher dimensions is described in Section 3. Section 4 contains numerical experiments.
2 One dimensional SIMO problem
We begin this work with a one dimensional problem to demonstrate the approach. In the first subsection we describe the problem setup one source and two receivers, one of which is at the source. We note the transfer function for this problem has the SISO transfer function as a component, which we use to construct the ROM. We review this briefly in the second subsection, along with tridiagonalization of the ROM and generation of the internal solutions. In the last subsection we show how to use the full transfer function and the internal solutions in the Lippmann-Schwinger equation in order to solve the fully nonlinear SIMO inverse problem.
2.1 Description of the SIMO problem
We start by considering the following single input multiple output (SIMO) inverse problem in one dimension
where . The source is assumed to be a compactly supported real distribution localized near the origin, for example, roughly speaking, with small . We also use a similar distribution near the right boundary to represent the second receiver. Consider , , since for nonnegative the above resolvent is well defined for off of the negative real axis. We note that the above formulation can be obtained via the Laplace transform of the one dimensional diffusion or wave problem. The SIMO transfer function is then
where throughout the paper we use to denote the continuous Hermitian inner product
which in the one dimensional case is .
For simplicity of exposition we consider real spectral data points that is, we consider the data
In this case all solutions are real and the conjugates are unnecessary. For complex data points, all of the following holds, see [DrMoZa]. The SIMO inverse problem is then to determine in (4) from the data (7).
2.2 Construction of the data-driven ROM, orthogonalization and internal solutions
The first component of the transfer function (5) is a SISO transfer function,
is the symmetric portion of , and can be used to construct the ROM exactly as in previous works e.g. [DrMoZa]. We describe this now briefly. The exact solutions to (4) , , form a basis for the projection subspace
We define the data-driven ROM as the Galerkin system corresponding to
where are symmetric positive definite matrices with the stiffness matrix given by
and mass matrix given by
The right hand side
is a column vector with components
and the Galerkin solution for the system is determined by the vector depending on . Note that corresponds to a column vector of coefficients of the solution with respect to the above basis of exact solutions. The matrices and can be obtained directly from the data, without knowing the exact solutions, from the formulas
Due to the matching conditions, the ROM transfer function corresponding to this Galerkin system
matches (the symmetric part of) the data exactly. This is well known; a proof of this for real is in [BoDrMaMoZa], and for complex in [DrMoZa]. Furthermore, for any , the solution to (4) is close to its Galerkin projection
where represents the row vector of basis functions ,
Next we perform a change of basis, that is, we orthogonalize by using the Lanczos algorithm. More precisely, we run steps of the -symmetric Lanczos algorithm corresponding to matrix and initial vector . This yields tridiagonal matrix and -orthonormal Lanczos vectors , such that
The new basis is orthonormal in and is given by the vector
The Galerkin solutions and transfer function can then be written in this new basis as
where is the first coordinate column vector in .
Now, we use this orthogonalized ROM to produce internal solutions. Of course, we do not know , so we don’t know the original basis of exact solutions . What we do instead is to replace the unknown orthogonalized internal solutions with orthogonalized background solutions corresponding to background . Here is the row vector of background solutions
to (4) corresponding to and the same spectral points . A ROM for this background problem is computed in the same way, and is computed from its Lanczos orthogonalization. That is, one can compute an approximation to the unknown internal solution using
which is obtained from data only.
2.3 Nonlinear inverse problem
We now use the Lippmann-Schwinger formulation to solve the nonlinear inverse problem for all of the data, not just the symmetric part, which is the SISO transfer function in this case. From (5) and its background counterpart we obtain the Lippmann-Schwinger equation
for where is the solution to the background problem
with zero Neumann conditions. Correspondingly, is the background transfer function
For real we then have equations
for . If we put the background solutions in the vector
for . As usual, the internal solutions and their derivatives with respect to are unknown, and depend on , so the system (15-18) is nonlinear with respect to . Now, just as in [DrMoZa] we replace in (15-18) with its approximation
which we obtained from the ROM constructed from the symmetric part of the data. We can write the new system for as
is a -dimensional vector of functions on . Recall is computed directly from the symmetric part of the data without knowing , by using (14), making nonlinear system (22) linear. We continue to refer to (22) as a Lippmann-Schwinger-Lanczos system. Setup (22) corresponds to the full array in (3).
3 Multidimensional MIMO problem
In this section we consider a more general multidimensional case where we may have any number of sources and receivers. The symmetric part of the transfer function will consist of all of the data from coinciding source/receiver pairs, that is, it corresponds to the left upper block of (3). Again, the symmetric part will be used to build the ROM as in [DrMoZa], while the remaining data will be used later in the Lippmann-Schwinger formulation. We will assume that the nonsymmetric, or remaining part of the transfer function consists of receivers, that is, they will be located in the right upper block of (3). As it was already mentioned in the introduction, this is without loss of generality, since any remaining sources can equivalently be viewed as receivers by reciprocity.
3.1 Description of the MIMO problem
We consider a formulation that can be obtained via the Laplace transform of the diffusion or wave problem, which can be written as the following boundary value problem on :
where are localized sources, e.g., boundary charge distributions, supported near or at an accessible part of . Consider also a set of distributions representing our receivers, which will include all of the sources distributions. That is, let
where . Let
be our solution functions corresponding to the sources. Then the multiple-input multiple output (MIMO) transfer function is a matrix valued function of
where represents the continuous inner product, that is,
is matrix valued with
We consider the inverse problem with data given by real symmetric matrices, that is,
for real spectral points , while we note again that complex data will follow in the same way [DrMoZa].
3.2 Construction of the data-driven ROM, block Lanczos, and internal solutions
We now construct the ROM based on the symmetric part of the transfer function:
and correspondingly define
We consider the dimensional projection subspace
and define the MIMO data-driven ROM as the system
where are symmetric positive definite matrices, is the symmetric part of the data
and the system solution is a matrix valued function of , again corresponding to coefficients of the solution with respect to the above basis of exact solutions. Stiffness and mass matrices are block versions of (8), with blocks given by
Again and are obtained by imposing the conditions that the ROM transfer function matches the data. Next, we perform block Lanczos tridiagonalization, with matrix
and initial block vector . From this we obtain the block-tridiagonal matrix
with blocks, and block-vectors which are orthonormal with respect to the inner product, which again corresponds to continuous orthonormality. From this we obtain the block counterpart of (11)
The Galerkin projection of the true solution or state solution can then be written in the new Lanczos basis as
where consists of the first
columns of identity matrix. We express the orthogonalized basis as
where are the blocks of matrix and are the corresponding vectors of solutions. Note that this basis depends on the unknown internal solutions, so as before we replace the basis in (26) with its background counterpart . This yields the data generated vector of internal solutions
3.3 Nonlinear MIMO inverse problem
We now solve the fully nonlinear inverse problem using the entire non-symmetric data set. From (24) and the equation for the background solutions we obtain the Lippman-Schwinger formulation
Here represents all background solutions with . For we will have
This setup corresponds to the full upper triangular (3) . However, as already mentioned, removal of any element of the data matrix with the column number larger than K does not affect the computation of , so we can use any sparsity patterns in the right upper blocks of and . Note that these patterns can also vary for different matching frequencies. This would result in the reduction of the range of the forward linear operator without changing the remaining columns.
4 Numerical experiments
4.1 One dimensional SIMO problem
We first demonstrate numerical results here for a one dimensional reconstruction problem as a proof of concept. The setup is exactly as in Section 2, although computations are done on a square, using finite elements and elements to compute the synthetic data. The LSL system is solved here using the background solutions as a basis for the unknown . We use positive spectral values . Note in Figure 1 that the SISO data, which is used to create the internal solutions, reconstructs the bar closer to the source quite well. Using the same internal solutions, we add data from the receiver on the right hand side, and we see that the second bar is captured. The Born approximation does not work well for this high contrast example.
4.2 Two dimensional examples
Here we show our numerical reconstructions of a 2D two-bump media. We assumed that the data is noiseless. The construction of the data-driven ROM may become unstable in the presence of noise, however, one can regularize it using the SVD-based algorithm of [borcea2019robust]. For simplicity, we omitted that part in this paper. We chose frequencies that provide enough sensitivity to recover the unknown scatterers, but so that the data is still non-redundant. In our first experiment, we considered the medium shown on the top left in Fig.2. The data acquisition setup mimics one from surface geophysical exploration, that is, we consider sources receivers on the top boundary of the domain. We use positive spectral values . The forward problem was discretized on a regular triangular grid using a finite-element method with elements. Equation (30) was approximated using nodal quadrature on a
grid. The obtained linear system is ill-conditioned, and we solved it via projecting onto its dominant eigenvectors. On the top right of Fig.2 we plotted the reconstruction using the LSL IE for the square data matrix when . Red crosses and green circles on plots show source(s) and receiver(s) locations, respectively. Then, we expanded the data matrix to the rectangular one by adding 10 receivers. For this scenario, the reconstruction using Born linearization, i.e. when in (3.3) is replaced by with ), is shown on bottom left of Fig.2. For the same data, the image produced by the extended LSL IE is plotted on the bottom right of Fig.2. As one can observe, it improves the image compared to both Born and the square LSL IE. In fact, in Fig.3 we plotted the solutions of (3.3) assuming the exact internal solution is available, which we call ”Cheated IE”. This represents the best possible result if one were to apply an iterative (aka distorted) Lippmann-Schwinger algorithm which converged. (Iterative Lippman-Schwinger would require multiple solutions of forward problems, which can be costly even if such an algorithm converges.) The cheated reconstructions are identical to the corresponding results for the LSL IE in Fig.2, and, obviously, in both cases improvements are obtained compared to the rectangular data just because of the additional receivers.
In the second experiment we considered a SIMO medical imaging setup with source and receivers located on all four sides of square domain (see top left in Fig.4) and compared it to SISO scenario . We discretized the problem using finite elements on a uniform triangular grid with elements. Equation 30 was approximated using quadrature on grid. In this example we considered spectral values , that is, we added a negative one to the set from the previous example. In Fig.4 we plotted the images obtained using the LSL IE for SIMO (top right) and SISO (bottom right) scenarios as well as the Born reconstruction for SIMO (bottom left). In this case the SISO LSL image is qualitatively wrong. Both SIMO LSL IE and born images are qualitatively correct, however the LSL IE image is quantitatively more accurate. To strengthen this argument we compared a slice (the red line in Fig 4 top left ) of the true with its Born and LSL IE SIMO images in Fig. 5.
5 Conclusion and discussion
Using the previously developed Lippmann-Schwinger-Lanczos (LSL) algorithm, we are able to extend successfully the model reduction approach to arrays with nonreciprocal subsets, including SIMO data sets. The main novelty of this algorithm is the exploitation of the product structure of the Lippmann-Schwinger equation. It allows us to plug in an approximate internal solution computed using only the data appropriate for the ROM construction (symmetric) subset of the data, and then using all of the given data to solve the enlarged linear LSL system. Even though all of the derivations here are made for the second order Schrodinger equation in the Laplace domain, and assuming exact data, we plan to extended them to first order systems and the time domain, as well as regularized noisy ROM formulations of [borcea2019robust, borcea2020reduced]. Future directions also include using iterative methods for very sparse data sets. We have already obtained promising preliminary results for the monostatic (SAR) problem, which is an important case of such data sets.
Acknowledgements We are grateful Liliana Borcea, Alex Mamonov and Jörn Zimmerling for productive discussions that inspired this research. V. Druskin was partially supported by AFOSR grant FA 955020-1-0079 and NSF grant DMS-2110773. S. Moskow was partially supported by NSF grants DMS-1715425 and DMS-2008441.