1 Problem formulation

The inverse scattering problem formulated for the Schrödinger operators arises in various fields, including quantum mechanics, radars, viscoelasticity, Biot problems, remote sensing, geophysical, and medical imaging. The goal of imaging is to find medium properties in the domain using near-field measured data. The model based nonlinear optimization which is the method of choice for the solution of the inverse problems can be unreliable and particularly expensive for such problems. Data driven nonlinear transforms can be an opening, however it was recently shown that the ReLU networks are intractable for reliable solution of the inverse problems in continuum using conventional digital computers. In the present work, following the success of data-driven reduced-order models (ROMs) developed in recent years, we propose a robust direct method for solving inverse scattering problems for the Schrödinger equation. Our approach is based on a Lippmann-Schwinger algorithm with a crucial component composed of adaptive data-driven ROMs in the frequency domain and efficient learning the internal solutions. Below we discuss the details of the algorithms as well as some bottlenecks.

We consider first-order formulation of frequency-domain wave problem in lossy medium

\begin{eqnarray} \nabla \cdot v+\frac {1}{2}\nabla (\ln (\sigma ))\cdot v+ru+i\omega u=f\\ \nabla u-\frac {1}{2}\nabla (\ln (\sigma )))u+i\omega v=0 \end{eqnarray}

in a bounded domain \(\Omega \) for \(m\) sources \(f\in \RR ^{\infty \times m}\). Here \(u,v\in \CC ^{\infty \times m}\) with columns being solutions for the corresponding sources. We assume that the measured data is given by \(f^Tu\) where columns of \(f\) and \(u\) are multiplied with respect to \(L_2\) nner product \((g;h)_{L_2}=\int _\Omega {ghdV}\). After spatial discretization we obtain MIMO port-Hamiltonian LTI dynamical system [2] in the frequency-domain

\begin{equation} \label {LTI} (A+P)w+i\omega w=F \end{equation}

with skew-symmetric matrix \(A=-A^T\in \RR ^{N\times N}\), symmetric matrix \(P=P^T\in \RR ^{N\times N}\), \(F\in \RR ^{N\times m}\). The measured data is given by \(D(\omega )=F^Tw\). We note that the obtained system (3) is symmetric with respect to indefinite (pseudo-)inner product \((x;y)_J=x^HJx\) where \(^H\) denotes Hermitian conjugate and

\begin{equation} J=\begin{pmatrix}I&0\\ 0&-I \end {pmatrix} \end{equation}

In the inverse scattering problem in lossy medium the goal is to recover unknown damping term \(r\) and reflectivity \(\ln (\sigma )\) under known data \(D(\omega ))\). In the LTI system (3) the discrete counterparts of both unknowns compose matrix \(P\).

Lippmann-Schwinger approach has been proven to be a powerful tool for solving inverse scattering problem [3]. It can be formulated in terms of integral equation involving solution \((u;v)\) for unknown medium as well as solution \((u^0;v^0)\) for background with some a priori known parameters (say, \(r=\ln (\sigma )=0\)):

\begin{eqnarray} \nabla \cdot v^0+i\omega u^0=f\\ \nabla u^0+i\omega v^0=0. \end{eqnarray}

After spatial discretization we obtain

\begin{equation} \label {LTIbckg} Aw^0+i\omega w^0=F \end{equation}

with background data given by \(D^0(\omega )=F^Tw^0\). The discrete counterpart of Lippmann-Schwinger equation can be formulated then as

\begin{equation} \label {lseq} (w^0)^H(\omega )JPw(\omega )=D^0(\omega )-D(\omega ). \end{equation}

This is a system of nonlinear equations with respect to unknown parameter matrix \(P\) because \(w(\omega )\) itself depends on \(P\). Tradional way to linearize this problem is so-called Born approximation, i.e. \(w(\omega )\approx w^0(\omega )\), however it is known to work for small \(P\) only. Below we will show how to exploit the measured data \(D(\omega )\) to construct a better approximant of \(w(\omega )\) for further linearizarion of Lippmann-Schwinger equation (8).

2 Adaptive data-driven ROMs

Consider Galerkin projection of (3) onto the rational Krylov subspace spanned on columns of \(V=\{w_1=w(\omega _1),\ldots , w_n=w(\omega _n)\}\) with respect to \((;)_J\) inner product. Here \(w(\omega _i),~i=1,\ldots ,n\) are solutions of (3) for \(\omega =\omega _1,\ldots ,\omega _n\), respectively. The projected system (3) has a form

\begin{equation} \label {eq:gal} \cS \cW +i\omega \cM \cW =\cB \end{equation}

where \(w\approx V\cW \), \(\cS \) is Hermitian indefinite stiffness matrix with block elements

\begin{equation} \cS _{pq}=w^H_qJ(A+P)w_p\in \CC ^{m\times m},~q,p=1,\ldots , n \end{equation}

, \(\cM \) is Hermitian indefinite mass matrix with block elements

\begin{equation} \cM _{pq}=w^H_qJw_p\in \CC ^{m\times m},~q,p=1,\ldots , n \end{equation}

and blocks of \(\cB \) are given by

\begin{equation} \cB _q=w^H_qF=\bar {D}_q=\bar {D}(\omega _q)\in \CC ^{m\times m}~q=1,\ldots , n. \end{equation}

The measured data in Galerkin formulation is given by

\begin{equation} \label {romdata}\cF =\cB ^H\cW \end{equation}

We note that although internal solutions \(w_p,~p=1,\ldots n\) are not accessible because \(P\) in (3) is unknown, blocks of mass and stiffness matrices can still be obtained directly from the data via Loewner framework:

\begin{equation} \label {massmtr} w^H_qJw_p=\frac {D_p-\bar {D}_q}{i\omega _q+i\omega _p} \end{equation}

and

\begin{equation} \label {stifmtr} w^H_qJ(A+P)w_p=\frac {\omega _q\bar {D}_q+D_p\omega _p}{\omega _q+\omega _p} \end{equation}

(with derivatives of the data \(D\) and \(\omega D\) for mass and stiffness matrices when \(\omega _p=-\omega _q\), respectively).

To improve the efficiency of the constructed ROM we employ greedy algorithms for adaptive choice of interpolation frequencies \(\omega _1,\ldots ,\omega _n\) which is similar to AAA algorithm [4]:

We note that the obtained ROM is not structure-preserving, i.e. it may not inherit passivity and even stability of the original full-scale system (3).

3 Lippmann-Schwinger-Lanczos approach

In this section we show how to exploit the constructed data-driven ROMs to construct an approximant of internal solution \(w\) in (8). Similar to the unknown medium, we can construct background matrix pencil \((\cS ^0;\cM ^0)\) for the selected set of frequencies \(\omega _1,\ldots ,\omega _n\). Note that it can be performed in model-driven way because all the internal solutions for known background \(P=0\) are accessible. Background counterpart of ROM (9) has a form

\begin{equation} \label {eq:galbckg} \cS ^0\cW ^0+i\omega \cM ^0\cW ^0=\cB ^0 \end{equation}

where \(w^0\approx V^0\cW ^0\) and \(V^0=\{w^0_1=w^0(\omega _1),\ldots , w^0_n=w^0(\omega _n)\}\). Let’s perform Lanczos orthogonalization for matrix \((\cM )^{-1}\cS \) with respect to indefinite inner product \((;)_{\cM }\) and starting vector \((\cM )^{-1}\cB /\|(\cM )^{-1}\cB \|_{\cM }\) and do the same for background part \((\cM ^0)^{-1}\cS ^0\) with respect to \((;)_{\cM ^0}\) and starting vector \((\cM ^0)^{-1}\cB ^0/\|(\cM ^0)^{-1}\cB ^0\|\):

\begin{eqnarray} (\cM )^{-1}\cS Q=QT,~Q^H\cM Q=I\\ (\cM ^0)^{-1}\cS ^0 Q^0=Q^0T^0,~(Q^0)^h\cM ^0 Q^0=I. \end{eqnarray}

As has been noted in [1], although \(V\) is totally different from \(V_0\), we have \(VQ\approx V_0Q_0\). It has been explained in that paper for lossless case via drawing an analogy with causal time-domain solutions, however similar reasoning is applicable for lossy scenario and ROM we developed. Therefore, we can construct an approximant of internal solution as

\begin{equation} w\approx V(\cS +i\omega \cM )^{-1}\cB =VQ(T+i\omega I)^{-1}E_1|(\cM )^{-1}\cB \|_{\cM }\approx V^0Q^0(T+i\omega I)^{-1}E_1|(\cM )^{-1}\cB \|_{\cM }=\bw \end{equation}

Once plugged in into the Lippmann-Schwinger equation (8), we obtain a linear equation with respect to \(P\):

\begin{equation} \label {lsleq} (w^0)^H(\omega )JP\bw (\omega )=D^0(\omega )-D(\omega ). \end{equation}

We call this algorithm Lippmann–Schwinger–Lanczos to emphasize the crucial component of constructing internal solution that is based Lanczos orthogonalization. There are multiple parts of our approach that need to be addressed to improve its performance:

• • Efficient handling of overfitting that results in rank-deficient matrix \(\cM \) and may breakdown Lanczos algorithm

• • Fast and robust solution of (20) that is typically underdetermined

• • Data completion approach to handle missing data in square MIMO transfer function \(D(\omega )\)

• • Construction of passive ROMs

• • Convergence estimates

References

• [1]  Justin Baker and Elena Cherkaev and Vladimir Druskin and Shari Moskow and Mikhail Zaslavsky, Regularized Reduced Order Lippman-Schwinger-Lanczos Method for Inverse Scattering Problems in the Frequency Domain, submitted to Journal of Computational Physics, arXiv preprint arXiv:2311.16367, 2023

• [2]  Christopher A. Beattie, Volker Mehrmann, Paul Van Dooren, Robust port-Hamiltonian representations of passive systems, Automatica, 100, 2019: 182–186

• [3]  D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, Springer, Berlin, 2nd ed., 1998.

• [4]  Nakatsukasa, Yuji and Sète, Olivier and Trefethen, Lloyd N, The AAA Algorithm for Rational Approximation, SIAM Journal on Scientific Computing, 40 (3), 2018: A1494-A1522