# Integrated Sensing and Communication with Millimeter Wave Full Duplex Hybrid
Beamforming

Md Atiqul Islam2, George C. Alexandropoulos4, and Besma Smida2  2Department of
Electrical and Computer Engineering, University of Illinois at Chicago, USA
4Department of Informatics and Telecommunications, National and Kapodistrian
University of Athens, Greece
emails<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Integrated Sensing and Communication (ISAC) has attracted substantial
attraction in recent years for spectral efficiency improvement, enabling
hardware and spectrum sharing for simultaneous sensing and signaling
operations. In-band Full Duplex (FD) is being considered as a key enabling
technology for ISAC applications due to its simultaneous transmission and
reception capability. In this paper, we present an FD-based ISAC system
operating at millimeter Wave (mmWave) frequencies, where a massive Multiple-
Input Multiple-Output (MIMO) Base Station (BS) node employing hybrid Analog
and Digital (A/D) beamforming is communicating with a DownLink (DL) multi-
antenna user and the same waveform is utilized at the BS receiver for sensing
the radar targets in its coverage environment. We develop a sensing algorithm
that is capable of estimating Direction of Arrival (DoA), range, and relative
velocity of the radar targets. A joint optimization framework for designing
the A/D transmit and receive beamformers as well as the Self-Interference (SI)
cancellation is presented with the objective to maximize the achievable DL
rate and the accuracy of the radar target sensing performance. Our simulation
results, considering fifth Generation (5G) Orthogonal Frequency Division
Multiplexing (OFDM) waveforms, verify our approach’s high precision in
estimating DoA, range, and velocity of multiple radar targets, while
maximizing the DL communication rate.

###### Index Terms:

Full duplex, millimeter wave, direction estimation, range-Doppler estimation,
joint communication and sensing.

##  I Introduction

Integrated Sensing and Communication (ISAC) is an emerging concept for future
wireless networks, where the previously competing sensing and communication
operations are jointly optimized in the same hardware platform using a unified
signal processing framework [1, 2, 3, 4, 5]. Recently, Full Duplex (FD)
massive Multiple-Input Multiple-Output (MIMO) communications have been
considered a key enabler for ISAC applications due to their simultaneous
UpLink (UL) and DownLink (DL) transmission capability within the entire
frequency band [6, 7, 8, 9]. Furthermore, FD massive MIMO ISAC applications at
millimeter Wave (mmWave) frequencies have the potential to provide high
capacity communication links while simultaneously achieving high-resolution
sensing, e.g., Direction of Arrival (DoA), range, and relative speed of radar
targets/scatterers.

The performance of the FD ISAC systems relies on the in-band Self-Interference
(SI) signal suppression capability that stems from the Transmitter (TX) to the
Receiver (RX) side during FD operation. Recently in [8, 10, 11, 7], SI
cancellation was achieved for the FD massive MIMO systems operating at mmWave,
utilizing a combination of propagation domain isolation, analog domain
suppression, and digital SI cancellation techniques. To alleviate the hardware
cost in mmWave massive MIMO transceivers, hybrid Analog and Digital (A/D)
beamformers are usually employed, where large-scale antenna arrays are usually
connected to a small number of Radio Frequency (RF) chains via analog
preprocessing networks comprised of phase shifters [12]. Such systems require
appropriate beam selection for analog TX/RX beamformers, chosen from
predefined codebooks, to maximize DL rate and sensing accuracy. Moreover, in
the envisioned FD ISAC with massive MIMO radios, the transmit waveform will be
utilized for both DL data transmission and sensing of the radar targets [1].
Therefore, a joint design of the A/D beamformers and SI cancellation along
with sensing techniques is required for maximizing the performance of FD ISAC
systems.

Very recently in [13, 14], joint radar communication and sensing frameworks,
leveraging FD operation, were considered for single-antenna systems, where
both communication and radar waveforms were optimized for sensing performance.
In [15, 16], the FD ISAC operation was proposed for mmWave frequency bands
considering a massive MIMO FD Base Station (BS), where the signal power is
maximized in the radar target direction, while maintaining a threshold DL rate
performance. Although the considered FD ISAC approach estimated the DoA two
radar targets, the range was only calculated for one target due to its
disassociated DoA and range estimation technique.

Figure 1: The considered FD massive MIMO ISAC system: the FD hybrid A/D
beamforming Base Station (BS) communicates in the downlink with a mobile
single-antenna half-duplex User Equipment (UE), while its reflected
transmitted signals from radar targets/scatterers in the environment are
received via hybrid combining and processed for DoA, range, and relative
velocity estimations of the radar targets/scatterers.

In this paper, we present a novel FD massive MIMO ISAC system operating at
mmWave frequencies and realizing hybrid A/D beamforming, where Orthogonal
Frequency Division Multiplexing (OFDM) waveforms are utilized for both DL
communication and radar target sensing. Unlike state-of-the-art works, we
devise an ISAC optimization framework that is capable of estimating the DoA,
range, and relative velocity of multiple radar targets, while maximizing the
DL communication rate. Our numerical results, considering the fifth Generation
(5G) New Radio (NR) OFDM waveform, verify the high sensing accuracy and the
increased communication rate of the proposed optimization design.

##  II System and Signal Models

We consider an FD ISAC system comprising of an FD mmWave massive MIMO Base
Station (BS) node $b$ equipped with $N_{b}$ TX and $M_{b}$ RX antenna elements
communicating in the DL direction with an RX User Equipment (UE) node $u$ with
$M_{u}$ antenna elements, while the reflected DL signal is utilized to detect
radar targets/scatterers randomly distributed within the communication/sensing
environment at the RX of BS node $b$, as depicted in Fig. 1. To reduce the
number of RF chains and phase shifters, the BS node $b$ employs a partially-
connected Hybrid BeamForming (HBF) structure with A/D TX and RX beamformers,
where each of the $N_{b}^{\rm RF}$ and $M_{b}^{\rm RF}$ TX/RX RF chains is
connected to a Uniform Linear Arrays (ULA) of $N_{b}^{\rm A}$ and $M_{b}^{\rm
A}$ antenna elements, respectively, via phase shifters. Therefore, it holds
$N_{b}=N_{b}^{\rm A}N_{b}^{\rm RF}$ and $M_{b}=M_{b}^{\rm A}M_{b}^{\rm RF}$
for total number of TX and RX antennas, respectively, at the BS node $b$.
Since the size of the antenna array at the UE node $u$ is typically much
smaller than at the FD massive MIMO BS $b$, we assume that the UE adopts a
fully digital beamforming structure.

It is assumed that BS node $b$ transmits mmWave OFDM waveforms in the DL
direction containing $Q$ OFDM symbols with $P$ active subcarriers. In the
BaseBand (BB), the unit power frequency-domain symbol vector
$\mathbf{s}_{p,q}\in\mbox{$\mathbb{C}$}^{d_{b}\times 1}$ at the $p$th
subcarrier of $q$th OFDM symbol is precoded using digital beamforming matrix
$\mathbf{V}_{b}^{\mathrm{BB}}\in\mbox{$\mathbb{C}$}^{N_{b}^{\rm RF}\times
d_{b}}$, where $d_{b}\leq\min\\{N_{b}^{\mathrm{RF}},M_{u}\\}$. Following the
BB precoder, the DL signal is processed by the analog beamformer
$\mathbf{V}_{b}^{\mathrm{RF}}\in\mbox{$\mathbb{C}$}^{N_{b}\times N_{b}^{\rm
RF}}$ containing the configurations of the phase shifters as follows:

$\mathbf{V}_{b}^{\mathrm{RF}}\triangleq\left[\begin{matrix}\mathbf{v}_{1}&\mathbf{0}_{N_{b}^{{\rm
A}}\times 1}&\cdots&\mathbf{0}_{N_{b}^{{\rm A}}\times 1}\\\
\mathbf{0}_{N_{b}^{{\rm A}}\times
1}&\mathbf{v}_{2}&\cdots&\mathbf{0}_{N_{b}^{{\rm A}}\times 1}\\\
\vdots&\vdots&\ddots&\vdots\\\ \mathbf{0}_{N_{b}^{{\rm A}}\times
1}&\mathbf{0}_{N_{b}^{{\rm A}}\times 1}&\cdots&\mathbf{v}_{N_{b}^{({\rm
RF})}}\end{matrix}\right].$ (1)

The elements of each $\mathbf{v}_{i}$ are assumed to have constant magnitude,
i.e., $|[\mathbf{v}_{i}]_{n}|^{2}=1/N_{b}^{{\rm A}}$
$\forall$$n=1,2,\ldots,N_{b}^{{\rm A}}$. We also assume that
$\mathbf{v}_{i}\in\mathbb{F}_{\rm TX}$
$\forall$$i=1,2,\ldots,N_{b}^{\mathrm{RF}}$, which means that all analog TX
precoding vectors belong in a predefined beam codebook $\mathbb{F}_{\rm TX}$
including ${\rm card}(\mathbb{F}_{\rm TX})$ distinct vectors (or analog
beams). Applying both A/D beamforming, the TX frequency-domain symbols at the
antenna elements are expressed as

$\mathbf{x}_{p,q}\triangleq\mathbf{V}_{b}^{\mathrm{RF}}\mathbf{V}_{b}^{\mathrm{BB}}\mathbf{s}_{p,q},$
(2)

where $\mathbf{x}_{p,q}\in\mbox{$\mathbb{C}$}^{N_{b}\times 1}$. The signal
transmissions at node $b$ is power limited according to
$\mathbb{E}\\{\|\mathbf{V}_{b}^{\mathrm{RF}}\mathbf{V}_{b}^{\mathrm{BB}}\mathbf{s}_{p,q}\|^{2}\\}\leq{\rm
P}_{b}$.

###  II-A Radar Signal Reception Model

For sensing operation, we consider that a collection of $K$ radar
targets/scatterers are randomly distributed within the communication/sensing
environment and are to be detected by the BS node $b$. All the $K$ targets
reflect DL signal back to the RX of the BS node, while only a subset of them
($L$ out of $K$) contributes to the DL communication scattering paths between
the BS node $b$ and the RX UE node $u$. The purpose of the sensing operation
is to estimate the DoA, range, and relative velocity of each radar targets.

We consider that the DoAs of the $K$ targets are defined as
$\Theta=[\theta_{1},\theta_{2},\cdots,\theta_{K}]$, while the distance and the
relative speed of $k$th target correspond to a delay $\tau_{k}$ and a Doppler
shift $f_{D,k}$, respectively. Enabled by FD, the received signal
$\mathbf{y}_{p,q}\in\mathbb{C}^{M_{b}\times 1}$ at the RX of the BS node $b$
combining the SI and the reflected signal by the Radar targets is expressed as

$\begin{split}\mathbf{y}_{p,q}\triangleq&\sum\limits_{k=1}^{K}\alpha_{k}e^{j2\pi(qT_{s}f_{D,k}-p\tau_{k}\Delta
f)}\mathbf{a}_{M_{b}}(\theta_{k})\mathbf{a}_{N_{b}}^{\rm
H}(\theta_{k})\mathbf{x}_{p,q}\\\
&+\mathbf{H}_{b,b}\mathbf{x}_{p,q}+\mathbf{n}_{p,q},\end{split}$ (3)

where $\alpha_{k}\in\mbox{$\mathbb{C}$}$ and
$\mathbf{n}_{p,q}\sim\mathcal{CN}(0,\sigma_{b}^{2}\mathbf{I}_{M_{b}})$
represent the reflection coefficient of the $k$th radar target and the
receiver noise floor, respectively. Here, $\Delta f$ and $T_{s}\triangleq
1/\Delta f+T_{cp}$ denote the subcarrier spacing and the total OFDM symbol
duration (including the cyclic prefix). The propagation delay causes the phase
shift $e^{-j2\pi p\tau_{k}\Delta f}$ across subcarriers, while the Doppler
shift contributes a row-wise oscillation across different OFDM symbols [17].
Considering ULA, the steering vector $\mathbf{a}_{N_{b}}(\theta)$ for $N_{b}$
antenna elements and any DoA $\theta$ is formulated as [18]

$\begin{split}\mathbf{a}_{M_{b}}(\theta)\\!\triangleq\\!\frac{1}{\sqrt{N_{b}}}\big{[}1,e^{j\frac{2\pi}{\lambda}d\sin(\theta)},\ldots,e^{j\frac{2\pi}{\lambda}(N_{b}-1)d\sin(\theta)}]^{\rm
T},\end{split}$ (4)

where $\lambda$ is the propagation signal wavelength and $d$ denotes the
distance between adjacent antenna elements. Here,
$\mathbf{H}_{b,b}\in\mbox{$\mathbb{C}$}^{M_{b}\times N_{b}}$ is the SI channel
path at the BS node $b$, which is modeled as a Rician fading mmWave channel as
[19, eq. (9)].

The received signal at the node $b$ RX is first processed by the RF combiner
$\mathbf{W}_{b}^{\mathrm{RF}}\in\mbox{$\mathbb{C}$}^{M_{b}\times M_{b}^{\rm
RF}}$, where the structure of the combiner is formulated similarly as (1).
Here, the analog RX combining vectors belong in a predefined beam codebook
$\mathbb{F}_{\rm RX}$ including ${\rm card}(\mathbb{F}_{\rm RX})$ distinct
vectors. After RF combination and A/D SI cancellation, the received signal is
expressed as

$\displaystyle\widetilde{\mathbf{y}}_{p,q}\\!\triangleq$
$\displaystyle(\mathbf{W}_{b}^{\rm RF})^{\rm
H}\\!\sum\limits_{k=1}^{K}\\!\\!\alpha_{k}e^{j2\pi(qT_{s}f_{D,k}-p\tau_{k}\Delta
f)}\mathbf{a}_{M_{b}}\\!(\theta_{k})\mathbf{a}_{N_{b}}^{\rm
H}\\!(\theta_{k})\mathbf{x}_{p,q}$
$\displaystyle+(\widetilde{\mathbf{H}}_{b,b}\\!+\\!\mathbf{C}_{b}\\!+\\!\mathbf{D}_{b})\mathbf{V}_{b}^{\rm
BB}\mathbf{s}_{p,q}\\!+\\!(\mathbf{W}_{b}^{\rm RF})^{\rm H}\mathbf{n}_{p,q},$
(5)

where $\widetilde{\mathbf{H}}_{b,b}\triangleq(\mathbf{W}_{b}^{\rm RF})^{\rm
H}\mathbf{H}_{b,b}\mathbf{V}_{b}^{\rm RF}$ is the effective SI channel after
analog TX/RX beamforming. Here, $\mathbf{C}_{b}$ and $\mathbf{D}_{b}$
represent the analog and digital SI cancellation, respectively, that are
designed following the structure presented in [7].

###  II-B DL Signal Reception Model

As mentioned above, $L$ out of $K$ scatterers contributes to the DL channel
$\mathbf{H}_{u,b}\in\mbox{$\mathbb{C}$}^{M_{u}\times N_{b}}$. As the principal
focus of the paper is to estimate Radar target parameters, we ignore the delay
and Doppler shift parameters for the DL channel. Now, the received DL signal
vector $\mathbf{r}_{p,q}\in\mbox{$\mathbb{C}$}^{M_{u}\times d_{b}}$ at the UE
RX is expressed as

$\begin{split}\mathbf{r}_{p,q}&\triangleq\mathbf{W}_{u}^{\rm
H}\Big{(}\sum\limits_{\ell=1}^{L}\beta_{\ell}\mathbf{a}_{M_{u}}(\theta_{\ell})\mathbf{a}_{N_{b}}^{\rm
H}(\theta_{\ell})\mathbf{x}_{p,q}+\mathbf{z}_{p,q}\Big{)}\\\
&=\mathbf{W}_{u}^{\rm
H}\Big{(}\mathbf{H}_{u,b}\mathbf{x}_{p,q}+\mathbf{z}_{p,q}\Big{)},\end{split}$
(6)

where $\beta_{\ell}\in\mbox{$\mathbb{C}$}$ and
$\mathbf{z}_{p,q}\sim\mathcal{CN}(0,\sigma_{u}^{2}\mathbf{I}_{M_{u}})$
represent the complex reflection coefficient of $\ell$th scatter path and the
noise floor at RX node $u$, respectively.

The achievable DL rate of the FD ISAC system can be expressed as

$\displaystyle\mathcal{R}_{\rm DL}\triangleq\log_{2}\Big{(}{\rm det}\Big{(}$
$\displaystyle\mathbf{I}_{d_{b}}+\mathbf{W}_{u}^{\rm
H}\mathbf{H}_{u,b}\mathbf{V}_{b}^{\rm RF}\mathbf{V}_{b}^{\rm
BB}(\mathbf{V}_{b}^{\rm RF}\mathbf{V}_{b}^{\rm BB})^{\rm H}$
$\displaystyle\times\mathbf{H}_{u,b}^{\rm H}\mathbf{W}_{u}(\mathbf{W}_{u}^{\rm
H}\mathbf{W}_{u}\sigma_{u}^{2})^{-1}\Big{)}\Big{)},$ (7)

##  III DoA, Delay, and Doppler Shift Estimation

In this section, we present the estimations of the DoA, delay, and Doppler
shift parameters of the Radar targets, which are realized by the BS’s RX using
the reflected signals.

###  III-A Radar Target DoA Estimation

For DoA estimation, we deploy the MUltiple SIgnal Classification (MUSIC)
algorithm; other DoA estimation techniques can be used as well [20]. First, we
estimate the covariance matrix of the radar target reflected signal. Across
all subcarriers and OFDM symbols of the communication slot, the covariance
matrix $\mathbf{R}_{b}\in\mbox{$\mathbb{C}$}^{M_{b}^{\rm RF}\times M_{b}^{\rm
RF}}$ can be estimated as

$\begin{split}{\mathbf{R}}_{b}\\!\triangleq\\!\mathbb{E}\\{\widetilde{\mathbf{y}}_{p,q}\widetilde{\mathbf{y}}_{p,q}^{\rm
H}\\},\,\widehat{\mathbf{R}}_{b}\triangleq\frac{1}{PQ}\sum\limits_{q=0}^{Q-1}\sum\limits_{p=0}^{P-1}\widetilde{\mathbf{y}}_{p,q}\widetilde{\mathbf{y}}_{p,q}^{\rm
H}.\end{split}$ (8)

By taking the eigenvalue decomposition of the estimated sample covariance
matrix $\widehat{\mathbf{R}}_{b}$, it is deduced that:

$\begin{split}\widehat{\mathbf{R}}_{b}\triangleq\mathbf{U}{\rm
diag}\\{\eta_{1},\eta_{2},\ldots,\eta_{M_{b}}\\}\mathbf{U}^{\rm
H},\end{split}$ (9)

where $\eta_{1}\geq\eta_{2}\geq\ldots\geq\eta_{M_{b}^{\rm RF}}$ are the
eigenvalues of $\widehat{\mathbf{R}}_{b}$ and
$\mathbf{U}\in\mathbb{C}^{M_{b}^{\rm RF}\times M_{b}^{\rm RF}}$ contains their
corresponding eigenvectors. Since we are interested in estimating the DoAs of
$K$ radar targets, the matrix $\mathbf{U}$ can be partitioned as
$\mathbf{U}=[\mathbf{U}_{s}|\mathbf{U}_{n}]$, where the columns in
$\mathbf{U}_{n}\in\mathbb{C}^{M_{b}^{\rm RF}\times M_{b}^{\rm RF}-K}$ are the
eigenvectors spanning the noise subspace and
$\mathbf{U}_{s}\in\mathbb{C}^{M_{b}^{\rm RF}\times K}$ contains the signal
space eigenvectors. The MUSIC spectrum for the considered HBF architecture can
be thus formulated as:

$\begin{split}S(\theta)\triangleq\left(\mathbf{a}_{M_{b}}^{\rm
H}(\theta)\mathbf{W}_{b}^{\rm RF}\mathbf{U}_{n}\mathbf{U}_{n}^{\rm
H}(\mathbf{W}_{b}^{\rm RF})^{\rm
H}\mathbf{a}_{M_{b}}(\theta)\right)^{-1},\end{split}$ (10)

whose $K$ peaks correspond to the $K$ estimated DoAs
$\widehat{\theta_{k}},\forall k$.

Algorithm 1 Delay and Doppler Shift Estimation

0: $\mathbf{x}_{p,q}$, $\widetilde{\mathbf{y}}_{p,q}$,$\forall p,q$,
$\mathbf{W}_{b}^{\rm RF}$, and $\widehat{\theta}_{k}\,\,\forall k$.

0: $\widehat{\tau}_{k},\widehat{f}_{D,k}$, $\forall k$.

1: Set $n=0,\cdots,P-1$ and $m=-\frac{Q}{2},\cdots,\frac{Q}{2}-1$.

2: for $k=1,2,\ldots,K$ do

3: Set
$\mathbf{g}_{p,q}\triangleq\mathbf{a}_{M_{b}}(\widehat{\theta}_{k})\mathbf{a}_{N_{b}}^{\rm
H}(\widehat{\theta}_{k})\mathbf{x}_{p,q}$.

4: Set
$z_{p,q}\triangleq\frac{1}{M_{b}}\sum\limits_{i=1}^{M_{b}}[\mathbf{W}_{b}^{\rm
RF}\widetilde{\mathbf{y}}_{p,q}]_{i}/[\mathbf{g}_{p,q}]_{i},$ $\forall p,q$.

5: Set
$A(n,m)\triangleq\sum\limits_{p=0}^{P-1}\left(\sum\limits_{q=0}^{Q-1}z_{p,q}e^{-j2\pi\frac{qm}{Q}}\right)e^{j2\pi\frac{pn}{P}}$.

6: Set $({n}^{*},{m}^{*})=\underset{n,m}{\text{arg max}}\quad|A(n,m)|^{2}$.

7: Set the estimated delay $\widehat{\tau}_{k}=\frac{{n}^{*}}{P\Delta f}$.

8: Set the Doppler frequency of $k$th target
$\widehat{f}_{D,k}=\frac{{m}^{*}}{QT_{s}}$.

9: end for

###  III-B Delay and Doppler Shift Estimation

The next step is to estimate the delay and Doppler shift parameters associated
with the $K$ estimated DoAs. Using the estimate DoA $\widehat{\theta}_{k}$ and
the known transmit signal $\mathbf{x}_{p,q}$, we formulate a reference signal
in the DoA direction as

$\begin{split}\mathbf{g}_{p,q}\triangleq\mathbf{a}_{M_{b}}(\widehat{\theta}_{k})\mathbf{a}_{N_{b}}^{\rm
H}(\widehat{\theta}_{k})\mathbf{x}_{p,q}.\end{split}$ (11)

Now, we utilize the received signal $\widetilde{\mathbf{y}}_{p,q}$ to derive
the quotient averaged across all RX antennas that includes the effect of delay
and Doppler shift in the direction of $\widehat{\theta}_{k}$ as

$\begin{split}z_{p,q}\triangleq\frac{1}{M_{b}}\sum\limits_{i=1}^{M_{b}}[\mathbf{W}_{b}^{\rm
RF}\widetilde{\mathbf{y}}_{p,q}]_{i}/[\mathbf{g}_{p,q}]_{i},\,\,\forall
p,q.\end{split}$ (12)

To estimate $\widehat{\tau}_{k}$ and $\widehat{f}_{D,k}$, we formulate the
likelihood function:

$\begin{split}A(n,m)\triangleq\sum\limits_{p=0}^{P-1}\left(\sum\limits_{q=0}^{Q-1}z_{p,q}e^{-j2\pi\frac{qm}{Q}}\right)e^{j2\pi\frac{pn}{P}},\end{split}$
(13)

where $n=0,\cdots,P-1$ and $m=-Q/2,\cdots,Q/2-1$. Finally, we find the best
quantized delay and Doppler shift that maximizes the likelihood function norm.
The delay and Doppler shift estimation procedure is provided in Algorithm 1.

##  IV Proposed ISAC Optimization Framework

In this section, we present a joint optimization framework deriving A/D
beamformers and SI cancellation matrices with the objective to maximize the DL
rate and the radar estimation accuracy.

We consider a time division duplexing communication protocol, where the DoAs
estimated in one communication time slot is utilized to derive the beamformers
and SI cancellation matrices for the successive slot. To optimize the DL rate
and radar target parameter estimation accuracy, we propose to maximize the SNR
in both the DL and radar target direction. Given the estimated DoAs
$\widehat{\theta}_{k},\,\forall k$, the SNR in the direction of all radar
targets can be written as

$\begin{split}\widehat{{\Gamma}}_{\rm
Radar}\triangleq\Big{\|}(\mathbf{W}_{b}^{\rm RF})^{\rm
H}\sum\limits_{k=1}^{K}\mathbf{a}_{M_{b}}(\widehat{\theta}_{k})\mathbf{a}_{N_{b}}^{\rm
H}(\widehat{\theta}_{k})\mathbf{V}_{b}^{\rm RF}\mathbf{V}_{b}^{\rm
BB}\Big{\|}^{2}{\Sigma}_{b}^{-1},\end{split}$ (14)

where
${\Sigma}_{b}=\|(\widehat{\widetilde{\mathbf{H}}}_{b,b}+\mathbf{C}_{b}+\mathbf{D}_{b})\mathbf{V}_{b}^{\rm
BB}\|^{2}+\left\|\mathbf{W}_{b}^{\rm RF}\right\|^{2}\sigma_{b}^{2}$ is the
interference plus noise covariance matrix at the RX of node $b$. Similarly,
the estimated DL SNR is expressed as

$\begin{split}\widehat{{\Gamma}}_{\rm DL}=\Big{\|}\mathbf{W}_{u}^{\rm
H}\sum\limits_{\ell=1}^{L}\mathbf{a}_{M_{u}}(\widehat{\theta}_{\ell})\mathbf{a}_{N_{b}}^{\rm
H}(\widehat{\theta}_{\ell})\mathbf{V}_{b}^{\rm RF}\mathbf{V}_{b}^{\rm
BB}\Big{\|}^{2}{\Sigma}_{u}^{-1},\end{split}$ (15)

where ${\Sigma}_{u}=\left\|\mathbf{W}_{u}\right\|^{2}\sigma_{u}^{2}$ is the
noise covariance matrix at the RX node $u$.

The optimization problem to maximize the Radar target and DL SNR can be
written as

$\displaystyle\mathcal{OP}$
$\displaystyle:\underset{\begin{subarray}{c}\mathbf{V}_{b}^{\rm
RF},\mathbf{V}_{b}^{\rm BB},\mathbf{W}_{b}^{\rm RF}\\\
\mathbf{C}_{b},\mathbf{D}_{b},\mathbf{W}_{u}\end{subarray}}{\max}\quad\widehat{{\Gamma}}_{\rm
Radar}+\widehat{{\Gamma}}_{\rm DL}$
$\displaystyle\text{\text{s}.\text{t}.}\quad\left\|[(\widehat{\widetilde{\mathbf{H}}}_{b,b}\\!+\\!\mathbf{C}_{b})\mathbf{V}_{b}^{\rm
BB}]_{(j,:)}\right\|^{2}\\!\\!\\!\leq\\!\lambda_{b},\\!\forall\\!j\\!=1,\ldots,M_{b}^{\rm
RF},$ $\displaystyle\quad\quad\mathbb{E}\\{\|\mathbf{V}_{b}^{\rm
RF}\mathbf{V}_{b}^{\rm BB}\|^{2}\\}\leq{\rm P}_{b},$ (16)
$\displaystyle\quad\quad\mathbf{w}_{j}\in\mathbb{F}_{\rm RX}\,\,\forall
j\,\,{\rm and}\,\,\mathbf{v}_{n}\in\mathbb{F}_{\rm TX}\,\,\forall
n=1,2,\ldots,N_{b}^{({\rm RF})}$

Algorithm 2 FD ISAC Optimization

0: $\widehat{\mathbf{H}}_{b,b}$, $N$, ${\rm P}_{b}$ and
$\widehat{\theta}_{k}\,\,\forall k$.

0: $\mathbf{V}_{b}^{\rm RF},\mathbf{V}_{b}^{\rm BB},\mathbf{W}_{b}^{\rm
RF},\mathbf{C}_{b},\mathbf{D}_{b}$, and $\mathbf{W}_{u}$.

1: Set $\widehat{\mathbf{H}}_{\rm
R}\triangleq\sum\limits_{k=1}^{K}\mathbf{a}_{M_{b}}(\widehat{\theta}_{k})\mathbf{a}_{N_{b}}^{\rm
H}(\widehat{\theta}_{k})$.

2: Set $\widehat{\mathbf{H}}_{\rm
u,b}\triangleq\sum\limits_{\ell=1}^{L}\mathbf{a}_{M_{u}}(\widehat{\theta}_{\ell})\mathbf{a}_{N_{b}}^{\rm
H}(\widehat{\theta}_{\ell})$.

3: Set $\mathbf{W}_{u}$ as the $d_{b}$ left-singular vectors of
$\widehat{\mathbf{H}}_{\rm u,b}$ corresponding to the singular values in
descending order.

4: Set $\mathbf{V}_{b}^{\rm
RF}\triangleq\underset{\mathbf{v}_{j}\in\mathbb{F}_{\rm TX}}{\text{arg
max}}\quad\|\widehat{\mathbf{H}}_{\rm R}\mathbf{V}_{b}^{\rm RF}\|^{2}$.

5: Set $\mathbf{W}_{b}^{\rm
RF}\triangleq\underset{\mathbf{w}_{j}\in\mathbb{F}_{\rm TX}}{\text{arg
max}}\quad\frac{\|(\mathbf{W}_{b}^{\rm RF})^{\rm H}\widehat{\mathbf{H}}_{\rm
R}\mathbf{V}_{b}^{\rm RF}\|^{2}}{\|(\mathbf{W}_{b}^{\rm RF})^{\rm
H}\widehat{\mathbf{H}}_{\rm b,b}\mathbf{V}_{b}^{\rm RF}\|^{2}}$.

6: Set $\widehat{\widetilde{\mathbf{H}}}_{b,b}=(\mathbf{W}_{b}^{\rm RF})^{\rm
H}\widehat{\mathbf{H}}_{b,b}\mathbf{V}_{b}^{\rm RF}$ and
$\widehat{\widetilde{\mathbf{H}}}_{u,b}=\widehat{\mathbf{H}}_{u,b}\mathbf{V}_{b}^{\rm
RF}$.

7: Set
$\mathbf{C}_{b}=\begin{bmatrix}-[\widehat{\widetilde{\mathbf{H}}}_{b,b}]_{:,1:\frac{N}{{M}_{b}^{\rm
RF}}}&\mathbf{0}_{:,({N}_{b}^{\rm RF}-\frac{N}{{M}_{b}^{\rm RF}}:{N}_{b}^{\rm
RF})}\end{bmatrix}$ and
$\mathbf{D}_{b}=-(\widehat{\widetilde{\mathbf{H}}}_{b,b}+C_{b})$.

8: Obtain $\mathbf{B}$ with the $N_{b}^{\rm RF}$ right-singular vectors of
$(\widehat{\widetilde{\mathbf{H}}}_{b,b}+\mathbf{C}_{b})$ corresponding to the
singular values in descending order.

9: for $\alpha={N}_{b}^{\rm RF},N_{b}^{\rm RF}-1,\ldots,2$ do

10: Set $\mathbf{F}=[\mathbf{B}]_{(:,N_{b}^{\rm RF}-\alpha+1:N_{b}^{\rm
RF})}$.

11: Set $\mathbf{G}$ as the optimum precoding for the effective DL MIMO
channel $\widehat{\widetilde{\mathbf{H}}}_{u,b}\mathbf{F}$ given ${\rm
P}_{b}$.

12: if
$\left\|[(\widehat{\widetilde{\mathbf{H}}}_{b,b}\\!+\\!\mathbf{C}_{b})\mathbf{F}\mathbf{G}]_{(j,:)}\right\|^{2}\\!\\!\\!\leq\\!\lambda_{b},\\!\forall\\!j\\!=1,\ldots,M_{b}^{\rm
RF}$, then

13: Output $\mathbf{V}_{b}^{\mathrm{BB}}=\mathbf{F}\mathbf{G}$ and stop the
algorithm.

14: else

15: Output that the $\mathbf{C}_{b}$ realization does not meet the residual SI
constraint.

16: end if

17: end for

The optimization problem in (IV) is a non-convex problem with coupling
variables, hence, quite difficult to tackle. In this work, we solve it
suboptimally using alternating optimization, leaving other possibilities for
future work.

First, using the estimated DoAs, we formulate a virtual channel of in the
radar target direction as $\widehat{\mathbf{H}}_{\rm
R}\triangleq\sum\limits_{k=1}^{K}\mathbf{a}_{M_{b}}(\widehat{\theta}_{k})\mathbf{a}_{N_{b}}^{\rm
H}(\widehat{\theta}_{k})$. Now to maximize the radar SNR, we find the TX
analog beams solving the following suboptimization problem:

$\displaystyle\mathcal{OP}1:$ $\displaystyle\quad\mathbf{V}_{b}^{\rm
RF}\triangleq\underset{\mathbf{v}_{j}\in\mathbb{F}_{\rm
TX}}{\operatornamewithlimits{argmax}}\quad\|\widehat{\mathbf{H}}_{\rm
R}\mathbf{V}_{b}^{\rm RF}\|^{2}$ (17)

The solution of $\mathcal{OP}1$ requires a simple search through the beam
codebook $\mathbb{F}_{\rm TX}$. Using the TX analog beamformer
$\mathbf{V}_{b}^{\rm RF}$, we derive the analog combiner solving the
suboptimization problem $\mathcal{OP}1$, where we simultaneously maximize the
radar SNR and suppress SI signal as follows:

$\displaystyle\mathcal{OP}2:$ $\displaystyle\quad\mathbf{W}_{b}^{\rm
RF}\triangleq\underset{\mathbf{w}_{j}\in\mathbb{F}_{\rm
TX}}{\operatornamewithlimits{argmax}}\quad\frac{\|(\mathbf{W}_{b}^{\rm
RF})^{\rm H}\widehat{\mathbf{H}}_{\rm R}\mathbf{V}_{b}^{\rm
RF}\|^{2}}{\|(\mathbf{W}_{b}^{\rm RF})^{\rm H}\widehat{\mathbf{H}}_{\rm
b,b}\mathbf{V}_{b}^{\rm RF}\|^{2}}$ (18)

Similar to $\mathcal{OP}1$, the solution of $\mathcal{OP}2$ requires a simple
search through the available beam codebook $\mathbb{F}_{\rm RX}$. Given the
analog SI cancellation taps $N\leq{N}_{b}^{\rm RF}{M}_{b}^{\rm RF}$ and the
analog TX/RX beamformer at the BS, we follow a similar procedure to [7],
aiming to find the digital beamforming matrix $\mathbf{V}_{b}^{\rm BB}$ and
the SI cancellation matrices maximizing the DL rate and suppressing the SI
signal power below the RF saturation level of $\lambda_{b}$. The latter will
ensure proper reception of the Radar target reflected signal. Our solution for
the optimization problem (IV) is summarized in Algorithm 2.

Figure 2: Sensing parameter estimation performance for $6$ radar targets with
$N_{b}^{\rm RF}=M_{b}^{\rm RF}=8,N_{b}^{\rm A}=N_{b}^{\rm A}=16$ and transmit
power of $30$dBm. (a) DoA and range estimation, (b) Range and relative
velocity estimation.

##  V Numerical Results

In this section, we present numerical results for the radar sensing and DL
rate performance of the proposed FD massive MIMO ISAC system operating at
mmWave frequencies.

###  V-A Simulation Parameters

We perform an extensive waveform simulation following the FD massive MIMO
architecture illustrated in Fig. 1 when operating at mmWave frequencies, where
a $128\times 128$ FD massive MIMO node $b$ is communicating in the DL
direction with $4$ antenna RX UE node $u$. The BS node $b$ employs $N_{b}^{\rm
RF}=M_{b}^{\rm RF}=8$ TX/RX RF chains with each of them connected to a ULA of
$N_{b}^{\rm A}=M_{b}^{\rm A}=16$ antenna elements via phase shifters. The
communication is performed using mmWave frequency of $28$GHz and a 5G NR OFDM
waveform with $100$MHz BandWidth (BW) and $\Delta f=120$KHz subcarrier
spacing. According to the 5G NR specifications, we have $66$ Physical Resource
Blocks (PRBs) resulting in $P=792$ active subcarriers and $Q=14$ OFDM symbols
in each communication slot. Total OFDM symbol duration is defined as
$T_{s}=8.92\mu$s. We have considered a radio subframe of $1$ms for DL
communication. The RX noise floors at all nodes were assumed to be $-90$dBm
for $100$MHz BW OFDM signal. To this end, the RXs have an effective dynamic
range of $60$dB provided by $14$-bit Analog-to-Digital Converters (ADC) for a
Peak-to-Average-Power-Ratio (PAPR) of $10$ dB. Therefore, the residual SI
power after analog SI cancellation at the input of each RX RF chain has to be
below $-30$dBm to avoid signal saturation. The pathloss of the clustered DL
channel is assumed to be $100$dB, whereas the SI channels are modeled as
Rician fading channels with a $\kappa$-factor of $35$dB and pathloss $40$dB
[21]. For the BS analog TX/RX beamformer, we consider a $5$-bit beam codebook
based on the Discrete Fourier Transform (DFT) matrix. We have used $1000$
independent Monte Carlo simulation runs to calculate the Radar sensing and DL
rate performance.

###  V-B Radar Target Parameters

We have considered $K=6$ radar targets randomly distributed in the
sensing/communication environment with DoAs $\theta_{k}\in[-90^{\circ}\quad
90^{\circ}],\forall k$. For communication scatters, $L=2$ out of $K=6$ targets
are chosen randomly. Each of the radar targets the range and relative velocity
is selected randomly with a maximum range of $80$m and maximum velocity of
$100$km/h.

Figure 3: Sensing parameter estimation performance for $6$ radar targets with
$N_{b}^{\rm RF}=M_{b}^{\rm RF}=8,N_{b}^{\rm A}=N_{b}^{\rm A}=16$ and DL
transmit power of $10$dBm. (a) DoA and range estimation, (b) Range and
relative velocity estimation.

###  V-C Radar Target Sensing performance

In Fig. 2, We have depicted the sensing performance of the proposed FD ISAC
system with a $128\times 128$ massive MIMO node transmitting DL signal with a
transmit power of $30$dBm. The DoA and range estimation is plotted in contrast
to the true target parameters in Fig. 2(a), where it is evident that the
proposed FD ISAC system can detect all $6$ targets successfully with high
precision. Even for really close target ($10^{\circ}$ and $12^{\circ}$) with
only $2^{\circ}$ angle and less than $5$m range difference, the estimation
performance is highly accurate. The superior sensing performance is provided
by the proposed associated delay estimation approach with high-resolution
MUSIC DoA estimation unlike previous FD ISAC work in [16], where range
estimation for such close targets was not possible. In Fig. 2(b), the relative
velocity is plotted with respect to the range estimation for $30$dBm DL
transmit power. The figure shows that the proposed FD ISAC system is capable
of estimating the relative velocity of all $6$ targets with less than $1.5\%$
estimation error.

In Fig. 3, radar target sensing performance of the proposed FD ISAC system is
presented for $10$dBm DL transmit power. It is evident from Fig. 3(a) that the
DoA estimation is almost accurate even at a low transmit power of $10$dBm.
However, for the target at $3^{\circ}$, the estimated range is around $3$m
away the actual value of $80$m. This is due to the low transmit power and
higher path loss of the furthest target. In Fig. 3(b), the relative velocity
estimation is showcased, where the proposed approach achieved sensing
performance with less than $5\%$ estimation error for a low transmit power of
$10$dBm.

Figure 4: DL rate with respect to transmit power in dBm for the $128\times
128$ massive MIMO FD BS communicating with a $4$ antenna UE RX node.

The DL rate performance of the proposed FD ISAC system with a $128\times 128$
massive MIMO BS node transmitting to a $4$ antenna UE RX node is depicted in
Fig. 4 with respect to transmit power. It is evident from the figure that the
proposed FD ISAC approach is capable of providing DL rate very close
($1.5$bps/Hz) to the ideal rate in addition to the high radar sensing
performance at the RX of BS node for $N=16$ ($25\%$) analog SI cancellation
taps at the FD massive MIMO node. For even smaller hardware complexity
($12.5\%$) taps, the DL rate performance is comparable up to $20$ dB transmit
power. However, as transmit power increases the impact of SI signal worsens,
for DL transmit power of $30$dBm, the FD ISAC system is capable of providing
$75\%$ of the ideal DL rate with only $N=8$ SI cancellation taps.

##  VI Conclusion

In this paper, we presented an FD-based ISAC optimization framework, where an
FD massive MIMO BS node is transmitting DL signals and concurrently performing
radar target sensing utilizing the reflected signals. We devised a DoA, delay,
and Doppler shift estimation algorithm for multiple radar target sensing
considering hybrid A/D beamforming at the BS node. Adopting a limited
complexity analog SI cancellation architecture, we presented a joint design of
the A/D beamformer and SI cancellation that maximizes the DL rate together
with the target sensing performance. Our performance results for a mmWave
channel model demonstrated the high precision DoA, range, velocity estimation
of multiple radar targets while providing maximized DL rate. In future work,
the proposed FD ISAC approach will be considered for sensing and communication
demanding practical applications, such as autonomous vehicles and flight
control systems.

## Acknowledgments

This work was partially funded by the National Science Foundation CAREER award
#1620902.

## References

  * [1] F. Liu, C. Masouros, A. P. Petropulu, H. Griffiths, and L. Hanzo, “Joint radar and communication design: Applications, state-of-the-art, and the road ahead,” _IEEE Trans. Commun._ , vol. 68, no. 6, pp. 3834–3862, Feb. 2020\.
  * [2] I. Alamzadeh, G. C. Alexandropoulos, N. Shlezinger, and M. F. Imani, “A reconfigurable intelligent surface with integrated sensing capability,” _Scientific Reports_ , vol. 11, no. 20737, pp. 1–10, Oct. 2021.
  * [3] G. C. Alexandropoulos, N. Shlezinger, I. Alamzadeh, M. F. Imani, H. Zhang, and Y. C. Eldar, “Hybrid reconfigurable intelligent metasurfaces: Enabling simultaneous tunable reflections and sensing for 6G wireless communications,” _[online] https://arxiv.org/abs/2104.04690_ , 2021.
  * [4] K. V. Mishra, M. B. Shankar, V. Koivunen, B. Ottersten, and S. A. Vorobyov, “Toward millimeter-wave joint radar communications: A signal processing perspective,” _IEEE Signal Process. Mag._ , vol. 36, no. 5, pp. 100–114, Sep. 2019.
  * [5] B. Paul, A. R. Chiriyath, and D. W. Bliss, “Survey of RF communications and sensing convergence research,” _IEEE Access_ , vol. 5, pp. 252–270, Jan. 2016.
  * [6] H. Alves, T. Riihonen, and H. A. Suraweera, _Full-Duplex Communications for Future Wireless Networks_. Singapore: Springer, 2020.
  * [7] G. C. Alexandropoulos, M. A. Islam, and B. Smida, “Full duplex hybrid A/D beamforming with reduced complexity multi-tap analog cancellation,” in _Proc. IEEE SPAWC_ , Atlanta, USA, May 2020.
  * [8] M. A. Islam, G. C. Alexandropoulos, and B. Smida, “Direction-assisted beam management in full duplex millimeter wave massive MIMO systems,” in _Proc. IEEE GLOBECOM_ , Madrid, Spain, Dec. 2021.
  * [9] ——, “Simultaneous data communication and channel estimation in multi-user full duplex MIMO systems,” in _Proc. IEEE ASILOMAR_ , Pacific Grove, USA, Nov. 2020.
  * [10] Z. Xiao, P. Xia, and X.-G. Xia, “Full-duplex millimeter-wave communication,” _IEEE Wireless Commun._ , vol. 24, no. 6, pp. 136–143, Dec. 2017.
  * [11] I. P. Roberts, H. B. Jain, and S. Vishwanath, “Equipping millimeter-wave full-duplex with analog self-interference cancellation,” in _Proc. IEEE ICC_ , Dublin, Ireland, Jun. 2020.
  * [12] V. Venkateswaran and A.-J. van der Veen, “Analog beamforming in MIMO communications with phase shift networks and online channel estimation,” _IEEE Trans. Signal Process._ , vol. 58, no. 8, pp. 4131–4143, Aug. 2010\.
  * [13] C. B. Barneto, T. Riihonen, M. Turunen, L. Anttila, M. Fleischer, K. Stadius, J. Ryynänen, and M. Valkama, “Full-duplex OFDM radar with LTE and 5G NR waveforms: Challenges, solutions, and measurements,” _IEEE Trans. Microw. Theory Techn._ , vol. 67, no. 10, pp. 4042–4054, Aug. 2019.
  * [14] S. D. Liyanaarachchi, T. Riihonen, C. B. Barneto, and M. Valkama, “Optimized waveforms for 5G–6G communication with sensing: Theory, simulations and experiments,” _IEEE Trans. Wireless Commun._ , Jun. 2021.
  * [15] C. B. Barneto, S. D. Liyanaarachchi, T. Riihonen, M. Heino, L. Anttila, and M. Valkama, “Beamforming and waveform optimization for OFDM-based joint communications and sensing at mm-waves,” in _Proc. IEEE ASILOMAR_ , Pacific Grove, USA, Nov. 2020, pp. 895–899.
  * [16] S. D. Liyanaarachchi, C. B. Barneto, T. Riihonen, M. Heino, and M. Valkama, “Joint multi-user communication and MIMO radar through full-duplex hybrid beamforming,” in _IEEE Int. Symp. Joint Commun. & Sensing (JC&S)_, Feb. 2021, pp. 1–5.
  * [17] M. Braun, C. Sturm, and F. K. Jondral, “Maximum likelihood speed and distance estimation for OFDM radar,” in _IEEE Radar Conf._ , May 2010, pp. 256–261.
  * [18] J. P. González-Coma _et al._ , “Channel estimation and hybrid precoding for frequency selective multiuser mmWave MIMO systems,” _IEEE J. Sel. Topics Signal Proces._ , vol. 12, no. 2, pp. 353–367, Mar. 2018\.
  * [19] K. Satyanarayana _et al._ , “Hybrid beamforming design for full-duplex millimeter wave communication,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 2, pp. 1394–1404, Dec. 2018.
  * [20] H. Krim and M. Viberg, “Two decades of array signal processing research: the parametric approach,” _IEEE Signal Process. Mag._ , vol. 13, no. 4, pp. 67–94, Jul. 1996.
  * [21] G. C. Alexandropoulos and M. Duarte, “Joint design of multi-tap analog cancellation and digital beamforming for reduced complexity full duplex MIMO systems,” in _Proc. IEEE ICC_ , Paris, France, May 2017.


#  Lippmann-Schwinger-Lanczos algorithm for inverse scattering problems

V. Druskin111Worcester Polytechnic Institute, Department of Mathematical
Sciences, Stratton Hall, 100 Institute Road, Worcester MA, 01609
<EMAIL_ADDRESS>S. Moskow222Department of Mathematics, Drexel
University, Korman Center, 3141 Chestnut Street, Philadelphia, PA 19104
<EMAIL_ADDRESS>M. Zaslavsky333Schlumberger-Doll Research Center, 1
Hampshire St., Cambridge, MA 02139-1578<EMAIL_ADDRESS>
###### Abstract

Data-driven reduced order models (ROMs) are combined with the Lippmann-
Schwinger integral equation to produce a direct nonlinear inversion method.
The ROM is viewed as a Galerkin projection and is sparse due to Lanczos
orthogonalization. Embedding into the continuous problem, a data-driven
internal solution is produced. This internal solution is then used in the
Lippmann-Schwinger equation, thus making further iterative updates
unnecessary. We show numerical experiments for spectral domain domain data for
which our inversion is far superior to the Born inversion and works as well as
when the true internal solution is known.

##  1 Introduction

This work extends the reduced order model (ROM) approach to inverse impedance,
scattering and diffusion problems that was developed in the sequel of papers
[5, 7, 14, 11, 8, 9, 6, 2]. This approach to solve multidimensional inverse
problems using a ROM framework can be summarized by the following:

  1. 1.

The boundary data, (discrete partial Dirichlet-to-Neumann maps as in [5] or
their transient variant [14, 11, 8, 9, 6]) is matched by data-driven ROMs. The
ROMs are sparse networks that can be written as a tridiagonal or block-
tridiagonal matrix for one-dimensional and multi-dimensional problems
respectively. Matching is performed via a direct layer stripping algorithm or
a sequence of such algorithms.

  2. 2.

The data-driven networks implicitly embed boundary data back into the
interior. The main property of such embeddings is that the network’s
coefficients can be approximated by localized averages of the PDE
coefficients. Approximate linear maps between these averages and the network
coefficients can be learned via so-called ”optimal grid” or ”‘finite-
difference Gaussian quadrature” [5, 14], or this map can be used as a
preconditioner in optimization algorithms [7, 2].

In short, the main nonlinearity of the inverse problem is absorbed during the
first stage thanks to the layer stripping, which is intrinsically nonlinear.
This allows one to treat the transformed data as almost linear with respect to
PDE coefficients. The layer stripping in essence is just the Euclidean
polynomial division algorithm, which was further developed in the seminal
works of Stieltjes and Lanczos, and later in the electrical circuit community.
It also related to seminal works from the Soviet school (Marchenko, Gelfand,
Levitan and Krein) on inverse spectral problems. The origin of the embedding
concept can be be traced to Krein [15], and to the idea of the spectrally
matched second order staggered finite-difference grids first introduced in
[13]. These ”optimal grids” or ”finite-difference Gaussian quadrature” rules
provided a clear geometric interpretation of the data driven network’s
coefficients [4].

It is well known that data-driven ROM can be equivalently rewritten in
projection form, for example as a Galerkin system, with the same matrix $T$
for state variable realization in the interior, see for example [1]. The
relationship between these projections and finite-difference Gaussian
quadrature rules was first studied in [12] and further developed in [6]. To
see the idea, let $A_{q}$ be the PDE operator, for example

$A_{q}=-\Delta+q,$ (1)

and let us denote by $V_{q}$ the row vector of orthonormal basis functions
${v_{q}}_{j}$, $j=1,\ldots,m$, such that

$T_{q}=\int V^{*}_{q}A_{q}V_{q}$ (2)

where $T_{q}$ is our sparse data driven network. A critical property of the
projection solution in the interior, which was first noticed in [14], is that

$V_{q}\approx V_{0},$ (3)

even for large $q$, where $V_{0}$ is the row vector of orthogonalized basis
corresponding to $q=0$. This property of the sparse realization $T_{q}$ would
not be possible, for example, if one were using the full stiffness and mass
matrices appearing in the Loewner product formulation normally used for data-
driven ROM construction. Because of the sparse structure of $T_{q}$,
(tridiagonal finite-difference in 1D problems), the basis functions
${v_{q}}_{j}$ are localized, that is, they behave somehow similarly to
piecewise linear finite element basis functions, and depend only weakly on the
media perturbations. Rigorous analysis of this property can be obtained using
already mentioned Marchenko, Gelfand, Levitan and Krein approach and is in
progress at the moment [3].

Numerical experiments have demonstrated that the weak dependence of the basis
functions on $q$ also holds in the multidimensional setting, and this gave
rise to so-called nonlinear back-projection imaging algorithm [11], that uses
the approximation

$A_{q}w\approx V_{0}T_{q}\int V_{0}^{*}w.$ (4)

The backprojection algorithm worked very well for large $m$ yielding a high
resolution basis for time domain wave problems, but failed in diffusion
problems, allowing only small $m$ due to inherent ill-posedness. To overcome
this problem, in [6] the authors with collaborators introduced the idea of
generating internal solutions. To see the idea of this method, we first
consider the one dimensional problem with source term $g$. Then the frequency
domain solution can be written as

$u_{q}(x,s)=(A_{q}+sI)^{-1}g.$ (5)

Then using (4) we obtain

$u_{q}(x,s)\approx\tilde{u}_{q}(x,s)=V_{0}(T_{q}+sI)^{-1}\int V_{0}^{*}g.$ (6)

Thanks to (4), $\tilde{u}_{q}-u_{q}<<u_{q}-u_{0}$, even when $q$ is rather
large. Then the approximate solution of inverse problem can be computed as

$q\approx\frac{\Delta\tilde{u}_{q}(x,s)-s\tilde{u}_{q}(x,s)+g(x)}{\tilde{u}_{q}(x,s)}.$
(7)

This formulation worked better for small $m$ where back-projection failed,
however, image quality was not high for large $q$. The reason was that
multiplication by $\Delta$ in (7) amplified the approximation error of the
internal solution $u_{q}(x,s)-\tilde{u}_{q}(x,s)$.

To overcome this artifact here we consider the Lippmann-Schwinger integral
equation framework. Consider data

$F_{q}(s)=\int gu_{q}=\int g(A_{q}+sI)^{-1}g,$

and similarly

$F_{0}(s)=\int gu_{0}=\int g(A_{0}+sI)^{-1}g.$

Then the Lippmann-Schwinger integral equation with respect to the unknown $q$
can be written as

$F_{q}(s)-F_{0}(s)=-\langle u_{0},qu_{q}\rangle$ (8)

where $\langle,\rangle$ is the $L_{2}$ inner product on the PDE domain.
Because of the dependence of $u_{q}$ on $q$, equation (8) is generally
nonlinear. The standard linearization is given by the famous Born
approximation in which one assumes that $u_{q}\approx u_{0}$, which is
accurate only for small $q$. In this paper we suggest to use the more accurate
approximation (6). This yields

$F_{q}(s)-F_{0}(s)\approx-\langle u_{0},q\tilde{u}_{q}\rangle.$ (9)

Indeed, the internal solution $\tilde{u}_{q}$ depends on $q$, however, it can
be precomputed directly from the data without knowing $q$. After that, (9)
becomes linear with respect to $q$. As we shall see, formula (9) relies on the
Lanczos algorithm for the computation of $T_{q}$, which is why we call (9) the
Lippmann-Schwinger-Lanczos equation.

As a by-product, the Lippmann-Schwinger-Lanczos approach resolves a number of
problems arising in the data-driven ROM framework. Recently in [2], the map
between $T_{q}-T$ and $q$ was computed via training and required multiple
solutions of forward problems for different $q$. The Lippmann-Schwinger-
Lanczos approach yields an explicit expression for this map. Furthermore, this
approach also allows for a more rigorous derivation of the back-projection
algorithm, and suggests natural extensions to more general data sets.

Finally, we should point out, that there are known approaches using internal
solutions in the Lippmann-Schwinger framework, of which Marchenko redatuming
is the closest to the suggested method here, for example, see [10]. The main
difficulty with Marchenko redatuming lies in the accurate approximation of the
inverse scattering transform in the continuous setting, and in the evaluation
of some integrals. The Lanczos based ROM approach here can be interpreted as a
linear-algebraic realization of the Marchenko-Gelfand-Levitan method [14] with
a data-driven spectral discretization, and as such promises the best possible
accuracy for a given data-set. As we shall see in our numerical experiments,
with as little as 6 Laplace frequencies and 8 source-receiver positions, our
approach produces results which are indistinguishable from the Lippmann-
Schwinger formulation using the exact internal solutions.

This paper is organized as follows. In Section 2 we describe the entire
process in detail for a one dimensional, single input single output (SISO)
problem. This includes the construction of the ROM from the data, the Lanczos
orthogonalization process, the generation of the internal solution and its use
in the Lippmann-Schwinger equation. The generalization of this process to
multiple input multiple output (MIMO) problems in higher dimensions is
described in Section 3. Section 4 contains numerical experiments, and in the
appendix we describe how the Lippmann-Schwinger Lanczos algorithm is related
to other approaches to inversion using the data driven ROM.

##  2 One dimensional SISO problem

We begin this work with the one dimensional problem, since the presentation is
simpler, and the ideas extend naturally to higher dimensions. In the first
subsection we describe the problem setup and in the second subsection we show
how one constructs the ROM from the data. In the third and fourth subsections
we discuss tridiagonalization of the ROM and generation of the internal
solution from data only. In the last subsection we show how to use the
internal solutions in the Lippmann-Schwinger equation in order to solve the
fully nonlinear inverse problem.

###  2.1 Description of the SISO problem

We start by considering the single input single output (SISO) inverse problem
in one dimension

$-\frac{d^{2}u(x,\lambda)}{dx^{2}}+q(x)u(x,\lambda)+\lambda
u(x,\lambda)=g(x)\quad\frac{du}{dx}|_{x=0}=0,\ \frac{du}{dx}|_{x=L}=0,$ (10)

where $0<L\leq\infty$. The source $g(x)$ is assumed to be a compactly
supported real distribution localized near the origin, for example, roughly
speaking, $g=\delta(x-\epsilon)$ with small $\epsilon>0$. (Note that when $g$
is a delta function at the origin this corresponds to an inhomogeneous Neumann
boundary condition.) We can write the solution formally as

$u=\left(-\frac{d^{2}}{dx^{2}}+{q}I+\lambda I\right)^{-1}g$ (11)

where the inverse is understood to correspond to Neumann boundary conditions.
Consider $\lambda_{j}\in{\mathbb{C}}\setminus{\mathbb{R}}_{-}$,
$j=1,\ldots,m$, with $\Im\lambda_{j}\geq 0$. The SISO transfer function is
then

$F(\lambda)=\int_{0}^{L}g(x)u(x,\lambda)dx=\langle g,u\rangle=\langle
g,\left(-\frac{d^{2}}{dx^{2}}+{q}I+\lambda I\right)^{-1}g\rangle$ (12)

where

$\langle w,v\rangle=\int_{0}^{L}\bar{w}(x)v(x)dx$

is the Hermitian inner product on $L^{2}(0,L)$. For $2m$ real data points,
that is, for $\Im\lambda_{j}=0$, we consider the data

$F(\lambda)|_{\lambda=\lambda_{j}}\in{\mathbb{R}},\ \ \
\frac{dF(\lambda)}{d\lambda}|_{\lambda=\lambda_{j}}\in{\mathbb{R}}\ \ \
\mbox{for}\ \ j=1,\ldots,m.$ (13)

For complex data points, we consider just

$F(\lambda)|_{\lambda=\lambda_{j}}\in{\mathbb{C}}\ \ \ \mbox{for}\ \
j=1,\ldots,m.$ (14)

Note that the complex case is equivalent to also having data at
$\bar{\lambda}_{j}$, since $F(\overline{\lambda})=\overline{F}(\lambda)$ from
the fact that $g$ and $q$ are real. From this one can see that for
$\lambda_{j}$ close to real,

$\displaystyle\frac{dF(\lambda)}{d\lambda}|_{\lambda=\Re{\lambda_{j}}}$
$\displaystyle=$ $\displaystyle\lim_{\Im\lambda_{j}\to
0}{F(\lambda_{j})-F(\overline{\lambda_{j}})\over{\lambda_{j}-\overline{\lambda_{j}}}}$
(15) $\displaystyle=$ $\displaystyle\lim_{\Im\lambda_{j}\to 0}{\Im
F(\lambda_{j})\over{\Im{\lambda_{j}}}}$

and hence the real data can be viewed as the natural extension of complex
data. The SISO inverse problem is then to determine $q(x)$ in (10) from the
data (14) or (13).

###  2.2 Construction of the data-driven ROM

We will treat $u(x,\lambda_{j})$ as a continuous basis function (or it can be
viewed as an infinite dimensional vector column) and consider the projection
subspace

${\bf
V}=\mbox{span}\\{u_{1}(x)=u(x,\lambda_{1}),\ldots,u_{m}(x)=u(x,\lambda_{m})\\}.$

We define the data-driven ROM as the Galerkin system for this subspace

$Sc(\lambda)+\lambda Mc(\lambda)=b$ (16)

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

$S_{ij}=\langle{u}^{\prime}_{i},u^{\prime}_{j}\rangle+\langle
q{u}_{i},u_{j}\rangle$

and mass matrix $M$ given by

$M_{ij}=\langle{u}_{i},u_{j}\rangle.$

The right hand side $b\in{\mathbb{C}}^{m}$ has components

$b_{j}=\langle{u}_{j},g\rangle,$

and the Galerkin solution for the system is determined by the vector valued
function of $\lambda$, $c(\lambda)\in{\mathbb{C}}^{m}$. Note that $c(\lambda)$
corresponds to coefficients of the solution with respect to the above basis of
exact solutions. The matrices $S$ and $M$ are obtained from matching
conditions that ensure that the ROM transfer function

$\tilde{F}(\lambda)=b^{*}c$

matches the data. For real $\lambda_{i}$ this is the same as in [6]. For
complex $\\{\lambda_{j}\\}^{m}_{j=1}$ multiplying (10) for
$\lambda=\lambda_{i}$ by $\bar{u}_{j}$ in the inner product
$\langle\cdot,\cdot\rangle$ and integrating by parts we obtain

$\bar{S}_{ij}+\lambda_{i}\bar{M}_{ij}=\bar{F}(\lambda_{j}).$ (17)

Similarly, multiplying the conjugate of (10) for $\lambda=\bar{\lambda}_{j}$
by $u_{i}$ in the inner product $\langle\cdot,\cdot\rangle$, we obtain

$S_{ij}+\lambda_{j}M_{ij}=b_{i},$ (18)

where $b_{i}=\bar{F}(\lambda_{i})$. Subtracting the complex conjugate of (17)
from (18) we obtain the expression for the mass matrix

$M_{ij}=\frac{\bar{F}(\lambda_{i})-F(\lambda_{j})}{\lambda_{j}-\bar{\lambda}_{i}}.$
(19)

Similarly, the elements of stiffness matrix have the form

$S_{ij}=\frac{F(\lambda_{j})\lambda_{j}-\bar{F}(\lambda_{i})\bar{\lambda}_{i}}{\lambda_{j}-\bar{\lambda}_{i}}.$
(20)

Furthermore, the solution to (10) is close to its Galerkin projection

$u(\lambda)\approx\tilde{u}(\lambda)={V}c(\lambda)={V}(S+\lambda M)^{-1}b$

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

${V}=(u_{1},\ldots,u_{m}).$

Then, using that $\bar{\tilde{F}}(\lambda)=\tilde{F}(\bar{\lambda})$, the
following proposition follows immediately.

######  Proposition 1.

Assume that $\Im\lambda_{j}\neq 0$ for $j=1,\ldots,m$. Then the Galerkin
projection of the solution of (10)

$\tilde{u}(\lambda)=Vc(\lambda)=V(S+\lambda M)^{-1}b$

is exact at $\lambda=\lambda_{j}$ ,

$\tilde{u}(\lambda_{j})=u(\lambda_{j}),$

and hence

$\tilde{F}(\lambda_{j})=b^{*}(S+\lambda_{j}M)^{-1}b=F(\lambda_{j})$

and

$\tilde{F}(\bar{\lambda}_{j})=b^{*}(S+\bar{\lambda}_{j}M)^{-1}b=F(\bar{\lambda}_{j})$

for $j=1,\ldots,m$.

The corresponding proposition for $\lambda_{j}$ real was shown in [6].

###  2.3 Lanczos tridiagonalization

In the next step, we orthogonalize the above basis of exact solutions 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{\mathbb{R}}^{m\times m}$ and
$M$-orthonornal Lanczos vectors $q_{i}\in{\mathbb{R}}^{m}$, such that

$AQ=QT,\qquad Q^{*}MQ=I,$ (21)

where

$Q=[q_{1},q_{2},\ldots,q_{m}]\in{\mathbb{R}}^{m\times m},$

and

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

This resulting new basis $VQ$ will be orthonormal with respect to the
$L^{2}(0,L)$ inner product. The Galerkin projection of (10) and its transfer
function can then be written in Lanczos coordinates as

$\tilde{u}(\lambda)=\sqrt{b^{*}M^{-1}b}VQ(T+\lambda I)^{-1}e_{1},$ (22)
$\tilde{F}(\lambda)=(b^{*}M^{-1}b)e_{1}^{*}(T+\lambda I)^{-1}e_{1}.$ (23)

It is known that when $\Sigma\in{\mathbb{C}}\setminus{\mathbb{R}}_{-}$ is
compact, for any sequence $\lambda_{j}\in\Sigma$, $j=1,\ldots,m$ we have that
for any $\lambda\in\Sigma$ the corresponding Galerkin solution converges
exponentially $\tilde{u}\rightarrow u$ as $m\rightarrow\infty$ with a uniform
linear rate.

###  2.4 Internal solutions

We assume that we do not know $q$, yet we want to compute $u(x,\lambda)$
directly from the data. Recall that the output of the Lanczos algorithm,
tridiagonal $T$ and change of basis $Q$, were obtained from data only.
However, we do not know the original basis of exact solutions $V$. We propose
to approximate $u$ internally, as was done in [6], by replacing the unknown
orthogonalized internal solutions $VQ$ with orthogonalized background
solutions $V_{0}Q_{0}$ corresponding to background $q_{0}=0$. Here $V_{0}$ is
the row vector containing the basis of background solutions

${V}=(u^{0}_{1},\ldots,u^{0}_{m})$

to (10) corresponding to $q=q_{0}=0$ and the same spectral points
$\lambda=\lambda_{1},\ldots\lambda_{m}$, and $Q_{0}$ is computed from Lanczos
orthogonalization of the background ROM. That is, one can compute an
approximation to $u(x,\lambda)$ using

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

which is obtained from data only.

######  Conjecture 2.

Assume that $V_{0}$ and $Q_{0}$, are the solution basis and Lanczos matrix
corresponding to the background $q_{0}=0$. Then for all $\lambda\in\Sigma$ we
have that

$u=\lim_{m\to\infty}{\mathbf{u}}=\sqrt{b^{*}M^{-1}b}V_{0}Q_{0}(T+\lambda
I)^{-1}e_{1}.$ (24)

that is, for any fixed $\lambda\in\Sigma$, as the number of data points in
$\Sigma$ approaches infinity, the data generated internal solution
${\mathbf{u}}$ converges to the true solution $u$.

###  2.5 Nonlinear inverse problem

We now use the Lippmann-Schwinger formulation for the solutions $u$ to solve
the nonlinear inverse problem by using the internal solutions described above.
From (12) we obtain

$\displaystyle F_{0}(\lambda_{j})-F(\lambda_{j})=\langle
g,\left(-\frac{d^{2}}{dx^{2}}+\lambda I\right)^{-1}g\rangle-\langle
g,\left(-\frac{d^{2}}{dx^{2}}+{q}I+\lambda I\right)^{-1}g\rangle$
$\displaystyle=\langle\left(-\frac{d^{2}}{dx^{2}}+{q}I+\bar{\lambda}I\right)^{-1}g,q\left(-\frac{d^{2}}{dx^{2}}+\lambda
I\right)^{-1}g\rangle$ (25)

that is,

$F_{0}(\lambda_{j})-F(\lambda_{j})=\int
u_{0}^{*}(x,\lambda_{j})u(x,\lambda_{j})q(x)dx,\\\ \qquad j=1,\ldots,m$ (26)

which can also be seen as a direct consequence of the usual Lippmann-Schwinger
formulation with the Green’s function. Here $u_{0}$ corresponds to the
solution to the background problem with $q_{0}=0$. If all $\lambda_{j}$ are
complex, then (26) yields $2m$ real equations for $q$, which is the same as
the number of data points. For real $\lambda_{j}\in\mathbb{R}$ we also have
$2m$ real equations

$\displaystyle F_{0}(\lambda_{j})-F(\lambda_{j})=\langle
u_{0},{q}u\rangle=\int u^{*}_{0}(x,\lambda_{j})u(x,\lambda_{j})q(x)dx,$ (27)
$\displaystyle\frac{d}{d\lambda}(F_{0}-F)|_{\lambda=\lambda_{j}}=\int\frac{d}{d\lambda}[u^{*}_{0}(x,\lambda)u(x,\lambda_{j})]_{\lambda=\lambda_{j}}q(x)dx$

for $j=1,\ldots,m$. Of course, the internal solutions $u(x,\lambda_{j})$ and
their derivatives with respect to $\lambda$ are unknown, and they depend on
$q$, so the system (26-27) is nonlinear with respect to $q$. A Born
linearization replaces $u_{0}$ with $u$, however, it is only accurate for
small $q$. Using Conjecture 2, we replace $u(x,\lambda)$ in (26-27) with its
approximation

$u\approx{\mathbf{u}}=\sqrt{b^{*}M^{-1}b}V_{0}Q_{0}(T+\lambda I)^{-1}e_{1}.$

For example, for the case of complex $\lambda_{j}$ we can write the new system
for $q$ as

$\delta{\mathbf{F}}=\langle W,q\rangle$ (28)

where

$\delta{\mathbf{F}}=[F_{0}(\lambda_{1})-F(\lambda_{1}),\ldots,F_{0}(\lambda_{m})-F(\lambda_{m})]^{*}\in{\mathbb{C}}^{m},$

and

$W=[{\mathbf{u}}^{*}(x,\lambda_{1})u_{0}(x,\lambda_{1}),\ldots,{\mathbf{u}}^{*}(x,\lambda_{m})u_{0}(x,\lambda_{m})]$

is an $m$-dimensional vector complex valued functions on $(0,L)$. By
construction, ${\mathbf{u}}$ can be directly computed from the data without
knowing $q$, by using (24), thus making nonlinear system (28) linear! A
reasonable regularization of (28) is to restrict $q$ to the dominant left
singular vectors of $W$. We will refer to (28) as a Lippmann-Schwinger-Lanczos
system.

##  3 Multidimensional MIMO problem

###  3.1 Description of the MIMO problem

We consider the boundary value problem on $\Omega\in{\mathbb{R}}^{d}$ for

$-\Delta u^{(r)}+qu^{(r)}+\lambda
u=g^{(r)},\quad\frac{du}{d\nu}\large|_{\partial\Omega}=0,\ r=1,\ldots,K,$ (29)

where $g^{(r)}$ are localized sources, e.g., boundary charge distributions,
supported near or at an accessible part $S$ of $\partial\Omega$. Let

$G=[g^{(1)},g^{(2)},\ldots,g^{(K)}]\in{\mathbb{R}}^{\infty\times k}$

and

$U=[u^{(1)},u^{(2)},\ldots,u^{(K)}]\in{\mathbb{R}}^{\infty\times k},$

again understood as vectors of continuous functions. Then the multiple-input
multiple output (MIMO) transfer function is a matrix valued function of
$\lambda$

$F(\lambda)=\langle G,U\rangle=\langle
G,\left(-\frac{d^{2}}{dx^{2}}+\mathop{\operator@font diag}\nolimits{q}+\lambda
I\right)^{-1}G\rangle\in{\mathbb{C}}^{K\times K},$ (30)

For real positive $\lambda$ we will have symmetric positive definite
$F(\lambda)$. As in the SISO case, we consider the inverse problem with data
given by $2m$ real symmetric $K\times K$ matrices, that is, for
$\Im\lambda_{j}=0$ our data are

$F(\lambda)|_{\lambda=\lambda_{j}}\in{\mathbb{R}},$

and

$\frac{F(\lambda)}{d\lambda}|_{\lambda=\lambda_{j}}\in{\mathbb{R}},$

and $F(\lambda)|_{\lambda=\lambda_{j}}\in{\mathbb{C}}$ otherwise.

###  3.2 Construction of the MIMO data-driven ROM

We consider the $mK$ dimensional projection subspace

$\bf{V}=\mbox{span}\\{U_{1}(x)=U(x,\lambda_{1}),\ldots,U_{m}(x)=U(x,\lambda_{m})\\}.$

We then define the MIMO data-driven ROM as

$(S+\lambda M)C(\lambda)=B$ (31)

where $S,M\in{\mathbb{C}}^{mK\times mK}$ are Hermitian positive definite
matrices, $B\in{\mathbb{R}}^{mK\times K}$, and $C\in{\mathbb{C}}^{mK\times K}$
is a matrix valued function of $\lambda$, again corresponding to coefficients
of the solution with respect to the above basis of exact solutions. Stiffness
and mass matrices can be written as block extensions of their counterparts in
(16), with blocks given by

$S=(S_{ij}=\langle{U^{\prime}}_{i},U^{\prime}_{j}\rangle)+\langle
q{U}_{i},U_{j}\rangle)$

and

$M=(M_{ij}=\langle{U}_{i},U_{j}\rangle).$

Once again $S$ and $M$ are obtained from imposing the conditions that the ROM
transfer function

$\tilde{F}(\lambda)=B^{*}C(\lambda)$

matches the data.

###  3.3 Block-Lanczos tridiagonalization

We then run $m$ steps of the $M$ symmetric block-Lanczos algorithm with matrix

$A=M^{-1}S$

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

$T\in{\mathbb{R}}^{Km\times Km}$

with $K\times K$ blocks, and $M$-orthonormal Lanczos block-vectors
$q_{i}\in{\mathbb{R}}^{mK\times K}$. From this we obtain the block counterpart
of (21)

$Q=[q_{1},q_{2},\ldots,q_{m}]\in{\mathbb{R}}^{mK\times mK},$

where

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

The state solution then can be written in Lanczos coordinates as

$\tilde{u}(\lambda)=\sqrt{B^{*}M^{-1}B}VQ(T+\lambda I)^{-1}E_{1},$ (32)

where $E_{1}\in{\mathbb{R}}^{mK\times K}$ consists of the first $K$ columns of
$I\in{\mathbb{R}}^{mK\times mK}$.

###  3.4 Internal solutions

The SISO result is generalizable to MIMO case, and can be vaguely formulated
as following conjecture

######  Conjecture 3.

Assume we have sequence of input sets

$G_{k}=[g_{k}^{(1)},g_{k}^{(2)},\ldots,g_{k}^{(K)}],$

such that as $k\rightarrow\infty$ the range of $G_{k}$ approximates surface
distributions at the observable boundary $S\subset\partial\Omega$. Then for
all $\lambda\in\Sigma$

$U=\lim_{m,K\to\infty}{\mathbf{U}}=\sqrt{B^{*}M^{-1}B}V_{0}Q_{0}(T+\lambda
I)^{-1}E_{1}.$ (33)

###  3.5 Nonlinear MIMO inverse problem

From (12) we obtain

$F_{0}(\lambda_{j})-F(\lambda_{j})=\langle G,\left(-\Delta+\lambda
I\right)^{-1}G\rangle-\langle G,\left(-\Delta+qI+\lambda
I\right)^{-1}G\rangle\\\ =\langle\left(-\Delta+{q}I+\lambda
I\right)^{-1}G,q\left(-\Delta+\lambda I\right)^{-1}G\rangle\\\ =\int
U^{*}_{0}(x,\lambda_{j})q(x)U(x,\lambda_{j})dx,$ (34)

that is,

$F_{0}(\lambda_{j})-F(\lambda_{j})=\int
U^{*}_{0}(x,\lambda_{j})q(x)U(x,\lambda_{j})dx,\qquad j=1,\ldots,m.$ (35)

Here again the subscript $0$ corresponds to background solutions with $q=0$.
If all $\lambda_{j}$ are complex, then (26) gives $2m$ real equations for $q$,
equal to the number of data points. For $\lambda_{j}\in R$ we will have

$\displaystyle(F_{0}-F)|_{\lambda_{j}}$ $\displaystyle=$ $\displaystyle\int
U^{*}_{0}(x,\lambda_{j})q(x){U}(x,\lambda_{j})dx,$
$\displaystyle\frac{d}{d\lambda}(F_{0}-F)|_{\lambda_{j}}$ $\displaystyle=$
$\displaystyle\int\frac{d}{d\lambda}(U^{*}_{0}(x,\lambda){U}(x,\lambda))|_{\lambda=\lambda_{j}}q(x)dx.$
(36)

Similar to the SISO case, by using Conjecture 3 we replace $U(x,\lambda)$ with
its approximation ${\mathbf{U}}$ in (35) and (3.5). Again, like in the SISO
case, precomputing ${\mathbf{U}}$ via the data-driven algorithm (33) will
yield the linear system for $q$:

$\displaystyle(F_{0}-F)|_{\lambda_{j}}$ $\displaystyle=$ $\displaystyle\int
U^{*}_{0}(x,\lambda_{j})q(x)\mathbf{U}(x,\lambda_{j})dx,$
$\displaystyle\frac{d}{d\lambda}(F_{0}-F)|_{\lambda_{j}}$ $\displaystyle=$
$\displaystyle\int\frac{d}{d\lambda}(U^{*}_{0}(x,\lambda)\mathbf{U}(x,\lambda))|_{\lambda=\lambda_{j}}q(x)dx.$
(37)

##  4 Numerical Experiments

In this section we present numerical results for reconstructing a 2D two-bump
media. For simplicity, we considered noiseless case only. In the presence of
noise, constructing the data-driven ROM may become unstable and requires
regularization. This problem was resolved in [9] by optimal SVD-based
truncation of the unstable part. However, to avoid additional complications,
we skipped that part in this work. In the first experiment we considered two
low contrast bump inclusions (see Fig.1, top left). The measurement setup
mimics one from medical imaging, i.e. we had $K=8$ sources located on the
boundary, two on each side (see red crosses at Fig.1, top left). We used $s=6$
positive spectral values. Here and below, to obtain $q(x)$, we solve ill-
conditioned linear system (3.5) by projecting it onto its dominant
eigenvectors. On the top right of Fig.1 we plotted the reconstruction when
actual true internal solution ${U}(x,\lambda))$ is available. This is
typically not the case, so we called it ’Cheated IE’. The usual Born
linearization result (when we replace ${U}(x,\lambda))$ in (3.5) with
${U}_{0}(x,\lambda))$) is shown on bottom left of Fig.1. Finally, the
reconstruction using our approach is plotted on the bottom right of Fig.1. As
one can observe, the Born approximation results in multiple artifacts and
produces a significantly worse reconstruction compared to our approach and
’Cheated IE’.

In our second experiment we considered the same two bumps but with increased
contrast (see Fig.2, top left). The source locations were the same as in
previous experiment, however, out of $s=6$ spectral values we took one
negative and 5 positive. Similar to our first experiments we plotted the
results obtained using the ’Cheated IE’, the Born approximation and our
approach (see on the top right, bottom left and bottom right of Fig.2,
respectively). Because of higher contrast, the Born linearization failed to
produce a meaningful reconstruction. However, both the ’Cheated IE’ and our
approach still performed well.

Fig. 1: Experiment 1: True medium (top left) and its reconstructions using
’Cheated IE’ (top right), Born linearization (bottom left) and our approach
(bottom right) Fig. 2: Experiment 2: True medium (top left) and its
reconstructions using ’Cheated IE’ (top right), Born linearization (bottom
left) and our approach (bottom right)

In our third experiment we again considered a medium with 2 bumps with
increased contrast (see Fig.3, top left), however, with data acquisition
similar to surface geophysics or radars. That is, we probed the medium with
$K=12$ sources located on one (upper) side $y=-1$ only (see red crosses at
Fig.3, top left). For better aperture, the lateral extent of the acquisition
boundary was three times larger than the depth of the domain. In this example
we considered $s=6$ positive frequencies. We compared the results of the
’Cheated IE’, the Born approximation and our approach (see on the top right,
bottom left and bottom right of Fig.3, respectively). Similar to the previous
example, the reconstruction via the Born approximant is totally unsatisfactory
because of high contrast of the inclusions. At the same time, our approach,
along with ’Cheated IE’, performs well.

Fig. 3: Experiment 3: True medium (top left) and its reconstructions using
’Cheated IE’ (top right), Born linearization (bottom left) and our approach
(bottom right)

##  Appendix A Generalizations and extensions

In this appendix we show how the system (28) is connected to other algorithms
for extracting the unknown $q$ from the ROM, and describe how it lends itself
naturally to extensions to other data sets.

###  A.1 Connection to the back-projection algorithm

The back-projection algorithm is based on the idea that the perturbed and
background operators differ only in the lower order term, so that their
difference should approximate $q$. Using the reduced models and replacing $VQ$
with $V_{0}Q_{0}$, the difference

$A_{q}-A_{0}=qI$

is approximated by the operator $\mathcal{Q}$

$A_{q}-A_{0}\approx\mathcal{Q}$ (38)

where

$\mathcal{Q}w=V_{0}{Q}_{0}(T-T_{0})\int(V_{0}{Q}_{0})^{*}w.$ (39)

Recall that $V_{0}Q_{0}$ is the row vector of orthogonalized background
solutions. Here $T$ and $T_{0}$ are obtained via data matching conditions, for
example via time-domain matching, or the frequency matching as described
above.

To see how this relates to Lippmann Schwinger, suppose that one uses the
operator approximation (39) for $q$ and replaces the right hand side of (28)
with $\langle u_{0},\mathcal{Q}{\mathbf{u}}\rangle$. We get that

$\langle
u_{0},\mathcal{Q}{\mathbf{u}}\rangle=\int(V_{0}Q_{0})(x)(T-T_{0})\int(V_{0}Q_{0})^{*}(x^{\prime}){\mathbf{u}}(x^{\prime})dx^{\prime}u^{*}_{0}(x)dx\\\
=(b_{0}^{*}M_{0}^{-1}b_{0})^{1/2}(b^{*}M^{-1}b)^{1/2}\int\int
e_{1}^{*}(T_{0}+\lambda_{j}I)^{-1}\cdot\\\
\cdot(V_{0}Q_{0})^{*}(x){(V_{0}Q_{0})}(x)(T-T_{0})(V_{0}Q_{0})^{*}(x^{\prime})(V_{0}Q_{0})(x^{\prime})(T+\lambda_{j}I)^{-1}e_{1}dx^{\prime}dx\\\
=(b_{0}^{*}M_{0}^{-1}b_{0})^{1/2}(b^{*}M^{-1}b)^{1/2}e_{1}^{*}(T_{0}+\lambda_{j}I)^{-1}(T-T_{0})(T+\lambda_{j}I)^{-1}e_{1}$
(40)

from othogonality of the basis. At the same time, the left hand side of (28)
is

$\displaystyle F_{0}(\lambda_{j})-F(\lambda_{j})$ $\displaystyle=$
$\displaystyle
g^{T}(\Delta+\lambda_{j}I)^{-1}g-g^{T}(\Delta+{q}I+\lambda_{j}I)^{-1}g$
$\displaystyle=$
$\displaystyle{b_{0}^{*}M_{0}^{-1}b_{0}}e_{1}^{*}(T_{0}+\lambda_{j}I)^{-1}e_{1}-{b^{*}M^{-1}b}e_{1}^{*}(T+\lambda_{j}I)^{-1}e_{1}$

since the ROM matches the data exactly. By definition,
$b_{0}^{*}M_{0}^{-1}b_{0}$ and $b^{*}M^{-1}b$ are the norms of projections of
$g$ on $V_{0}$ and $V$ respectively, which should become close as the
dimension grows, that is, we expect that

$\lim_{m\to\infty}b_{0}^{*}M_{0}^{-1}b_{0}=b^{*}M^{-1}b.$ (41)

We then obtain

$\displaystyle F_{0}(\lambda_{j})-F(\lambda_{j})$ $\displaystyle\approx$
$\displaystyle(b^{*}M^{-1}b)e_{1}^{*}(T_{0}+\lambda_{j}I)^{-1}(T-T_{0})(T+\lambda_{j}I)^{-1}e_{1}$
(42) $\displaystyle\approx$ $\displaystyle\langle
u_{0},\mathcal{Q}{\mathbf{u}}\rangle$ (43)

from (40) and (41). Thus we have that, modulo (41), (39) satisfies the part of
system (28) corresponding to $F_{0}(\lambda_{j})-F(\lambda_{j})$ for
$j=1,\ldots,m$. Hence the system (28) would be valid for
$\lambda_{j},\bar{\lambda}_{j}$, $j=1,\ldots,m$ if we replace $q$ with its
operator approximation (39). The extension to the derivative data in the real
case can be obtained via limiting transition. In summary, the operator
$\mathcal{Q}$ essentially satisfies the Lippmann-Schwinger-Lanczos equation.

The back projection algorithm extracts a scalar version of (39) in a natural
way. Indeed, the simplest back projection reconstruction step is to take

$q(x)=\int
q(x^{\prime})\delta(x-x^{\prime})dx^{\prime}\approx\int\mathcal{Q}\delta(x-x^{\prime})dx^{\prime}\\\
=\int V_{0}(x^{\prime}){Q}_{0}(T-T_{0})(V_{0}(x){Q}_{0})^{*}dx^{\prime}.$ (44)

It can be show that the probing function (44) uses approximations of
$\delta(x-x^{\prime})$ on the solution snapshot space. Thanks to Conjecture 2,
this can be defined via projection

$\tilde{\delta}(x-x^{\prime})=(V_{0}(x){Q}_{0})(V_{0}(x^{\prime}){Q}_{0})^{*},$

and is related to a known approach of imaging diagonal operators, e.g., [16].

The back-projection approach of [11] instead uses a sharper point spread
function, obtained by the squaring and re-normalization of the delta function:

$\frac{\tilde{\delta}(x-x^{\prime})^{2}}{\sqrt{\int\tilde{\delta}(x-x^{\prime})^{2}dx^{\prime}}},$

which yields

$q(x)\approx\frac{\int\tilde{\delta}(x-x^{\prime})\mathcal{Q}\tilde{\delta}(x-x^{\prime})dx^{\prime}}{\int\tilde{\delta}(x-x^{\prime})^{2}dx^{\prime}}=\frac{V_{0}(x){Q}_{0}(T-T_{0})(V_{0}(x){Q}_{0})^{*}}{\int\tilde{\delta}(x-x^{\prime})^{2}dx^{\prime}}.$
(45)

###  A.2 Linear mapping from ROM to potential

From (40), (42) and Conjecture 2 we obtain

$s_{0}(\lambda_{j})^{*}(T-T_{0})s(\lambda_{j})\approx
F_{0}(\lambda_{j})-F(\lambda_{j})\approx\int
u^{*}_{0}(x,\lambda_{j}){\mathbf{u}}(x,\lambda_{j})q(x)dx,$ (46)

where

$s_{0}(\lambda_{j})=(b_{0}^{*}M_{0}^{-1}b_{0})^{1/2}e_{1}^{*}(T_{0}+\lambda_{j}I)^{-1}\in{\mathbb{C}}^{m},$
$s(\lambda_{j})=(b^{*}M^{-1}b)^{1/2}e_{1}^{*}(T+\lambda_{j}I)^{-1}\in{\mathbb{C}}^{m},$

and where ${\mathbf{u}}$ is given by (24). Note that $s$ and $s_{0}$ are data
driven.

######  Theorem 4.

Assuming that $s,s_{0}$ and ${\mathbf{u}}$ are precomputed from the data, (46)
yields a linear relationship between $q$ and $(T-T_{0})$ that uniquely defines
the latter for complex $\lambda_{j}$, $j=1,\ldots,m$. For real $\lambda_{j}$
for uniqueness we need to add derivatives at $\lambda_{j}$

###### Proof.

The result follows from uniqueness of [m-1/m] Padé from $2m$ matching
conditions. ∎

###  A.3 Extension to general ROMs

The linear mapping approach is not limited to the ROM obtained via frequency
matching. For any reduced order model we can assume the relationship

$s_{0}(\lambda)^{*}(T-T_{0})s(\lambda)\approx(F_{0}-F)(\lambda)\approx\\\ \int
u^{*}_{0}(x,\lambda){\mathbf{u}}(x,\lambda)q(x)dx,$ (47)

for any $\lambda\in\Sigma\subset\mathbb{C}\setminus{\mathbb{R}}_{-}$, where
$\Sigma$ is a compact complex domain. Then $(T-T_{0})$ can be uniquely defined
via matching (47) for any $m$ complex $\lambda_{j}$ in that domain, or $m$
real frequencies and their derivatives. Furthermore, (47) can be transformed
to the time domain where one will obtain convolution equations.

Acknowledgements V. Druskin was partially supported by AFOSR grant FA
955020-1-0079. S. Moskow was partially supported by NSF grants DMS-1715425 and
DMS-2008441.

## References

  * [1] Christopher Beattie and Serkan Gugercin.  Model Reduction and Approximation, chapter 7: Model Reduction by Rational Interpolation, pages 297–334.  SIAM, 2017.
  * [2] L. Borcea, V. Druskin, A. Mamonov, M. Zaslavsky, and J. Zimmerling.  Reduced order model approach to inverse scattering.  SIAM Journal on Imaging Sciences, 13(2):685–723, 2020.
  * [3] L. Borcea and J. Zimmerling.  Personal communication.
  * [4] Liliana Borcea and Vladimir Druskin.  Optimal finite difference grids for direct and inverse Sturm-Liouville problems.  Inverse Problems, 18(4):979–1001, 2002.
  * [5] Liliana Borcea, Vladimir Druskin, Fernando Guevara Vasquez, and Alexander V. Mamonov.  Resistor network approaches to electrical impedance tomography.  Inverse Problems and Applications: Inside Out II, Math. Sci. Res. Inst. Publ, 60:55–118, 2011.
  * [6] Liliana Borcea, Vladimir Druskin, Alexander V. Mamonov, Shari Moskow, and Mikhail Zaslavsky.  Reduced order models for spectral domain inversion: embedding into the continuous problem and generation of internal data.  Inverse Problems, 36(5), 2020.
  * [7] Liliana Borcea, Vladimir Druskin, Alexander V. Mamonov, and Mikhail Zaslavsky.  A model reduction approach to numerical inversion for a parabolic partial differential equation.  Inverse Problems, 30(12):125011, 2014.
  * [8] Liliana Borcea, Vladimir Druskin, Alexander V Mamonov, and Mikhail Zaslavsky.  Untangling nonlinearity in inverse scattering with data-driven reduced order models.  Inverse Problems, 34(6):065008, 2018.
  * [9] Liliana Borcea, Vladimir Druskin, Alexander V. Mamonov, and Mikhail Zaslavsky.  Robust nonlinear processing of active array data in inverse scattering via truncated reduced order models.  Journal of Computational Physics, 381:1–26, 2019.
  * [10] Leon Diekmann and I. Vasconcelos.  Imaging with the exact linearised Lippmann-Schwinger integral by means of redatumed in-volume wavefields.  SEG Technical Program Expanded Abstracts, 2020.
  * [11] V. Druskin, A.V. Mamonov, and M. Zaslavsky.  A nonlinear method for imaging with acoustic waves via reduced order model backprojection.  SIAM Journal on Imaging Sciences, 11(1):164–196, 2018.
  * [12] V. Druskin and S. Moskow.  Three-point finite-difference schemes, padé and the spectral galerkin method. i. one-sided impedance approximation.  Mathematics of computation, 71(239):995–1019, 2002.
  * [13] Vladimir Druskin and Leonid Knizhnerman.  Gaussian spectral rules for the three-point second differences: I. A two-point positive definite problem in a semi-infinite domain.  SIAM Journal on Numerical Analysis, 37(2):403–422, 1999.
  * [14] Vladimir Druskin, Alexander V. Mamonov, Andrew E. Thaler, and Mikhail Zaslavsky.  Direct, nonlinear inversion algorithm for hyperbolic problems via projection-based model reduction.  SIAM Journal on Imaging Sciences, 9(2):684–747, 2016.
  * [15] IS Kac and MG Krein.  On the spectral functions of the string.  Amer. Math. Soc. Transl, 103(2):19–102, 1974.
  * [16] Howard W Levinson and Vadim A Markel.  Solution of the nonlinear inverse scattering problem by t-matrix completion. i. theory.  Physical Review E, 94(4):043317, 2016.