On extension of the data driven ROM inverse scattering framework to partially nonreciprocal arrays

1 Introduction

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 $u_{p}$ (corresponding to unknown coefficient $p$ ) introduced in [BoDrMaMoZa] were used in the Lippmann-Schwinger integral equation framework. If we consider the data $F_{p}$ for the unknown coefficient $p$ and background data $F_{0}$ corresponding to a known background, then the Lippmann-Schwinger integral equation with respect to the unknown $p$ can be written as

F_{p} - F_{0} = - ⟨ u_{0}, p u_{p} ⟩

(1)

where $⟨, ⟩$ is the continuous $L^{2}$ inner product on the PDE domain, and where $u_{p}$ and $u_{0}$ are the unknown and background internal solutions respectively. For the Lippmann-Schwinger Lanczos integral equation (LSL IE), we use the data generated internal solution $u_{p}$ in place of $u_{p}$ :

F_{p} - F_{0} \approx - ⟨ u_{0}, p u_{p} ⟩ .

(2)

Since $u_{p}$ is precomputed directly from the data without knowing $p$ , (2) becomes linear with respect to $p$ .

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.

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} . & . & . & . & . . & . & . & . . & . & . . & . . \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ & \begin{matrix} . & . . . & . . & . . & . \end{matrix} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

(3)

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

- \frac{d^{2} u}{d x^{2}} (x, λ) + p (x) u (x, λ) + λ u (x, λ) = g^{(1)} (x) \frac{d u}{d x} |_{x = 0} = 0, \frac{d u}{d x} |_{x = L} = 0,

(4)

where $0 < L \leq \infty$ . The source $g^{(1)} (x)$ is assumed to be a compactly supported real distribution localized near the origin, for example, roughly speaking, $g^{(1)} = δ (x - ϵ)$ with small $ϵ > 0$ . We also use a similar distribution $g^{(2)}$ near the right boundary $x = L$ to represent the second receiver. Consider $λ_{j} \in C ∖ R_{-}$ , $j = 1, \dots, m$ , since for nonnegative $p$ 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

	$F (λ)$	$=$	$[\int_{0}^{L} g^{(1)} (x) u (x, λ) d x, \int_{0}^{L} g^{(2)} (x) u (x, λ) d x]$		(5)
		$=$	$[⟨ g^{(1)}, u ⟩, ⟨ g^{(2)}, u ⟩]$		(6)

where throughout the paper we use $⟨, ⟩$ to denote the continuous Hermitian $L^{2}$ inner product

⟨ w, v ⟩ = \int_{0}^{L} ¯ w (x) v (x) d x,

which in the one dimensional case is $L^{2} (0, L)$ .

For simplicity of exposition we consider $2 m$ real spectral data points that is, we consider the data

F (λ) |_{λ = λ_{j}} \in R, \frac{d F (λ)}{d λ} |_{λ = λ_{j}} \in R for j = 1, \dots, m .

(7)

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 $p (x)$ 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,

F^{σ} := (F (λ))_{1}

is the symmetric portion of $F$ , 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) , ${u (x, λ_{j})}$ , form a basis for the projection subspace

U = span {u_{1} (x) = u (x, λ_{1}), \dots, u_{m} (x) = u (x, λ_{m})} .

We define the data-driven ROM as the Galerkin system corresponding to $U$

(S + λ M) c = b

(8)

where $S, M \in R^{m \times m}$ are symmetric positive definite matrices with the stiffness matrix $S$ given by

S_{i j} = ⟨ u_{i}^{'}, u_{j}^{'} ⟩ + ⟨ p u_{i}, u_{j} ⟩

and mass matrix $M$ given by

M_{i j} = ⟨ u_{i}, u_{j} ⟩ .

The right hand side $b \in R^{m}$

is a column vector with components

b_{j} = ⟨ u_{j}, g ⟩,

and the Galerkin solution for the system is determined by the vector $c \in R^{m}$ depending on $λ$ . Note that $c = c (λ)$ corresponds to a column vector of coefficients of the solution with respect to the above basis of exact solutions. The matrices $S$ and $M$ can be obtained directly from the data, without knowing the exact solutions, from the formulas

M_{i j} = \frac{F^{σ} (λ_{i}) - F^{σ} (λ_{j})}{λ_{j} - λ_{i}}, M_{i i} = - \frac{d F^{σ}}{d λ} (λ_{i}) .

(9)

and

S_{i j} = \frac{F^{σ} (λ_{j}) λ_{j} - F^{σ} (λ_{i}) λ_{i}}{λ_{j} - λ_{i}}, S_{i i} = \frac{d (λ F^{σ})}{d λ} (λ_{i}) .

(10)

Due to the matching conditions, the ROM transfer function corresponding to this Galerkin system

{~ F}^{σ} (λ) = b^{⊤} c

matches (the symmetric part of) the data exactly. This is well known; a proof of this for real $λ_{i}$ is in [BoDrMaMoZa], and for complex $λ$ in [DrMoZa]. Furthermore, for any $λ$ , the solution to (4) is close to its Galerkin projection

u (λ) \approx ~ u (λ) = V c = V (S + λ M)^{- 1} b

where $V$ represents the row vector of basis functions $u_{i}$ ,

V = [u_{1}, \dots, u_{m}] .

Next we perform a change of basis, that is, we orthogonalize by using the Lanczos algorithm. More precisely, we run $m$ steps of the $M$ -symmetric Lanczos algorithm corresponding to matrix $A = M^{- 1} S$ and initial vector $M^{- 1} b$ . This yields tridiagonal matrix $T \in R^{m \times m}$ and $M$ -orthonormal Lanczos vectors $q_{i} \in C^{m}$ , such that

A Q = Q T, Q^{⊤} M Q = I,

(11)

where

Q = [q_{1}, q_{2}, \dots, q_{m}] \in C^{m \times m},

and

q_{1} = M^{- 1} b / \sqrt{b^{⊤} M^{- 1} b} .

The new basis is orthonormal in $L^{2} (0, L)$ and is given by the vector

V Q = [m \sum j = 1 q_{j 1} u_{j}, \dots, m \sum j = 1 q_{j m} u_{j}] .

The Galerkin solutions and transfer function can then be written in this new basis as

~ u (λ) = \sqrt{b^{⊤} M^{- 1} b} V Q (T + λ I)^{- 1} e_{1},

(12)

{~ F}^{σ} (λ) = (b^{⊤} M^{- 1} b) e_{1}^{⊤} (T + λ I)^{- 1} e_{1}

(13)

where $e_{1} = (1, 0, \dots, 0)^{T}$ is the first coordinate column vector in $R^{m}$ .

Now, we use this orthogonalized ROM to produce internal solutions. Of course, we do not know $p$ , so we don’t know the original basis of exact solutions $V$ . What we do instead is to replace the unknown orthogonalized internal solutions $V Q$ with orthogonalized background solutions $V_{0} Q_{0}$ corresponding to background $p_{0} = 0$ . Here $V_{0}$ is the row vector of background solutions

V_{0} = [u_{1}^{0}, \dots, u_{m}^{0}]

to (4) corresponding to $p = p_{0} = 0$ and the same spectral points $λ = λ_{1}, \dots λ_{m}$ . A ROM for this background problem is computed in the same way, and $Q_{0}$ is computed from its Lanczos orthogonalization. That is, one can compute an approximation $u$ to the unknown internal solution $u (x, λ)$ using

u \approx u = \sqrt{b^{⊤} M^{- 1} b} V_{0} Q_{0} (T + λ I)^{- 1} e_{1}

(14)

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

F_{0} (λ_{j}) - F (λ_{j}) = [\int u_{0}^{(1)} (x, λ_{j}) u (x, λ_{j}) p (x) d x, \int u_{0}^{(2)} (x, λ_{j}) u (x, λ_{j}) p (x) d x]

(15)

for $j = 1, \dots, m$ where $u_{0}^{(i)}$ is the solution to the background problem

- \frac{d^{2} u_{0}^{(i)}}{d x^{2}} (x, λ) + λ u_{0}^{(i)} (x, λ) = g^{(i)} (x)

(16)

with zero Neumann conditions. Correspondingly, $F_{0}$ is the background transfer function

F_{0} = [⟨ g^{(1)}, u_{0}^{(1)} ⟩, ⟨ g^{(2)}, u_{0}^{(1)} ⟩] .

(17)

For real $λ_{j} \in R$ we then have $2 \times 2 m$ equations

F_{0} (λ_{j}) - F (λ_{j}) = [\int u_{0}^{(1)} (x, λ_{j}) u (x, λ_{j}) p (x) d x, \int u_{0}^{(2)} (x, λ_{j}) u (x, λ_{j}) p (x) d x],

(18)

\frac{d}{d λ} (F_{0} - F) |_{λ = λ_{j}} = [\int \frac{d}{d λ} [u_{0}^{(1)} (x, λ) u (x, λ)]_{λ = λ_{j}} p (x) d x, \int \frac{d}{d λ} u_{0}^{(2)} (x, λ) u (x, λ)]_{λ = λ_{j}} p (x) d x]

(19)

for $j = 1, \dots, m$ . If we put the background solutions in the vector

U_{0} = [u_{0}^{(1)}, u_{0}^{(2)}]

we can write (18) and (19) as

F_{0} (λ_{j}) - F (λ_{j}) = \int U_{0} u (x, λ_{j}) p (x) d x,

(20)

(21)

for $j = 1, \dots, m$ . As usual, the internal solutions $u (x, λ_{j})$ and their derivatives with respect to $λ$ are unknown, and depend on $p$ , so the system (15-18) is nonlinear with respect to $p$ . Now, just as in [DrMoZa] we replace $u (x, λ)$ in (15-18) with its approximation

u \approx u = \sqrt{b^{*} M^{- 1} b} V_{0} Q_{0} (T + λ I)^{- 1} e_{1},

which we obtained from the ROM constructed from the symmetric part of the data. We can write the new system for $p$ as

δ F = \int W (x) p (x) d x

(22)

where

δ F = [(F_{0} - F) (λ_{1}), \dots, (F_{0} - F) (λ_{m}), \frac{d}{d λ} (F_{0} - F) (λ_{1}), \dots, \frac{d}{d λ} (F_{0} - F) (λ_{m})] \in R^{4 m},

and

W = [u U_{0} (λ_{1}), \dots, u U_{0} (λ_{m}), \frac{d}{d λ} (u U_{0}) |_{λ = λ_{1}}, \dots, \frac{d}{d λ} (u U_{0}) |_{λ = λ_{m}}]

is a $4 m$ -dimensional vector of functions on $(0, L)$ . Recall $u$ is computed directly from the symmetric part of the data without knowing $p$ , 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 $1 \times 2$ 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 $Ω \in R^{d}$ :

- Δ u^{(r)} + p u^{(r)} + λ u = g^{(r)}, \frac{d u}{d ν} |_{\partial Ω} = 0, r = 1, \dots, K,

(23)

where $g^{(r)}$ are $K$ localized sources, e.g., boundary charge distributions, supported near or at an accessible part $S$ of $\partial Ω$ . Consider also a set of distributions representing our receivers, which will include all of the sources distributions. That is, let

G = [g^{(1)}, g^{(2)}, \dots, g^{(K)}, \dots, g^{(L)}],

where $L \geq K$ . Let

U = [u^{(1)}, u^{(2)}, \dots, u^{(K)}]

be our solution functions corresponding to the $K$ sources. Then the multiple-input multiple output (MIMO) transfer function is a $K \times L$ matrix valued function of $λ$

F (λ) = ⟨ G, U ⟩ \in R^{K \times L}

(24)

where $⟨, ⟩$ represents the continuous $L^{2} (Ω)$ inner product, that is,

⟨ G, U ⟩ = \int_{Ω} U^{⊤} G d x,

is matrix valued with

(F (λ))_{i j} = \int_{Ω} u^{(i)} g^{(j)} d x .

We consider the inverse problem with data given by $2 m$ real symmetric $K \times L$ matrices, that is,

F (λ) |_{λ = λ_{j}} \in R^{K \times L},

and

\frac{F (λ)}{d λ} |_{λ = λ_{j}} \in R^{K \times L},

for real spectral points ${λ_{j}}_{j = 1, \dots m}$ , 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:

F^{σ} (λ) := (F (λ)) |_{i, j = 1 \dots K}

and correspondingly define

G^{σ} = [g^{(1)}, g^{(2)}, \dots, g^{(K)}] .

We consider the $m K$ dimensional projection subspace

U = span {U_{1} (x) = U (x, λ_{1}), \dots, U_{m} (x) = U (x, λ_{m})}

and define the MIMO data-driven ROM as the system

(S + λ M) C = B

(25)

where $S, M \in R^{m K \times m K}$ are symmetric positive definite matrices, $B \in R^{m K \times K}$ is the symmetric part of the data

B = ⟨ G^{σ}, U ⟩,

and the system solution $C = C (λ) \in R^{m K \times K}$ 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

S = (S_{i j} = ⟨ {U^{'}}_{i}, U_{j}^{'} ⟩) + ⟨ p U_{i}, U_{j} ⟩)

and

M = (M_{i j} = ⟨ U_{i}, U_{j} ⟩) .

Again $S$ and $M$ are obtained by imposing the conditions that the ROM transfer function ${~ F}^{σ} (λ) = B^{⊤} C$ matches the data. Next, we perform block Lanczos tridiagonalization, with matrix

A = M^{- 1} S

and initial block vector $M^{- 1} B$ . From this we obtain the block-tridiagonal matrix

T \in R^{K m \times K m}

with $K \times K$ blocks, and block-vectors $q_{i} \in R^{m K \times K}$ which are orthonormal with respect to the $M$ inner product, which again corresponds to continuous $L^{2}$ orthonormality. From this we obtain the block counterpart of (11)

Q = [q_{1}, q_{2}, \dots, q_{m}] \in R^{m K \times m K},

where

q_{1} = M^{- 1} B (B^{⊤} M^{- 1} B)^{- 1 / 2} .

The Galerkin projection of the true solution or state solution can then be written in the new Lanczos basis as

~ U (λ) = \sqrt{B^{⊤} M^{- 1} B} V Q (T + λ I)^{- 1} E_{1},

(26)

where $E_{1} \in R^{m K \times K}$ consists of the first $K$

columns of identity matrix

I \in R^{m K \times m K}

. We express the orthogonalized basis as

V Q = [m \sum j = 1 U_{j} q_{j 1}, \dots, m \sum j = 1 U_{j} q_{j m}]

where ${q_{j i} \in R^{K \times K}}_{j i}^{K \times K}$ are the blocks of matrix $Q$ and $U_{j} = [u_{j}^{(1)} (x), \dots u_{j}^{(K)} (x)]$ are the corresponding vectors of solutions. Note that this basis depends on the unknown internal solutions, so as before we replace the basis $V Q$ in (26) with its background counterpart $V_{0} Q_{0}$ . This yields the data generated vector of internal solutions

U := \sqrt{B^{⊤} M^{- 1} B} V_{0} Q_{0} (T + λ I)^{- 1} E_{1} .

(27)

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

F_{0} (λ_{j}) - F (λ_{j}) = \int U^{⊤} (x, λ_{j}) U_{0} (x, λ_{j}) p (x) d x, j = 1, \dots, m .

(28)

Here $U_{0}$ represents all $L$ background solutions with $p = 0$ . For $λ_{j} \in R$ we will have

	$(F_{0} - F) \|_{λ_{j}}$	$=$	$\int U^{⊤} (x, λ_{j}) U_{0} (x, λ_{j}) p (x) d x,$
	$\frac{d}{d λ} (F_{0} - F) \|_{λ_{j}}$	$=$	$\int \frac{d}{d λ} (U^{⊤} (x, λ_{j}) U_{0} (x, λ)) \|_{λ = λ_{j}} p (x) d x .$		(29)

Now we replace $U (x, λ)$ with its approximation $U$ in (28) and (3.3). Again, precomputing $U$ from the data only (27) will yield the linear system for $p$ , which is given by the upper triangular part of (30) and (31):

	$(F_{0} - F) \|_{λ_{j}}$	$=$	$\int U^{⊤} (x, λ_{j}) U_{0} (x, λ_{j}) p (x) d x,$		(30)
	$\frac{d}{d λ} (F_{0} - F) \|_{λ_{j}}$	$=$	$\int \frac{d}{d λ} (U^{⊤} (x, λ_{j}) U_{0} (x, λ)) \|_{λ = λ_{j}} p (x) d x .$		(31)

This setup corresponds to the full upper triangular $K \times L$ (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 $U$ , so we can use any sparsity patterns in the right upper blocks of $(F_{0} - F)$ and $\frac{d}{d λ} (F_{0} - F)$ . 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 $N = 100 \times 100$ elements to compute the synthetic data. The LSL system is solved here using the background solutions as a basis for the unknown $p$ . We use $m = 6$ positive spectral values $λ = 1, 2, 14, 50, 128, 262.2672$ . 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.

Figure 1: True one dimensional medium (top left) and its reconstructions using SIMO data and the Born approximation (top right), Lippmann Schwinger Lanczos using SISO data only (bottom left), and Lippmann Schwinger Lanczos using SIMO data (bottom right)

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 $K$ sources $L$ receivers on the top boundary of the domain. We use $m = 5$ positive spectral values $λ = 1, 2, 14, 50, 128$ . The forward problem was discretized on a regular triangular grid using a finite-element method with $N = 451 \times 151$ elements. Equation (30) was approximated using nodal quadrature on a $901 \times 301$

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

K = L = 2

. 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

U (x, λ))

in (3.3) is replaced by with

U_{0} (x, λ))

), 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

U (x, λ))

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 $K = 1$ source and $L = 20$ receivers located on all four sides of square domain (see top left in Fig.4) and compared it to SISO scenario $K = L = 1$ . We discretized the problem using finite elements on a uniform triangular grid with $N = 151 \times 151$ elements. Equation 30 was approximated using quadrature on $451 \times 451$ grid. In this example we considered $m = 6$ spectral values $λ = - 24, 1, 2, 14, 50, 128$ , 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 $p$ with its Born and LSL IE SIMO images in Fig. 5.

Figure 2: Experiment 1: True medium (top left) and its reconstructions using the LSL IE for square data $K = L = 2$ (top right), Born linearization for rectangular data $K = 2, L = 12$ (bottom left) and extended LSL IE for the same rectangular data(bottom right)

Figure 3: Experiment 1: LSL IE image with exact internal solution $U (x, λ))$ (Cheated IE) for square data $K = L = 2$ (left) and for rectangular data $K = 2, L = 12$ (right)

Figure 4: True p with its LSL SIMO, Born SIMO, and LSL SISO images. The SIMO reconstructions are qualitatively correct, while the SISO reconstruction is not. The LSL SIMO (IE) is a quantitatively better approximation to the true $p$ than Born.

Figure 5: Slices of true $p$ , the SIMO Born and the SIMO (LSL) IE. The LSL IE reconstruction is obviously closer to the true model.

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.

On extension of the data driven ROM inverse scattering framework to partially nonreciprocal arrays

Authors

Lippmann-Schwinger-Lanczos algorithm for inverse scattering problems

Reduced Order Model Approach to Inverse Scattering

Data completion algorithms and their applications in inverse acoustic scattering with limited-aperture backscattering data

Windowed Green function method for wave scattering by periodic arrays of 2D obstacles

Imaging based on Compton scattering: model-uncertainty and data-driven reconstruction methods

Reduced order models for spectral domain inversion: Embedding into the continuous problem and generation of internal data

1 Introduction