# Oversampled Adaptive Sensing via
a Predefined Codebook

Ali Bereyhi, Saba Asaad, and Ralf R. Müller  Friedrich-Alexander Universität
Erlangen-Nürnberg, Germany
{ali.bereyhi, saba.asaad<EMAIL_ADDRESS>This work has been presented
in 2021 IEEE International Symposium on Joint Communications & Sensing (JC&S).
The link to the final version in the proceedings will be available later.

###### Abstract

Oversampled adaptive sensing (OAS) is a Bayesian framework recently proposed
for effective sensing of structured signals in a time-limited setting. In
contrast to the conventional blind oversampling, OAS uses the prior
information on the signal to construct posterior beliefs sequentially. These
beliefs help in constructive oversampling which iteratively evolves through a
sequence of time sub-frames.

The initial studies of OAS consider the idealistic assumption of full control
on sensing coefficients which is not feasible in many applications. In this
work, we extend the initial investigations on OAS to more realistic settings
in which the sensing coefficients are selected from a predefined set of
possible choices, referred to as the codebook. We extend the OAS framework to
these settings and compare its performance with classical non-adaptive
approaches.

###### Index Terms:

Oversampled adaptive sensing, sparse recovery, Bayesian inference, stepwise
regression.

OAS

    oversampled adaptive sensing
AWGN

    additive white Gaussian noise
i.i.d.

    independent and identically dis-tributed
r.h.s.

    right hand side
l.h.s.

    left hand side
w.r.t.

    with respect to
MSE

    mean squared error
MMSE

    minimum mean squared error
SNR

    signal-to-noise ratio
SINR

    signal to interference and noise ratio
MF

    match filtering

##  I Introduction

Consider the following classical problem which raises in several sensing
scenarios: A set of signal samples $x_{1},\ldots,x_{N}$ are to be sensed in a
noisy environment via $K$ sensors within a limited time frame. Each sensor is
tunable and can observe various linear combinations of signal samples. The
ultimate goal is to collect some observations from which the signal samples
are recovered with minimum distortion. There are in general two approaches to
address this goal:

  1. 1.

The undersampled non-adaptive approach, in which the sensors are tuned once at
the beginning of the time frame and kept fixed for the whole duration. In this
case, $K$ high quality observations, i.e., observations with high signal-to-
noise ratio (SNR), are collected. These are then given to an estimator to
recover the signal samples.

  2. 2.

The oversampled adaptive sensing (OAS) framework [1, 2, 3], in which the given
time frame is divided into multiple subframes. In each subframe, $K$
observations are collected. These observations, along with those collected in
previous subframes, are used to decide for a new tuning strategy in the next
subframe. Clearly in this case, a larger number of observations are collected.
This increase in the number of observations is obtained at the expense of
lower quality for each particular observation111Note that since subframes are
shorter than the whole time frame, the noise variance is higher. This results
in lower SNR for a particular observation. This point is illustrated in the
Section II..

Despite initial studies on adaptive approaches for signal sensing, e.g., [4,
5, 6], it was commonly believed in the literature that adaptation does not
result in a significant performance enhancement, assuming that the trade-off
between the quality and quantity of the observations leads to no performance
gain. In [1], this myth was shown to be wrong. OAS illustrates this trade-off
between the number of observations and their quality as follows: For cases
with some prior information on the signal samples, e.g., sparse recovery, the
adaptation can result on a significant enhancement. Nevertheless, in lack of
prior information, the conventional belief holds.

The superiority of OAS framework is intuitively understood as follows: In
presence of prior information, several samples222For instance, zero samples in
the example of sparse recovery. are effectively recovered in the primary
subframes. By omitting those samples, the space of unknowns shrinks, and hence
the recovery can be performed effectively.

###  I-A Contributions

The initial studies of OAS consider no restrictions on the sensor tuning. This
is an idealistic assumption which is not necessarily satisfied in practice;
see for instance the example given in Section I-B. It is hence more realistic
to assume that the tuning of sensors is controlled by a predefined codebook
which includes available choices of tuning. Given this new restriction, this
study develops an OAS algorithm which collects observations in each subframe
using a predefined codebook. The performance of the proposed algorithm is
investigated in various respects, and the impacts of limited codebook size are
clarified through several experiments.

###  I-B Motivation and Applications

To give an intuition on sensing via a predefined codebook, let us consider the
following example: A signal is sampled via a preinstalled sensor network
within a limited time frame. At a given time, only a subset of sensors are
active in the network. To apply OAS, the time frame is divided into several
subframes, and the active sensors are altered in each subframe. Clearly, the
sensing procedure in this setting is restricted by the preinstalled network
and cannot be changed arbitrarily.

Similar to the given example, in several practical scenarios, sensing is
performed by means of a preinstalled setting whose possible configurations for
observing a signal is restricted by a codebook. Examples of such scenarios are
found in radar and positioning systems. In these applications, it is often the
case that prior information is available on the target signal; for instance,
in many radar systems, the observed signal is known to be sparse. The main
motivation of this study is to extend the scope of OAS to these applications.

###  I-C Notation

Scalars, vectors and matrices are shown with non-bold, bold lower case and
bold upper case letters, respectively. $\mathbf{I}_{K}$ and
$\boldsymbol{0}_{K\times N}$ are the $K\times K$ identity matrix and $K\times
N$ all-zero matrix, respectively. $\mathbf{A}^{\mathsf{T}}$ denotes the
transpose of $\mathbf{A}$. The set of real numbers is shown by $\mathbb{R}$.
We use the shortened notation $[N]$ to represent
$\left\\{1,\ldots,N\right\\}$.

##  II Problem Formulation

Consider ${\boldsymbol{x}}\in\mathbb{R}^{N}$ containing $N$ signal samples. We
postulate333This postulation is not necessarily matched to the true
distribution. that the samples are independent and identically dis-tributed
(i.i.d.) with prior distribution $q\left(x\right)$. These samples are observed
via $K$ sensors within a restricted time interval of duration $T$. A
particular sensor $k\in\left[K\right]$ observes a noisy linear combination of
signal samples, i.e.,

$\displaystyle
y_{k}=\left.{\mathbf{a}}_{k}^{\mathsf{T}}\right.{\boldsymbol{x}}+z_{k}$ (1)

for additive noise $z_{k}$ and the vector of coefficients
${\mathbf{a}}_{k}\in\mathbb{R}^{N}$:

  1. 1.

The noise power is reversely proportional to the observation time: Assume that
sensor $k$ operates for a time frame of length $t\leq T$. Then, $z_{k}$ is
assumed to be zero-mean Gaussian with variance
$\sigma^{2}\left(t\right)={\sigma^{2}_{0}}/{t}$ for some $\sigma^{2}_{0}\geq
0$.

  2. 2.

The coefficient vectors ${\mathbf{a}}_{k}$ for $k\in\left[K\right]$ are
selected from a predefined codebook
$\mathcal{C}=\left\\{{\mathbf{c}}_{1},\ldots,{\mathbf{c}}_{S}\right\\}$ where
$S\geq K$ and ${\mathbf{c}}_{s}\in\mathbb{R}^{N}$ for $s\in\left[S\right]$.
This means that

$\displaystyle\left\\{{\mathbf{a}}_{1},\ldots,{\mathbf{a}}_{K}\right\\}\subseteq\mathcal{C}=\left\\{{\mathbf{c}}_{1},\ldots,{\mathbf{c}}_{S}\right\\}.$
(2)

There are two key points in the sensing model (1) which deviate from the
classical models:

  1. 1.

The sensing quality is time-dependent.

  2. 2.

The sensing is performed using a predefined codebook.

The time-dependent model follows the fact that the SNR of a particular
observation linearly scales with the duration of observation444See [1] for
detailed illustrations.. This is a typical assumption in OAS. The latter point
is the key difference of the considered system model, to the typical settings
considered for OAS.

##  III OAS Framework

In a nutshell, the OAS framework can be represented via the following steps:

  1. (a)

The time frame is divided into $M$ subframes.

  2. (b)

In subframe $m\in\left[M\right]$, the sensors observe

$\displaystyle{\boldsymbol{y}}_{m}=\left.\mathbf{A}_{m}\right.{\boldsymbol{x}}+{\boldsymbol{z}}_{m},$
(3)

where ${\boldsymbol{z}}_{m}\sim\mathcal{N}\left(\boldsymbol{0},\sigma^{2}_{\rm
s}\mathbf{I}_{K}\right)$ with $\sigma^{2}_{\rm s}=M\sigma^{2}\left(T\right)$,
and

$\displaystyle\mathbf{A}_{m}=\left[{\mathbf{a}}_{1}\left(m\right),\ldots,{\mathbf{a}}_{K}\left(m\right)\right]^{\mathsf{T}}$
(4)

with
$\left\\{{\mathbf{a}}_{1}\left(m\right),\ldots,{\mathbf{a}}_{K}\left(m\right)\right\\}\subseteq\mathcal{C}$
containing the vectors selected in subframe $m$.

  3. (c)

A processing unit collects the observations up to subframe $m$ in a matrix of
stacked observations

$\displaystyle\mathbf{Y}_{m}\coloneqq\left[{\boldsymbol{y}}_{1},\ldots,{\boldsymbol{y}}_{m}\right].$
(5)

It then uses a Bayesian estimator to calculate some soft information. In
general, the soft information is written as

$\displaystyle\hat{{\boldsymbol{x}}}_{m}=\mathbb{E}\hskip
1.42262pt\left\\{{\boldsymbol{x}}|\mathbf{Y}_{m},\mathbb{A}_{m}\right\\},$ (6)

where $\mathbb{A}_{m}$ denotes the collection of sensing matrices up to
subframe $m$, i.e.,

$\displaystyle\mathbb{A}_{m}=\left\\{\mathbf{A}_{1},\ldots,\mathbf{A}_{m}\right\\},$
(7)

and the expectation is taken with respect to the posterior distribution
$\mathrm{p}\left({\boldsymbol{x}}|\mathbf{Y}_{m},\mathbb{A}_{m}\right)$ which
is calculated from the likelihood using the postulated prior distribution
$q\left(x\right)$.

  4. (d)

The processing unit specifies the sensing matrix of subframe $m+1$ using the
soft information. This means that

$\displaystyle\mathbf{A}_{m+1}=f_{\rm
Adp}\left(\hat{{\boldsymbol{x}}}_{m}\right)$ (8)

for some adaptation function $f_{\rm Adp}\left(\cdot\right)$.

###  III-A Performance Characterization

The conventional metric which quantifies the recovery performance is the
average distortion: Let
$\Delta\left(\cdot;\cdot\right):\mathbb{R}\times\mathbb{R}\mapsto\mathbb{R}_{0}^{+}$
be a distortion function555For instance, the Euclidean distance.. The
distortion between two vectors
${\boldsymbol{x}},\hat{\boldsymbol{x}}\in\mathbb{R}^{N}$ with respect to
$\Delta\left(\cdot;\cdot\right)$ is determined as

$\displaystyle\Delta\left({\boldsymbol{x}};\hat{\boldsymbol{x}}\right)=\sum_{n=1}^{N}\Delta\left(x_{n};\hat{x}_{n}\right).$
(9)

In an OAS-based algorithm, the signal samples are finally recovered as
${{\boldsymbol{x}}^{\star}}=\mathrm{g}\left(\hat{{\boldsymbol{x}}}_{M}\right)$
where $\mathrm{g}\left(\cdot\right)$ is a decision function666For example, a
hard or soft thresholding function. used to recover the signal samples from
the soft estimate $\hat{{\boldsymbol{x}}}_{M}$. Consequently, the average
distortion is

$\displaystyle D=\mathbb{E}\hskip
1.42262pt\left\\{\Delta\left({\boldsymbol{x}};{\boldsymbol{x}}^{\star}\right)\right\\}.$
(10)

A non-adaptive recovery technique can be seen as an OAS-based algorithm with a
single time subframe. The performance in this case is hence characterized by
setting $M=1$.

##  IV OAS via a Predefined Codebook

In settings with no predefined codebook, designing OAS-based algorithms deals
with high degrees of freedom. In fact, sensing matrices are freely constructed
in such scenarios. With a predefined codebook, algorithm design is
significantly restricted, since degrees of freedom are limited via the
codebook.

In the sequel, we discuss the design via a sample algorithm which performs OAS
via a predefined codebook. The algorithm is given in Algorithm 1. In this
algorithm,

  * •

The selector operator $\mathrm{Sel}\left(\mathbb{X}\right)$, for an index set
$\mathbb{X}\subseteq\left[N\right]$ with $L$ indices, returns an $L\times N$
matrix. Setting $\mathbf{P}=\mathrm{Sel}\left(\mathbb{X}\right)$, the
$\ell$-th row of $\mathbf{P}$ is the standard basis of $\mathbb{R}^{N}$ whose
non-zero entry occurs at the $\ell$-th index of $\mathbb{X}$, e.g., for $N=4$
and $\mathbb{X}=\left\\{1,3\right\\}$, we have

$\displaystyle\mathrm{Sel}\left(\mathbb{X}\right)=\begin{bmatrix}1&0&0&0\\\
0&0&1&0\end{bmatrix}.$ (11)

  * •

The posterior distribution $\mathrm{p}_{0}\left(u|y,\sigma^{2}\right)$ is a
scalar distribution and is given by

$\displaystyle\mathrm{p}_{0}\left(u|y,\sigma^{2}\right)=\frac{\displaystyle\exp\left\\{-\frac{\left(y-u\right)^{2}}{2\sigma^{2}}\right\\}q\left(u\right)}{\displaystyle\int\exp\left\\{-\frac{\left(y-u\right)^{2}}{2\sigma^{2}}\right\\}q\left(u\right)\mathrm{d}u}.$

  * •

$\mathcal{A}\left(\mathcal{C},K,\mathbb{F},\hat{{\boldsymbol{x}}}\right)$ is a
selection algorithm selecting $K$ vectors from the codebook $\mathcal{C}$ in
order to sense signal samples which are indexed by $\mathbb{F}$ using the
current estimate $\hat{{\boldsymbol{x}}}$.

Algorithm 1 OAS via a Predefined Codebook

$K$, $L$, codebook $\mathcal{C}$ and postulated prior $q(x_{n})$.

For $n\in\left[N\right]$, set $d_{n}=+\infty$,
$\hat{x}_{n}=\mathrm{y}_{n}=\hat{\sigma}_{n}^{2}=0$ and
$\mathbb{M}_{n}\left(0\right)=\emptyset$.

for $m\in\left[M\right]$ do

  1. 1.

Sort $d_{1},\ldots d_{N}$ as $d_{i_{1}^{m}},\ldots d_{i_{N}^{m}}$ such that

$\displaystyle d_{i_{1}^{m}}\geq\ldots\geq d_{i_{N}^{m}}.$

  2. 2.

Set $\mathbb{F}_{m}$ by worst-case adaptation as

$\displaystyle\mathbb{F}_{m}=\left\\{{i_{1}^{m}},{i_{L}^{m}}\right\\},$

and update
$\mathbb{M}_{n}\left(m\right)=\mathbb{M}_{n}\left(m-1\right)\cup\left\\{m\right\\}$
for $n\in\mathbb{F}_{m}$.

  3. 3.

Set
$\mathbf{A}_{m}=\left[{\mathbf{a}}_{1}\left(m\right),\ldots,{\mathbf{a}}_{K}\left(m\right)\right]^{\mathsf{T}}$
where

$\displaystyle\left\\{{\mathbf{a}}_{1}\left(m\right),\ldots,{\mathbf{a}}_{K}\left(m\right)\right\\}=\mathcal{A}\left(\mathcal{C},K,\mathbb{F}_{m},\hat{{\boldsymbol{x}}}\right)$

  4. 4.

Sense the samples for duration $T/M$ via $\mathbf{A}_{m}$.

  5. 5.

Set $\mathbf{P}_{m}=\mathrm{Sel}\left(\mathbb{F}_{m}\right)$ and
$\mathbf{E}_{m}=\mathrm{Sel}\left(\left[N\right]\backslash\mathbb{F}_{m}\right)$.
Let $\tilde{{\boldsymbol{x}}}_{m}=\mathbf{E}_{m}\hat{{\boldsymbol{x}}}$,
$\mathbf{Q}_{m}=\mathbf{A}_{m}\mathbf{P}_{m}^{\mathsf{T}}$ and
$\mathbf{W}_{m}=\mathbf{A}_{m}\mathbf{E}_{m}^{\mathsf{T}}$, and set

$\displaystyle\mathbf{F}_{m}$
$\displaystyle=\left(\mathbf{Q}_{m}^{\mathsf{T}}\mathbf{Q}_{m}\right)^{-1}\mathbf{Q}_{m}^{\mathsf{T}}$

  6. 6.

Decouple the observations as

$\displaystyle\hat{{\boldsymbol{y}}}_{m}=\mathbf{F}_{m}\left({\boldsymbol{y}}_{m}-\mathbf{W}_{m}\tilde{{\boldsymbol{x}}}_{m}\right).$

  7. 7.

For $\ell\in\left[L\right]$, update

$\displaystyle\mathrm{y}_{i_{\ell}^{m}}$
$\displaystyle=\mathrm{y}_{i_{\ell}^{m}}+\hat{y}_{m,\ell}$
$\displaystyle\hat{\sigma}_{i_{\ell}^{m}}^{2}$
$\displaystyle=\hat{\sigma}_{i_{\ell}^{m}}^{2}+M\lVert\mathbf{f}_{m,\ell}\rVert^{2}\sigma^{2}$

with $\mathbf{f}_{m,\ell}^{\mathsf{T}}$ being the $\ell$-th row of
$\mathbf{F}_{m}$.

  8. 8.

For $n\in\mathbb{F}_{m}$, update

$\displaystyle\hat{x}_{n}$ $\displaystyle=\int
u\;\mathrm{p}_{0}\left(u\left|\frac{\mathrm{y}_{n}}{\lvert\mathbb{M}_{n}\left(m\right)\rvert},\frac{\hat{\sigma}_{n}^{2}}{\lvert\mathbb{M}_{n}\left(m\right)\rvert^{2}}\right.\right)\mathrm{d}u$
$\displaystyle{d}_{n}$
$\displaystyle=\int\left(u-\hat{x}_{n}\right)^{2}\mathrm{p}_{0}\left(u\left|\frac{\mathrm{y}_{n}}{\lvert\mathbb{M}_{n}\left(m\right)\rvert},\frac{\hat{\sigma}_{n}^{2}}{\lvert\mathbb{M}_{n}\left(m\right)\rvert^{2}}\right.\right)\mathrm{d}u$

end for

###  IV-A Derivation of Algorithm 1

In a nutshell, Algorithm 1 follows these steps: It first focuses on a subset
of $L\leq K$ signal samples via worst-case adaptation and selects $K$ vectors
from the codebook to sense them. Using the fact that recovery of this subset
is a determined problem, the algorithm uses linear filtering to decouple the
observations. It then cancels out the other $N-L$ signal samples using the
estimates obtained in previous subframes and estimates the $L$ samples via a
scalar Bayesian estimator assuming that the residual terms are Gaussian.

To derive the update rules given in the algorithm, consider subframe $m-1$ at
which the signal samples are estimated as $\hat{x}_{1},\ldots,\hat{x}_{N}$ and
posterior distortions $d_{1},\ldots,d_{N}$ are calculated for the signal
samples. The posterior distortion $d_{n}$ for the $n$-th signal sample is
determined using $y_{n}$, which is a statistic on $x_{n}$, and estimated noise
variance $\hat{\sigma}^{2}_{n}$ as

$\displaystyle{d}_{n}$
$\displaystyle=\int\left(u-\hat{x}_{n}\right)^{2}\mathrm{p}_{0}\left(u|y_{n},\hat{\sigma}^{2}_{n}\right)\mathrm{d}u.$
(12)

The statistic and estimated noise variance are calculated from
$\mathbf{Y}_{m-1}$ and $\mathbb{A}_{m-1}$ in a way that will be clarified in
the sequel.

For $L\leq K$, the algorithm finds $L$ samples whose posterior distortions are
maximum. The indices of these samples are collected in
$\mathbb{F}_{m}\subseteq\left[N\right]$. A selection algorithm
$\mathcal{A}\left(\mathcal{C},K,\mathbb{F}_{m},\hat{{\boldsymbol{x}}}\right)$
is then used to select $K$ vectors from the codebook $\mathcal{C}$. At this
point, we consider a generic selection algorithm, and leave the discussions on
the design of
$\mathcal{A}\left(\mathcal{C},K,\mathbb{F}_{m},\hat{{\boldsymbol{x}}}\right)$
for later.

Let $\mathbf{A}_{m}$ be the collection of vectors selected in subframe $m$.
The observations are hence given by (3) which read

$\displaystyle{\boldsymbol{y}}_{m}$
$\displaystyle=\mathbf{A}_{m}{\boldsymbol{x}}+{\boldsymbol{z}}_{m}$ (13a)
$\displaystyle=\left(\mathbf{A}_{m}\mathbf{P}_{m}^{\mathsf{T}}\mathbf{P}_{m}+\mathbf{A}_{m}\mathbf{E}_{m}^{\mathsf{T}}\mathbf{E}_{m}\right){\boldsymbol{x}}+{\boldsymbol{z}}_{m}$
(13b)

where $\mathbf{P}_{m}=\mathrm{Sel}\left(\mathbb{F}_{m}\right)$ and
$\mathbf{E}_{m}=\mathrm{Sel}\left(\left[N\right]\backslash\mathbb{F}_{m}\right)$.
Defining the matrices
$\mathbf{Q}_{m}=\mathbf{A}_{m}\mathbf{P}_{m}^{\mathsf{T}}$ and
$\mathbf{W}_{m}=\mathbf{A}_{m}\mathbf{E}_{m}^{\mathsf{T}}$, we can write

$\displaystyle{\boldsymbol{y}}_{m}=\mathbf{Q}_{m}\breve{{\boldsymbol{x}}}_{m}+\mathbf{W}_{m}\bar{{\boldsymbol{x}}}_{m}+{\boldsymbol{z}}_{m}$
(14)

where $\breve{{\boldsymbol{x}}}_{m}=\mathbf{P}_{m}{\boldsymbol{x}}$ contains
the $L$ samples selected by $\mathbb{F}_{m}$ and
$\bar{{\boldsymbol{x}}}_{m}=\mathbf{E}_{m}{\boldsymbol{x}}$ consists of the
remaining $N-L$ signal samples. Since the entries of
$\bar{{\boldsymbol{x}}}_{m}$ are not estimated in this subframe, we
approximate them by the estimates of the previous subframe, i.e.,
$\bar{{\boldsymbol{x}}}_{m}\approx\tilde{{\boldsymbol{x}}}_{m}=\mathbf{E}_{m}\hat{{\boldsymbol{x}}}$.
We hence can write

$\displaystyle{\boldsymbol{y}}_{m}-\mathbf{W}_{m}\tilde{{\boldsymbol{x}}}_{m}\approx\mathbf{Q}_{m}\breve{{\boldsymbol{x}}}_{m}+{\boldsymbol{z}}_{m}.$
(15)

Noting that $\mathbf{Q}_{m}\in\mathbb{R}^{K\times L}$ with $K\geq L$, we could
conclude that for a proper selection of codewords
$\mathbf{Q}_{m}^{\mathsf{T}}\mathbf{Q}_{m}$ is full-rank, and thus, we can
calculate the pseudo-inverse of $\mathbf{Q}_{m}$ as

$\displaystyle\mathbf{F}_{m}$
$\displaystyle=\left(\mathbf{Q}_{m}^{\mathsf{T}}\mathbf{Q}_{m}\right)^{-1}\mathbf{Q}_{m}^{\mathsf{T}}.$
(16)

Consequently, the observations can be decoupled into $L$ samples with additive
noise as

$\displaystyle\hat{{\boldsymbol{y}}}_{m}$
$\displaystyle=\mathbf{F}_{m}\left({\boldsymbol{y}}_{m}-\mathbf{W}_{m}\tilde{{\boldsymbol{x}}}_{m}\right)$
(17)
$\displaystyle\approx\breve{{\boldsymbol{x}}}_{m}+\mathbf{F}_{m}{\boldsymbol{z}}_{m}.$
(18)

We now consider this mismatched assumption that
$\mathbf{F}_{m}{\boldsymbol{z}}_{m}$ is a Gaussian vector with independent
entries777Such an assumption is asymptotically correct for large random
codebooks with i.i.d. entries.. As the result, the $\ell$-th entry of
$\hat{{\boldsymbol{y}}}_{m}$ can be written as

$\displaystyle\hat{y}_{m,\ell}\approx x_{i_{\ell}^{m}}+\hat{z}_{m,\ell}$ (19)

where $i_{\ell}^{m}$ denotes the $\ell$-th entry of $\mathbb{F}_{m}$, and
decoupled noise $\hat{z}_{m,\ell}$ is zero-mean Gaussian with variance
$M\lVert\mathbf{f}_{m,\ell}\rVert^{2}\sigma^{2}$. Here,
$\mathbf{f}_{m,\ell}^{\mathsf{T}}$ denotes the $\ell$-th row of
$\mathbf{F}_{m}$.

We now use $\hat{y}_{m,\ell}$ to update the statistic on $x_{i_{\ell}^{m}}$ by
adding the decoupled observation corresponding to $x_{i_{\ell}^{m}}$ to the
current statistic: Let $\mathrm{y}_{i_{\ell}^{m}}\left(m-1\right)$ denote the
statistic on $x_{i_{\ell}^{m}}$ collected in subframes $1,\ldots,m-1$. We
update this parameter as

$\displaystyle\mathrm{y}_{i_{\ell}^{m}}\left(m\right)$
$\displaystyle=\mathrm{y}_{i_{\ell}^{m}}\left(m-1\right)+\hat{y}_{m,\ell}.$
(20)

Let $\mathbb{M}_{n}\left(m\right)\subseteq\left[m\right]$ collects all
subframes at which $n\in\mathbb{F}_{m}$, i.e., $x_{n}$ is selected to be
sensed in these subframes. Moreover, assume that
$\hat{\sigma}_{n}^{2}\left(m\right)$ denotes the sum of all decoupled noise
variances in $\mathbb{M}_{n}\left(m\right)$ which correspond to index $n$.
After updating our statistic on $x_{i_{\ell}^{m}}$, we have

$\displaystyle\mathrm{y}_{i_{\ell}^{m}}\left(m\right)\approx\lvert\mathbb{M}_{n}\left(m\right)\rvert
x_{i_{\ell}^{m}}+\mathrm{z}_{i_{\ell}^{m}}\left(m\right)$ (21)

where $\mathrm{z}_{i_{\ell}^{m}}\left(m\right)$ is given by

$\displaystyle\mathrm{z}_{i_{\ell}^{m}}\left(m\right)=\mathrm{z}_{i_{\ell}^{m}}\left(m-1\right)+\hat{z}_{i_{\ell}^{m}}\left(m\right).$
(22)

Hence, the estimated noise variance is updated as

$\displaystyle\hat{\sigma}_{i_{\ell}^{m}}^{2}\left(m\right)=\hat{\sigma}_{i_{\ell}^{m}}^{2}\left(m-1\right)+M\lVert\mathbf{f}_{m,\ell}\rVert^{2}\sigma^{2}.$
(23)

Using the Bayesian estimator, we finally update the estimate $\hat{x}_{n}$ and
the posterior distortion $d_{n}$ via the updated statistics $\mathrm{y}_{n}$
for $n\in\mathbb{F}_{m}$. This concludes Algorithm 1.

##  V Example of Sparse Recovery

Initial investigations demonstrate that OAS can significantly outperform
classical sparse recovery approaches. This is due to the prior information on
the sparsity of signal samples. In this section, we employ our algorithm to
address this example.

For a sparse signal, the prior distribution of the samples is

$\displaystyle
q\left(x\right)=\left(1-\rho\right)\delta\left(x\right)+\rho\phi\left(x\right)$
(24)

for sparsity factor $\rho\in\left[0,1\right]$ and a distribution
$\phi\left(x\right)$ satisfying

$\displaystyle\int_{0^{-}}^{0^{+}}\phi\left(x\right)=0.$ (25)

The model indicates that samples are zero with probability $1-\rho$ and are
drawn from distribution $\phi\left(x\right)$ with probability $\rho$. In the
sequel, we consider an i.i.d. sparse-Gaussian signal for which the samples are
i.i.d. and

$\displaystyle\phi\left(x\right)=\frac{1}{\sqrt{2\pi}}\exp\left\\{-\frac{x^{2}}{2}\right\\}.$
(26)

Considering the system model, the main task in this example is to recover the
samples of a given sparse signal using the observations made within a time
frame of length $T$ via the $K$ active sensors. This task can be addressed in
two ways:

  * •

Invoking the theory of compressive sensing, one can select a sensing matrix
from the codebook, and collect $K$ observations. A classical sparse recovery
algorithm, e.g., LASSO, is then used to recover the signal samples.

  * •

An alternative approach is to use an OAS algorithm, e.g., Algorithm 1, to
adaptively estimate the signal samples.

We investigate these two approaches in the sequel.

###  V-A Sparse Recovery via OAS

For a sparse-Gaussian prior distribution, the Bayesian estimation in Algorithm
1 reduces to

$\displaystyle\hat{x}_{n}$
$\displaystyle=G\left(\frac{\mathrm{y}_{n}}{\lvert\mathbb{M}_{n}\left(m\right)\rvert}\left|\frac{\hat{\sigma}_{n}^{2}}{\lvert\mathbb{M}_{n}\left(m\right)\rvert^{2}}\right.\right),$
$\displaystyle d_{n}$
$\displaystyle=E\left(\frac{\mathrm{y}_{n}}{\lvert\mathbb{M}_{n}\left(m\right)\rvert}\left|\frac{\hat{\sigma}_{n}^{2}}{\lvert\mathbb{M}_{n}\left(m\right)\rvert^{2}}\right.\right).$

where $G\left(y|\sigma^{2}\right)$ and $E\left(y|\sigma^{2}\right)$ are given
by

$\displaystyle G\left(y|\sigma^{2}\right)$
$\displaystyle=\frac{I_{1}\left(y|\sigma^{2}\right)y}{J\left(y|\sigma^{2}\right)}$
(27a) $\displaystyle E\left(y|\sigma^{2}\right)$
$\displaystyle=\frac{I_{1}\left(y|\sigma^{2}\right)}{J\left(y|\sigma^{2}\right)}\left(\sigma^{2}+\frac{I_{0}\left(y|\sigma^{2}\right)}{J\left(y|\sigma^{2}\right)}y^{2}\right)$
(27b)

with $I_{0}\left(y|\sigma^{2}\right)$, $I_{1}\left(y|\sigma^{2}\right)$ and
$J\left(y|\sigma^{2}\right)$ being

$\displaystyle I_{0}\left(y|\sigma^{2}\right)$
$\displaystyle=\left(1-\rho\right)\sqrt{1+\sigma^{2}}\exp\left\\{-\frac{y^{2}}{2{\sigma^{2}}}\right\\}$
(28a) $\displaystyle I_{1}\left(y|\sigma^{2}\right)$
$\displaystyle=\rho\;\sigma\exp\left\\{-\frac{y^{2}}{2\left(1+\sigma^{2}\right)}\right\\}$
(28b) $\displaystyle J\left(y|\sigma^{2}\right)$
$\displaystyle=\left(1+\sigma^{2}\right)\left[I_{0}\left(y|\sigma^{2}\right)+I_{1}\left(y|\sigma^{2}\right)\right].$
(28c)

##  VI Numerical Experiments

$50$$100$$150$$200$$-25$$-20$$-15$$-10$(a)$L$MSE in
[dB]OASLASSOMMSE$20$$40$$60$$80$$100$$-25$$-20$$-15$$-10$(b)$L$MSE in
[dB]$10$$20$$30$$40$$50$$-20$$-15$$-10$(c)$L$MSE in
[dB]$10$$20$$30$$40$$-20$$-15$$-10$(d)$L$MSE in [dB] Figure 1: MSE versus $L$
for various compression rates: (a) $R=1$, (b) $R=2$, (c) $R=4$, and (d) $R=5$.
Here, the time frame is divided to $M=80$ subframes and $S=1000$.

For numerical investigations, we consider a sparse-Gaussian signal with
$N=200$ samples whose sparsity factor is $\rho=0.1$. The noise variance for an
interval of length $T$ is $\sigma^{2}=0.01$. This means that by dividing the
time frame into $M$ subframes, the noise variance in each subframe is $0.01M$.
The compression rate is further defined as $R={N}/{K}$.

The codebook is generated randomly from an i.i.d. Gaussian ensemble. This
means that the vectors in $\mathcal{C}$ are independent and have i.i.d. zero-
mean Gaussian entries with variance $1/\sqrt{K}$.

As benchmarks, we consider a classic compressive sensing setting in which the
$K$ observations are collected using $K$ randomly selected vectors. We
evaluate the performance for two recovery algorithms:  1. LASSOin which the
samples are recovered via regularized $\ell_{1}$-norm minimization, and  2.
mmseestimator in which the samples are recovered via the optimal Bayesian
estimator.  It is worth mentioning that the computational complexity of the
former algorithm is moderate while the latter approach is not numerically
tractable. We hence use the asymptotic characterizations for these schemes [7,
8].

The performance is evaluated in terms of the average mean squared error (MSE)
given by the average distortion when

$\displaystyle\Delta\left(x_{n};\hat{x}_{n}\right)=\frac{1}{N}\lvert
x_{n}-\hat{x}_{n}\rvert^{2}.$ (29)

The first experiment considers a codebook of size $S=1000$ and compares the
classical approaches with an OAS scheme in which $M=80$. The results are shown
in Fig. 1 for multiple compression rates. The OAS scheme uses Algorithm 1 with
a random selection algorithm, i.e., the $K$ vectors are selected randomly in
each subframe. As the result shows, the OAS-based algorithm performs close to
LASSO in low compression rates, while at $R=5$, it even outperforms the MMSE
bound.

The behavior illustrated via Fig. 1 is consistent with the prior
investigations on OAS with no predefined codebooks [1, 2, 3]. In fact, those
initial cases can be observed as a special case of OAS via a predefined
codebook whose size goes to infinity. To investigate the performance
degradation imposed by the codebook restriction, we further plot the MSE
achieved via the OAS-based algorithm against the codebook size in Fig. 2. In
this figure, we set $R=4$ and $M=60$. $L$ is further set to $L=K/2=25$. As the
figure shows, the MSE drops fast and converges to the asymptotic value which
specifies the OAS performance with no predefined codebooks.

$200$$400$$600$$800$$-20$$-15$$-10$$-5$$S$MSE in [dB]OASLASSOMMSE Figure 2:
MSE against the codebook size. Here, the compression rate is set to $R=4$,
number of subframes $M=60$ and $L=K/2=25$.

Similar behavior is observed in Fig. 3 in which the MSE is plotted against the
number of subframes for two different codebook sizes. Here, we set $R=4$ and
$L=K/2=25$. From the figure, it is observed that the OAS-algorithm converges
to its asymptotic relatively fast, and further increase in the number of
subframes improves the performance negligibly.

$20$$40$$60$$80$$100$$120$$-20$$-15$$-10$$-5$$M$MSE in
[dB]$S=500$$S=1000$LASSOMMSE Figure 3: MSE against the number of subframes for
different codebook sizes. Here, we set $R=4$ and $L=K/2=25$.

###  VI-A Impact of Selection Algorithm

Figs. 1-3 consider a random selection algorithm. The performance of the OAS
algorithm can be further improved by developing a more efficient selection
algorithm. Such an algorithm can be developed by defining a concept of
optimality. In the sequel, we give an instance of such algorithms.

From derivations in Section IV-A, we know that Algorithm 1 cancels out the
residual samples in each subframe using the estimates given in the previous
subframes; see (15). This means that the algorithm relies on this postulation
that estimated samples converge to good estimates as the algorithm evolves.
Nevertheless, the approximated approach leads to unwanted interference terms
in each subframe.

The interference in the subframe $m$ is ideally suppressed, if the residual
samples lie in the kernel of $\mathbf{W}_{m}$, i.e.
$\mathbf{W}_{m}\bar{{\boldsymbol{x}}}_{m}=0$. However, this constraint is not
necessarily fulfilled, since  1. the residual samples are unknown, and  2.
$\mathbf{A}_{m}$is constructed from the codebook.  The first issue is
addressed by following the approximation used in Algorithm 1, i.e., replacing
$\bar{{\boldsymbol{x}}}_{m}$ with $\tilde{{\boldsymbol{x}}}_{m}$. For the
second issue, an optimal approach is to search the codebook, such that the
sum-power of the interference terms is minimized, i.e., setting
$\mathbf{A}_{m}=\left[{}_{1}\left(m\right),\ldots,_{K}\left(m\right)\right]^{\mathsf{T}}$
where

${}_{1}\left(m\right),\ldots,_{K}\left(m\right)=\operatorname*{argmin}_{{\mathbf{u}}_{1},\ldots,{\mathbf{u}}_{K}\in\mathcal{C}}\lVert\mathbf{U}\mathbf{E}_{m}^{\mathsf{T}}\tilde{{\boldsymbol{x}}}_{m}\rVert^{2}$
(30)

with
$\mathbf{U}=\left[{\mathbf{u}}_{1},\ldots,{\mathbf{u}}_{K}\right]^{\mathsf{T}}$.

The selection algorithm in (30) deals with an exhaustive search whose number
of choices exponentially grows by the codebook size. For practical scenarios
with large codebooks, this selection algorithm is not numerically tractable.
One can hence approximate the solution with a greedy algorithm. An example is
Algorithm 2 which uses stepwise regression.

Algorithm 2 Adaptation Algorithm
$\mathcal{A}\left(\mathcal{C},K,\mathbb{F},\hat{{\boldsymbol{x}}}\right)$

$K$, codebook $\mathcal{C}$ and index set $\mathbb{F}\subset\left[N\right]$.

Set ${}_{1}={\mathbf{c}}_{i_{1}}$ with ${\mathbf{c}}_{i_{1}}$ being selected
at random from $\mathcal{C}$. Let
$\mathbb{I}_{1}=\left[S\right]\backslash\left\\{i_{1}\right\\}$,
$\mathbf{P}=\mathrm{Sel}\left(\mathbb{F}\right)$,
$\mathbf{E}=\mathrm{Sel}\left(\left[N\right]\backslash\mathbb{F}\right)$ and
$\tilde{{\boldsymbol{x}}}=\mathbf{E}\hat{{\boldsymbol{x}}}$.

for $k\in\left[2:K\right]$ do

  1. 1.

Set $\mathbf{A}_{k}=\left[{\mathbf{a}}_{1},\ldots,{\mathbf{a}}_{k-1}\right]$
and let

$\displaystyle\mathbf{v}_{k}=\mathbf{P}\mathbf{A}_{k}\mathbf{A}_{k}^{\mathsf{T}}\mathbf{E}^{\mathsf{T}}\tilde{{\boldsymbol{x}}}$

  2. 2.

For ${\mathbf{u}}\in\mathbb{R}^{N}$, let

$\displaystyle
f_{k}\left({\mathbf{u}}\right)=\lVert\left({\mathbf{u}}^{\mathsf{T}}\mathbf{E}^{\mathsf{T}}\tilde{{\boldsymbol{x}}}\right)\mathbf{P}{\mathbf{u}}+\mathbf{v}_{k}\rVert^{2}.$

  3. 3.

Set ${}_{k}={\mathbf{c}}_{i_{k}}$, where

$\displaystyle
i_{k}=\operatorname*{argmin}_{i\in\mathbb{I}_{k-1}}f_{k}\left({\mathbf{c}}_{i}\right).$

  4. 4.

Update $\mathbb{I}_{k}=\mathbb{I}_{k-1}\backslash\left\\{i_{k}\right\\}$.

end for

Fig. 4 compares the performance of Algorithm 2 with random selection. Here,
Algorithm 1 is run with both of the selection algorithms when $R=4$ and a
codebook of size $S=500$ is available. The results are plotted for $M=20$
subframes. As the figure shows, the greedy algorithm outperforms the random
selection which agrees with the intuition.

$10$$20$$30$$40$$50$$-20$$-15$$-10$$-5$$L$MSE in [dB]RandomStepwiseLASSOMMSE
Figure 4: MSE versus $L$ for random and stepwise adaptations. Here, the
compression rate is set to $R=4$, codebook size $S=500$ and number of
subframes is $M=20$.

##  VII Conclusions

In presence of a predefined codebook, design of OAS-based algorithms deal with
more challenges. This is due to the lower degrees of freedom provided by the
setting. Although this restriction leads to degraded performance, the same
behavior as the one, observed in the basic form of OAS framework, is reported.
The performance degradation is further compensated by effective design of
selection algorithms and higher number of subframes. These are feasible means
in many applications in which OAS seems to be a good candidate.

## References

  * [1] R. R. Müller, A. Bereyhi, and C. Mecklenbräuker, “Oversampled adaptive sensing,” _Proc. Inf. Th. & App. Work. (ITA)_, pp. 1–7, Feb. 2018, USA.
  * [2] R. R. Müller, A. Bereyhi, and C. F. Mecklenbräuker, “Oversampled adaptive sensing with random projections: Analysis and algorithmic approaches,” _Proc. IEEE In. Symp. on Signal Proc. and Inf. Tech. (ISSPIT)_ , pp. 336–341, December 2018, Louisville, USA.
  * [3] A. Bereyhi and R. R. Müller, “An adaptive bayesian framework for recovery of sources with structured sparsity,” _Proc. 8th Int. Work. on Comp. Adv. in Multi-Sensor Adap. Proc. (CAMSAP)_ , pp. 71–75, December 2019.
  * [4] J. D. Haupt, R. G. Baraniuk, R. M. Castro, and R. D. Nowak, “Compressive distilled sensing: Sparse recovery using adaptivity in compressive measurements,” _Proc. 43rd Asilomar Conf. on Signals, Systems and Computers_ , pp. 1551–1555, November 2009, Pacific Grove, CA, USA.
  * [5] J. Haupt, R. Nowak, and R. Castro, “Adaptive sensing for sparse signal recovery,” _Proc. 13th Digital Signal Proc. Work. and 5th IEEE Signal Proc. Edu. Work._ , pp. 702–707, January 2009, Marco Island, FL, USA.
  * [6] M. L. Malloy and R. D. Nowak, “Near-optimal adaptive compressed sensing,” _IEEE Trans. Inf. Th._ , vol. 60, no. 7, pp. 4001–4012, Jul. 2014.
  * [7] A. Bereyhi, R. R. Müller, and H. Schulz-Baldes, “Statistical mechanics of MAP estimation: General replica ansatz,” _IEEE Trans. Inf. Th._ , vol. 65, no. 12, pp. 7896–7934, Dec. 2019.
  * [8] A. Bereyhi, R. Müller, and H. Schulz-Baldes, “Rsb decoupling property of map estimators,” _Inf. Theory, International Workshop on (ITW)_ , 2016.

  *[SNR]: signal-to-noise ratio
  *[OAS]: oversampled adaptive sensing
  *[i.i.d.]: independent and identically dis-tributed
  *[MSE]: mean squared error


0  This is the author’s version of an article that has been published in this
journal. Changes were made to this version by the publisher prior to
publication.
The final version of record is available at
http://dx.doi.org/10.1109/LAWP.2018.2869926

# A Weak-Form Combined Source Integral Equation with Explicit Inversion of the
Combined-Source Condition

Jonas Kornprobst, , and Thomas F. Eibert Manuscript received February 2, 2018;
revised August 29, 2018; accepted August 29, 2018; date of this version August
29, 2018. _(Corresponding author: Jonas Kornprobst.)_ The authors are with the
Chair of High-Frequency Engineering, Department of Electrical and Computer
Engineering, Technical University of Munich, 80290 Munich, Germany (e-mail:
j.kornprobst@tum.de).Color versions of one or more of the figures in this
letter are available online at http://ieeexplore.ieee.org.Digital Object
Identifier 10.1109/LAWP.2018.2869926

###### Abstract

The combined source integral equation (CSIE) for the electric field on the
surface of a perfect electrically conducting scatterer can be discretized very
accurately with lowest-order Rao-Wilton-Glisson basis and testing functions if
the combined-source (CS) condition is enforced in weak form. We introduce a
technique to accelerate the iterative solution for this kind of CSIE. It is
demonstrated that the iterative solution of the equation system can be
performed very efficiently by explicitly inverting the weak-form CS condition
in any evaluation of the forward operator. This reduces the number of unknowns
and results in improved convergence behavior at negligible, linear cost.
Numerical results demonstrate that the new CSIE outperforms the classical CFIE
for high-accuracy simulations.

###### Index Terms:

electromagnetic scattering, Rao-Wilton-Glisson functions, combined source
integral equation, accuracy, low-order discretization

##  I Introduction

Electromagnetic radiation and scattering problems involving perfect
electrically conducting (PEC) objects are often treated by boundary integral
equations (IEs) due to their low discretization effort and good accuracy.
Common formulations for the tangential field components on the surface of the
scatterer are the electric field IE (EFIE) and the magnetic field IE (MFIE)
[1]. To obtain a system of equations with one unique solution, i.e. avoiding
the interior resonance problem, a combination of the two equations into the
combined field IE (CFIE) [2] is the most common solution.

However, a combination of electric and magnetic sources to the combined source
IE (CSIE) [3] is also possible; this involves the EFIE and MFIE operators as
well and is sometimes known under a different name, the Brackhage-Werner trick
[4, 5, 6, 7, 8]. A related concept (without physical currents) are the single-
source surface IEs (SSSIEs), which commonly only employ one type of current
unknowns [9, 10, 11, 12, 13, 14, 15].

With lowest-order div-conforming expansion functions on a triangular mesh, the
Rao-Wilton-Glisson (RWG) functions [16], the discretization of the MFIE
operator and also of the CFIE, is not accurate [17]. Therefore, the Buffa-
Christiansen (BC) functions, which reside in the dual space of the RWG
functions, have been successfully employed as testing functions for an
accurate mixed discretization of the MFIE [18]. For the CSIE as well, BC basis
functions have been demonstrated to produce accurate results [19]. Rather
recently, we demonstrated that an accurate and interior-resonance-free
formulation of the CSIE with pure RWG discretization is also possible with an
RWG-tested EFIE [20]. Thereby, the tangential electric field
$\bm{n}\times\bm{n}\times\bm{E}$ is tested with RWG functions and the combined
source (CS) side condition is also implemented with RWG and $\bm{n}\times$RWG
functions. Both electric and magnetic equivalent surface current densities are
modeled with the div-conforming RWG basis functions, which leads to double the
number of unknowns.

In this letter, we employ the CSIE discretization first in a dual-source
manner, but reduce it then to a discrete SSSIE: the weak-form CS condition is
explicitly solved in an iterative solver, leading to faster iterative solver
convergence. In addition, the number of unknowns is reduced by a factor of
two. First, we present the formulation and discretization of the CSIE. Next,
the discretization and implementation strategy of the CS condition is
discussed. Following that, simulation results demonstrate the excellent
accuracy and improved iterative solver convergence and conditioning of the
novel formulation of the RWG-discretized CSIE.

##  II The Combined Source Integral Equation

###  II-A Formulation

According to the Huygens principle, a PEC scattering scenario is described by
the boundary condition

$\bm{n}(\bm{r})\times\bm{E}(\bm{r})=\bm{n}(\bm{r})\times\bm{E}^{\mathrm{inc}}(\bm{r})+\bm{n}(\bm{r})\times\bm{E}^{\mathrm{sca}}(\bm{r})=\mathbf{0}\,,$
(1)

i.e. the superposition of the incident field $\bm{E}^{\mathrm{inc}}$ and the
scattered field $\bm{E}^{\mathrm{sca}}$ vanishes on the surface $S_{0}$ of the
scatterer. The CSIE employs both electric and magnetic surface current
densities $\bm{J}_{\mathrm{S}}(\bm{r})$ and $\bm{M}_{\mathrm{S}}(\bm{r})$ on
$S_{0}$ as equivalent sources to model the scattered fields. Expressing
$\bm{E}^{\mathrm{sca}}$ in terms of the equivalent currents for a suppressed
time factor $\mathrm{e}^{\,\mathrm{j}\omega t}$, the EFIE

$\displaystyle\bm{n}(\bm{r})\times\bm{E}^{\mathrm{inc}}(\bm{r})=\bm{n}(\bm{r})\times\oiint\nolimits_{S_{0}}\bm{\nabla}G_{0}(\bm{r},\bm{r}^{\prime})\times\bm{M}_{\mathrm{S}}(\bm{r}^{\prime})\mathop{}\mathrm{d}s^{\prime}\hskip
34.14322pt$
$\displaystyle+\frac{1}{2}\bm{M}_{\mathrm{S}}(\bm{r})+\mathrm{j}k_{0}Z_{0}\bm{n}(\bm{r})\times\left[\frac{1}{k_{0}^{2}}\bm{\nabla}\oiint\nolimits_{S_{0}}G_{0}(\bm{r},\bm{r}^{\prime})\bm{\nabla}^{\prime}\cdot\bm{J}_{\mathrm{S}}(\bm{r}^{\prime})\mathop{}\mathrm{d}s^{\prime}\right.$
$\displaystyle\left.+\oiint\nolimits_{S_{0}}G_{0}(\bm{r},\bm{r}^{\prime})\bm{J}_{\mathrm{S}}(\bm{r}^{\prime})\mathop{}\mathrm{d}s^{\prime}\right]$
(2)

is a more detailed version of (1). The involved scalar Green’s function of
free space is

$G_{0}(\bm{r},\bm{r}^{\prime})=\frac{\mathrm{e}^{-\mathrm{j}k_{0}\left|\bm{r}-\bm{r}^{\prime}\right|}}{4\uppi\left|\bm{r}-\bm{r}^{\prime}\right|}$
(3)

with source coordinate $\bm{r}^{\prime}$ and observation coordinate $\bm{r}$,
the free-space wave number $k_{0}$ and the free-space wave impedance $Z_{0}$.
Since the equivalent sources are not the Love currents, and thus ambiguous,
the CS constraint

$\bm{M}_{\mathrm{S}}(\bm{r})=\alpha
Z_{0}\bm{n}(\bm{r})\times\bm{J}_{\mathrm{S}}(\bm{r})$ (4)

can be enforced to obtain a unique combination of electric and magnetic
surface current densities, i.e. a quadratic system of equations with the same
number of unknowns as equations. The combination parameter $\alpha$ is
typically chosen as unity to achieve balanced contributions of electric and
magnetic currents to the radiated fields [3].

###  II-B Discretization

The discretization of (2) with RWG functions $\bm{\beta}_{n}$, defined on the
$n$th pair of neighboring triangles, is well-known. Suitable testing functions
in the dual space of the electric field $\bm{n}\times\bm{n}\times\bm{E}$ are
$\bm{\beta}$ functions. The MoM equation system for PEC objects is written as

$\left[-\frac{1}{2}\mathbfsfit A+\mathbfsfit
K\right]\small\textsf{{{v}}}+\mathrm{j}k_{0}Z_{0}\mathbfsfit
T\,\small\textsf{{{i}}}=\small\textsf{{{e}}}$ (5)

with the matrix elements of $\mathbfsfit K$ and $\mathbfsfit T$ defined in
[21] and the currents discretized as a summation of RWG basis functions

$\bm{J}_{\mathrm{S}}=\sum\nolimits_{n=1}^{N}\bm{\beta}_{n}\textit{{\small
i}}_{n}\,,\qquad\bm{M}_{\mathrm{S}}=\sum\nolimits_{n=1}^{N}\bm{\beta}_{n}\textit{{\small
v}}_{n}\,.$ (6)

Furthermore, we discretize the combined source condition (4) in the same way,
by expanding the currents according to (6) and testing the equation with RWG
functions. This gives the additional equation system

$\mathbfsfit{A^{\prime}}\small\textsf{{{v}}}-\alpha Z_{0}\mathbfsfit
A\small\textsf{{{i}}}=\small\textsf{{{0}}},$ (7)

where the Gram matrix $\mathbfsfit{A}$ linking $\bm{\beta}$ and
$\bm{n}\times\bm{\beta}$ functions has the entries

$\textit{{A}}_{mn}=\iint\nolimits_{S_{m}}\\!\bm{\beta}_{m}(\bm{r})\cdot\big{(}\bm{n}\times\bm{\beta}_{n}(\bm{r})\big{)}\mathop{}\\!\mathrm{d}s$
(8)

and the Gram matrix $\mathbfsfit{A}^{\prime}$ of the RWG functions is defined
with the entries

$\textit{{A}}_{mn}^{\prime}=\iint\nolimits_{S_{m}}\\!\bm{\beta}_{m}(\bm{r})\cdot\bm{\beta}_{n}(\bm{r})\mathop{}\\!\mathrm{d}s.$
(9)

Overall, we obtain the complete CSIE equation system

$\begin{bmatrix}\mathrm{j}k_{0}Z_{0}\mathbfsfit T&-\frac{1}{2}\mathbfsfit
A+\mathbfsfit K\\\ -\alpha Z_{0}\mathbfsfit
A&\mathbfsfit{A^{\prime}}\end{bmatrix}\begin{bmatrix}\small\textsf{{{i}}}\\\
\small\textsf{{{v}}}\end{bmatrix}=\begin{bmatrix}\small\textsf{{{e}}}\\\
\small\textsf{{{0}}}\end{bmatrix}.$ (10)

Since (10) has electric and magnetic current unknowns, it is referred to as
CSIE-JM in the following. In this formulation, as it has been presented in our
previous work [20], the discretized weak-form CS condition (7) augments the
EFIE as a side condition. As it will be demonstrated by numerical simulations,
this formulation has conditioning disadvantages compared to a formulation with
only one of the sources as unknowns. It is important to note that the
conditioning of (10) is improved considerably if the magnetic current unknowns
are replaced by the weighted version

$\small\textsf{{{v}}}^{\prime}=\small\textsf{{{v}}}/Z_{0}$ (11)

and the correct magnetic current coefficients are retrieved in the end by a
multiplication with $Z_{0}$. Alternatively, a diagonal preconditioner shows
similar benefits.

To improve the implementation of the CSIE, the Gram matrices are stored
separately and only the symmetric matrices $\mathbfsfit K$ and $\mathbfsfit T$
remain to be stored for the system matrix. The memory consumption of the Gram
matrices is negligible with linear complexity, since $\mathbfsfit{A}$ contains
$4N$ and $\mathbfsfit{A}^{\prime}$ $5N$ real-valued entries. To reduce the
number of unknowns, and halve the memory consumption of the iterative solver
of the whole system, which is a generalized minimum residual (GMRES) solver,
the side condition (7) is solved explicitly within each matrix vector product
according to

$\small\textsf{{{v}}}=\alpha Z_{0}\mathbfsfit{A^{\prime}}^{-1}\\!\mathbfsfit
A\mathbfsfit i\,,$ (12)

where ${\mathbfsfit{A}^{\prime}}^{-1}$ denotes the inverse of
${\mathbfsfit{A}^{\prime}}$ which is also obtained in an iterative manner by
using the conjugate gradient method. The Gram matrix $\mathbfsfit{A^{\prime}}$
is well-conditioned and the solver (with diagonal preconditioning) needs only
a few iterations for a low residual error. Due to the linear matrix size, the
convergence time of this inner iterative solver is negligible. We obtain the
CSIE with electric current unknowns only (CSIE-J)

$\left[\alpha Z_{0}\left(-\frac{1}{2}\mathbfsfit A+\mathbfsfit
K\right)\mathbfsfit{A^{\prime}}^{-1}\mathbfsfit
A+\mathrm{j}k_{0}Z_{0}\mathbfsfit
T\right]\small\textsf{{{i}}}=\small\textsf{{{e}}}\,.$ (13)

The magnetic currents are calculated in a post-processing step according to
(12) since it is not necessary to store them during the solution process.

In general, the easiest way to accelerate the iterative solver convergence is
diagonal preconditioning. For the CSIE-J, we investigate two possibilities.
One possibility is to neglect the influence of the magnetic currents on the
main diagonal and use only the diagonal matrix

$\mathbfsfit D=\mathrm{j}k_{0}Z_{0}\,\mathrm{diag}(\mathbfsfit{T})$ (14)

for the preconditioning of (13). Since the inverse of
$\mathbfsfit{A}^{\prime}$ is not directly available, only the main diagonal is
taken to estimate the inverse. Hence, the matrix

$\mathbfsfit D=\mathrm{j}k_{0}Z_{0}\mathrm{diag}(\mathbfsfit{T})-\alpha
Z_{0}\mathrm{diag}\left(\mathbfsfit{A}\,\mathrm{diag}(\mathrm{diag}(\mathbfsfit{A}^{\prime}))^{-1}\mathbfsfit{A}\right)$
(15)

can alternatively be employed for preconditioning.

##  III Numerical Results

For evaluating the performance of the two forms of the CSIE, they are compared
to the accurate but ill-conditioned EFIE and to the well-conditioned but
inaccurate MFIE. The common CFIE with combination parameter $0.5$ is also
considered, which is comparable to $\alpha=1$ in the CSIE.

First, the convergence behavior with mesh refinement is analyzed for a 1 m
square cube with plane wave illumination with a wavelength of 2 m. To analyze
the MFIE operator within the CSIE on its own, a combination parameter
$\alpha=10$ is chosen for the mesh refinement analysis. We consider three
different formulations of the CSIE: first the straight-forward version CSIE-JM
(10) as originally proposed in [20], then the modified CSIE-JM with weighted
magnetic currents according to (11), and finally the novel formulation (13).
In Fig. 1, the iterative solver convergence behavior to a residual error of
$10^{-4}$ is analyzed. The new CSIE-J formulation shows an acceptable
convergence, faster than the EFIE and almost as fast as the CFIE. For this
small example, in the absence of interior resonances, the CSIE-JM is much
slower than the EFIE and, thus, the formulation with the slowest iterative
solver convergence in this comparison. The block-diagonal preconditioning can
only overcome this in part. In Fig. 2, the averaged error of the scattered
far-field with respect to a higher-order EFIE solution on the finest mesh is
shown. As expected, both CSIE-J and CSIE-JM show almost exactly the same error
level, which is about 10 dB better than the error of the MFIE. The EFIE
solution is even about 5 dB more accurate. The performance difference of CFIE
and CSIE-J will be analyzed in the following more carefully in order to
highlight the benefits of the novel CSIE-J. It will be shown that only this
CSIE is able to outperform the CFIE.

Figure 1: Iterative solver convergence analysis for the mesh refinement of a 1
m square cube. Figure 2: Scattered far-field accuracy for the 1 m square cube.

A PEC sphere with a diameter of 1 m and 999 electric current unknowns is
considered. In the simulated frequency range, the average edge length of the
mesh is between $\lambda/25$ and $\lambda/7$. The first two interior
resonances, corresponding to the resonance frequencies of a spherical PEC
cavity with the same diameter, are observed for the EFIE and the MFIE. In Fig.
3, the iterative solver convergence of a GMRES solver is analyzed for a
residual error of $10^{-4}$. In the CSIE-J, the inversion of the Gram matrix
for the side condition is solved to a residual error of $10^{-5}$, which is
one order of magnitude lower. The EFIE, MFIE, CFIE and CSIE-J solutions are
obtained without preconditioning; the CSIE-JM (10) is solved without
preconditioning, with diagonal preconditioning and with block-diagonal
preconditioning according to (11). Both EFIE and MFIE suffer from the interior
resonances due to a null space at resonance frequencies. The CFIE and both
CSIE formulations show a stable convergence behavior, while only CSIE-J is
faster than the EFIE on average. In addition, the matrix condition numbers
without preconditioning (except for CSIE-JM) are shown in Fig. 4.

Figure 3: Iterative solver convergence analysis for the first two resonances
of a PEC sphere with a diameter of 1 m. Figure 4: Condition number around the
first two resonances of a 1 m sphere. Figure 5: Normalized singular value
spectra for the PEC sphere at $0.4$ GHz. For $J$ and $M$ unknowns, $M=1998$,
else $M=999$.

Again, the CFIE and CSIE formulations have the advantage of an absolutely
stable condition number at interior resonance frequencies. However, one
observation is remarkable: While the EFIE, CSIE-J and both preconditioned
CSIE-JM equations show almost the same condition number, the iterative solver
convergence of the CSIE-J is much better. The most reasonable cause is a
different eigenvalue distribution in the complex plane, e.g. by an improved
clustering more suitable for GMRES. In the singular-value spectrum, see Fig.
5, we can recognize a certain amount of singular values with larger values for
the CSIE-J as compared to the EFIE. The better iterative solver convergence of
CSIE-J has been verified for many different right-hand sides, including random
vectors and unit vectors of the standard basis. Furthermore, the accuracy as
compared to the analytical Mie series solution is shown in Fig. 6. It is
observed that the MFIE shows the poorest accuracy, with huge deteriorations at
the resonance frequencies. Compared to this behavior, the EFIE only shows
inaccuracies quite exactly at the resonance frequencies. The CFIE is
influenced by the MFIE error and, thus, its accuracy is worse than the CSIE.
Both CSIE formulations show absolutely the same error level, as expected,
which is also about the same as the EFIE error.

Figure 6: Accuracy analysis around the first two resonances of a PEC sphere.

As a further problem, the stealth object Flamme is simulated with $\lambda/10$
mean edge length discretization and $287\,922$ electric current unknowns [22,
23]. The MFIE takes $1267$ iterations for convergence to a residual error of
$10^{-5}$, the CFIE $315$, the CSIE-J $448$, and the CSIE-JM $1264$
iterations. The EFIE did not converge without preconditioning even after
thousands of iterations. With diagonal preconditioning, the iterative solution
of the MFIE is accelerated to $899$ iterations and the iterative solution of
the CFIE to $253$. The CSIE-J converges within $384$ iterations if only the
EFIE operator is taken for preconditioning and within $381$ iterations if the
approximated diagonal according to (15) is utilized. The average time per
matrix-vector product is $10.5$ seconds for EFIE, MFIE, and CFIE, and $20.7$
seconds for the CSIE-JM, where all computations have been performed single-
threaded on an Intel E5-1650 v4 processor with $3.6$ GHz clock speed. The
explicit inversion of the CS condition takes only additional $0.04$ seconds on
average, or about $10$ iterations, for a residual error of $5\cdot 10^{-7}$.
Results of the bistatic radar cross section are shown in Fig. 7. Compared to
an 3rd order EFIE solution on a refined mesh [21], the MFIE shows a maximum
error of $-37.3$ dB, the CFIE of $-44.2$ dB, both CSIE-J and CSIE-JM of
$-57.8$ dB, and the EFIE of $59.8$ dB.

With a difference in the iterative solver convergence, it is still not clear
whether the CSIE-J or the CFIE offer the better accuracy-convergence trade-
off. Therefore, we analyze the influence of the corresponding weighting or
combination factors for the cases of CSIE-J and CFIE with diagonal
preconditioning and CSIE-JM with block-diagonal preconditioning in Fig. 8. The
simulated object is again the Flamme with the same parameters as before. It
becomes clear that the CSIE-JM can provide the same error level only at the
cost of more solver iterations. Additionally, the CSIE-J is able to outperform
the CFIE for near-EFIE accuracies. For instance, the CSIE-J with $\alpha=2.7$
converges within 313 iterations and offers an error level of $56.6$ dB, while
the CFIE with combination parameter $0.9$ (90% EFIE) converges in 758
iterations for the same error. This still means less than half the solver
memory consumption for the CSIE-J, approximately the same memory for the
matrices, and about twice the time per matrix-vector-product (but less
solution time overall).

Figure 7: Side-lobe section of the bistatic radar cross section $\sigma$ of
the Flamme for the EFIE, CFIE, and CSIE solutions. Figure 8: Error level vs.
iterations with different weightings in CFIE and CSIE for Flamme.

##  IV Conclusion

A new form of the CSIE with weak-form CS condition has been introduced,
showing an improved iterative solver convergence and a reduced number of
unknowns. The magnetic current unknowns are explicitly computed from the
electric currents for each matrix-vector product within the iterative solver
and, therefore, eliminated. The new CSIE is able to offer a better
accuracy/convergence trade-off than the CFIE and the weighting factor in the
CSIE is not a critical choice (other than in the CFIE) in order to achieve
almost EFIE accuracy.

## References

  * [1] R. F. Harrington, “Boundary integral formulations for homogeneous material bodies,” _J. Electromagn. Waves Appl._ , vol. 3, no. 1, pp. 1–15, Jan. 1989\.
  * [2] J. R. Mautz and R. F. Harrington, “H-field, E-field, and combined-field solutions for conducting bodies of revolution,” _A.E.Ü. (Germany)_ , vol. 32, no. 4, pp. 157–164, Apr. 1978.
  * [3] J. Mautz and R. Harrington, “A combined-source solution for radiation and scattering from a perfectly conducting body,” _IEEE Trans. Antennas Propag._ , vol. 27, no. 4, pp. 445–454, Jul. 1979.
  * [4] H. Brakhage and P. Werner, “Über das dirichletsche Außenraumproblem für die Helmholtzsche Schwingungsgleichung,” _Archiv der Mathematik_ , vol. 16, no. 1, pp. 325–329, 1965.
  * [5] A. Buffa and R. Hiptmair, “Regularized combined field integral equations,” _Numerische Mathematik_ , vol. 100, no. 1, pp. 1–19, 2005.
  * [6] M. Darbas, “Generalized combined field integral equations for the iterative solution of the three-dimensional Maxwell equations,” _Appl. Math. Lett._ , vol. 19, no. 8, pp. 834–839, 2006.
  * [7] O. Steinbach and M. Windisch, “Modified combined field integral equations for electromagnetic scattering,” _SIAM J. Numer. Anal._ , vol. 47, no. 2, pp. 1149–1167, 2009.
  * [8] J. M. Melenk, “Mapping properties of combined field Helmholtz boundary integral operators,” _SIAM J. Math. Anal._ , vol. 44, no. 4, pp. 2599–2636, 2012.
  * [9] E. Marx, “Single integral equation for wave scattering,” _J. Math. Phys._ , vol. 23, no. 6, pp. 1057–1065, Jun. 1982.
  * [10] A. Glisson, “An integral equation for electromagnetic scattering from homogeneous dielectric bodies,” _IEEE Trans. Antennas Propag._ , vol. 32, no. 2, pp. 173–175, Feb 1984.
  * [11] M. S. Yeung, “Single integral equation for electromagnetic scattering by three-dimensional homogeneous dielectric objects,” _IEEE Trans. Antennas Propag._ , vol. 47, no. 10, pp. 1615–1622, Oct. 1999.
  * [12] A. Menshov and V. Okhmatovski, “New single-source surface integral equations for scattering on penetrable cylinders and current flow modeling in 2-D conductors,” _IEEE Trans. Microw. Theory Techn._ , vol. 61, no. 1, pp. 341–350, Jan. 2013.
  * [13] Y. Shi and C. Liang, “An efficient single-source integral equation

solution to EM scattering from a coated conductor,” _IEEE Antennas Wireless
Propag. Lett._ , vol. 14, pp. 547–550, 2015.

  * [14] U. R. Patel, P. Triverio, and S. V. Hum, “A novel single-source surface integral method to compute scattering from dielectric objects,” _IEEE Antennas Wireless Propag. Lett._ , vol. 16, pp. 1715–1718, 2017.
  * [15] F. S. H. Lori, A. Menshov, R. Gholami, J. B. Mojolagbe, and V. I. Okhmatovski, “Novel single-source surface integral equation for scattering problems by 3-D dielectric objects,” _IEEE Trans. Antennas Propag._ , vol. 66, no. 2, pp. 797–807, Feb. 2018.
  * [16] S. M. Rao and D. R. Wilton, “E-field, H-field, and combined field solution for arbitrarily shaped three-dimensional dielectric bodies,” _Electromagn._ , vol. 10, no. 4, pp. 407–421, Oct. 1990.
  * [17] Ö. Ergül and L. Gürel, “Investigation of the inaccuracy of the MFIE discretized with the RWG basis functions,” in _Proc. IEEE Antennas Propag. Soc. Symp. (APSURSI)_ , Monterey, CA, Jun. 2004, pp. 3393–3396.
  * [18] K. Cools, F. P. Andriulli, D. De Zutter, and E. Michielssen, “Accurate and conforming mixed discretization of the MFIE,” _IEEE Antennas Wireless Propag. Lett._ , vol. 10, pp. 528–531, 2011.
  * [19] P. Ylä-Oijala, S. P. Kiminki, K. Cools, F. P. Andriulli, and S. Järvenpää, “Stable discretization of combined source integral equation for scattering by dielectric objects,” _IEEE Trans. Antennas Propag._ , vol. 60, no. 5, pp. 2575–2578, May 2012.
  * [20] J. Kornprobst and T. J. Eibert, “A combined source integral equation with weak form combined source condition,” _IEEE Trans. Antennas Propag._ , vol. 66, no. 4, pp. 1251–2155, Apr. 2018.
  * [21] Ismatullah and T. F. Eibert, “Surface integral equation solutions by hierarchical vector basis functions and spherical harmonics based multilevel fast multipole method,” _IEEE Trans. Antennas Propag._ , vol. 57, no. 7, pp. 2084–2093, Jul. 2009.
  * [22] L. Gürel, H. Bağcı, J. C. Castelli, A. Cheraly, and F. Tardivel, “Validation through comparison: Measurement and calculation of the bistatic radar cross section of a stealth target,” _Radio Sci._ , vol. 38, no. 3, Jun. 2003.
  * [23] T. F. Eibert, “A diagonalized multilevel fast multipole method with spherical harmonics expansion of the k-space integrals,” _IEEE Trans. Antennas Propag._ , vol. 53, no. 2, pp. 814–817, Feb. 2005.