The following article has been submitted to the Journal of the Acoustical
Society of America.

After it is published, it will be found at
http://asa.scitation.org/journal/jas.

1 0 current page.north east -6.0cm -0.4cm

# Gain and phase calibration of sensor arrays from ambient noise by cross-
spectral measurements fitting

Charles Vanwynsberghe<EMAIL_ADDRESS>Technology Innovation
Institute, 9639 Masdar City, Abu Dhabi, United Arab Emirates  Simon Bouley
MicrodB, F-69134, Écully, France  Jérôme Antoni  Univ Lyon, INSA Lyon, LVA,
EA677, 69621, Villeurbanne, France

###### Abstract

We address the problem of blind gain and phase calibration of a sensor array
from ambient noise. The key motivation is to ease the calibration process by
avoiding a complex procedure setup. We show that computing the sample
covariance matrix in a diffuse field is sufficient to recover the complex
gains. To do so, we formulate a non-convex least-square problem based on
sample and model covariances. We propose to obtain a solution by low-rank
matrix approximation, and two efficient proximal algorithms are derived
accordingly. The first one solves the problem modified with a convex
relaxation to guarantee that the solution is a global minimizer, and the
second one directly solves the initial non-convex problem. We investigate the
efficiency of the proposed algorithms by both numerical and experimental
results according to different sensing configurations. These show that
efficient calibration highly depends on how the measurements are correlated.
That is, estimation is achieved more accurately when the field is spatially
over-sampled.

††preprint: Author, JASA

##  I Introduction

Calibration refers to the task of fixing the gain of a sensor accurately. The
task is conventionally done on a single element, but can also be performed on
a distributed network [14, 5] or a synchronous array. The latter configuration
is considered in this paper. Several acoustical problems rely on the use of
sensor arrays (e.g. source localization and identification) and are therefore
sensitive to potential model uncertainties [20]. Both amplitudes and phases of
sensor gains are a part of such critical uncertainties [18, 46, 11, 12].

Classically, gain calibration is done on a single element either by the
reciprocity or the comparison method [39, ch. 24]. The reciprocity method is
adopted in standardized protocols to calibrate reference microphones [17]. In
spite of its reliability, going for an element-wise calibration with a large
sensor array is cumbersome, if not infeasible in practice. Consequently,
attempts were made to calibrate the whole array in one shot. The overall
picture is as follows in the noiseless case: provided $1$ sensor reading, one
needs to discriminate the (unknown) sensor gain from the (unknown) actual
value of the physical field itself, resulting in $2$ unknown parameters for
$1$ reading. By extension, $M$ sensor readings correspond to $2M$ unknowns, so
the problem is under-determined for any array of $M$ sensors. However,
reducing the number of unknowns can be tackled in different ways. A first
example is to rely on the range space deficiency in which the physical field
lies [4, 3]. A sequence of well-calibrated measurements is then supposed to
lie in a signal subspace whose dimension is small compared with $M$. A second
example rather assumes that the physical field is generated by a small number
of sources, and is based on sparse approximation [7, 25]. Moreover, if sensor
gains are somehow correlated by the underlying physical phenomena (e.g.
temperature), they can also be constrained into a subspace [31, 30] (i.e. of
low dimension as compared to $M$).

One context that has rarely been explored in acoustics is in situ calibration
with an ambient noise physical field. It is of potential interest to several
applications, such as in underwater acoustics (e.g. source monitoring [34],
passive fathometer [43], geoacoustic inversion [49], Green function retrieval
[40], thermocline [29, 21] or temperature [47, 44] estimation), in seismic
sensing (e.g. tomography [42], detection [41]), or structural health
monitoring by tomography [15]. Other inverse problems applied to imaging are
listed in [19]. This non-exhaustive list depicts the diversity of applications
where sensor gain uncertainty surely causes estimation issues and, therefore,
should be mitigated.

Rare works port the attention on how to benefit from ambient noise to
calibrate sensors. In [37], a seismometer is calibrated with ground noise, and
placed in proximity to a supplementary seismometer of reference whose absolute
gain is known. This is the aforementioned comparison method: it relies on the
calculus of the transfer function between both seismometers. It assumes that
the latter record the same physical field, and it is a sensor-wise procedure.
More recent efforts were made by Abkar et al. to obtain the microphone gains
[2, 1] of a spherical array, with the help of a diffuse field. These authors
emphasize the practicality of this approach because it is easier to generate a
spatially scattered acoustic noise than a controlled source (or more) in free
field conditions. However, the suggested method is energy-based, and provides
the amplitude part of the gain only.

We propose a sensor gain and phase calibration where the “ambient noise”
precisely refers to a diffuse field. In such conditions, cross-spectral
densities between sensors can be expressed by a well-known second-order
statistic model [35, 27]. When the cross-spectral readings are arranged in the
form of a covariance matrix, we show that the problem can be turned into a
low-rank matrix approximation problem[13]. To the best of our knowledge, the
solution proposed in this paper is distinct from the state-of-the art by the
fact that it relies on second-order statistical measurements (the sample
covariance matrix) rather than on a collection of snapshot vectors. As a
consequence, the complexity is independent of the observation time duration
and, due to the averaging of noise, good performance is reached even with low
signal-to-noise ratios. Moreover, the low-rank matrix approximation problem is
often addressed in contexts where dimensions are large (e.g. in machine
learning [13]), and efficient methods were proposed accordingly. Leveraging
low-complexity methods is relevant here, because it allows the simultaneous
calibration of a large number of sensors.

The paper is organized as follows. Section II describes the sensing model
based on the diffuse field assumption. The trivial derivation of a gain
estimator is presented with its main limitations for a practical use in Sec.
III, which will serve as the baseline in experiments. Then, Sec. IV defines
the least-square problem in the original non-convex form, followed by a convex
relaxation to ensure that the obtained solution is the global minimizer. Since
the estimated gains are relative up to a constant, Sec. IV.3 discusses
different options for relative to absolute value conversion. Finally, two
efficient proximal algorithms [36] that solve the non-convex and convex
problems are derived in Sec. V. The performance of these solvers with regard
to signal-to-noise ratio, acquisition time, model mismatch and the density of
spatial sampling are studied via simulations in Sec. VI. At last, an
experimental validation is presented in an acoustic reverberant environment to
calibrate a $43$-microphone array in Sec. VII. A final conclusion summarizes
the whole results and suggests relevant practices to perform calibration in
ambient noise.

### Mathematical notations

The notations in the paper are as follows. Bold lowercase (resp. uppercase)
letters describe vectors (resp. matrices). The identity matrix of size
$N\times N$ is defined by $\mathbf{I}_{N}$; $\mathbf{0}_{N}$ stands for the
column vector of $N$ zeros, and $\mathbf{1}_{M\times N}$ stands for the
$M\times N$ matrix of ones. $\operatorname{diag}(\bm{v})$ is the diagonal
square matrix whose diagonal contains the elements of the vector $\bm{v}$.
$\odot$ and $\oslash$ respectively stands for the element-wise product and
division of two equal-size matrices. $\|\mathbf{M}\|_{\mathsf{F}}$ denotes the
Frobenius norm of the matrix $\mathbf{M}$. $\mathbf{M}^{\mathsf{H}}$ stands
for the transpose-conjugate of matrix $\mathbf{M}$.

##  II Covariance-based sensing model in ambient noise

Throughout this paper, the considered scenario involves an array of $M$
sensors providing pointwise and omnidirectional measurements in a homogeneous
media, meaning that the propagation speed $c_{0}$ is constant. One given
frequency $f$ relates to one wavelength value, denoted $\lambda$. The sensing
model is expressed in the frequency domain after applying the Discrete Fourier
Transform (DFT) on a collection of time series to obtain $L$ snapshots indexed
by $l$.

This paper considers the sensor calibration problem for which an unknown
frequency-dependent multiplicative and complex-valued gain affects the output
value of each sensor. In the noisy case, the considered model for the $l$-th
snapshot at frequency $f$ reads

$\mathbf{y}_{l}=\operatorname{diag}(\bm{\alpha})\mathbf{x}_{l}+{\mathbf{w}}_{l},$
(1)

where $\bm{\alpha}\in\mathbb{C}^{M\times 1}$ denotes the set of unknown sensor
gains, $\mathbf{y}_{l}\in\mathbb{C}^{M\times 1}$ are the values read by the
sensors, $\mathbf{x}_{l}$ are the true values of the field, and
${\mathbf{w}}_{l}$ denotes the additive sensor noise (frequency dependence
will be systematically omitted for ease of notation). By considering complex
gains, both amplitude and phase variations are taken into account.

The physical field and the sensor noise realizations are supposed to be two
independent, zero-mean, random processes. The covariance matrix of the
physical field is
$\operatorname{\mathrm{E}}\\{\mathbf{x}_{l}\mathbf{x}_{l}^{\mathsf{H}}\\}=\sigma_{x}^{2}\bm{\Sigma}_{x}$,
for any $l$, where $\sigma_{x}^{2}$ denotes the variance of the field and
$\bm{\Sigma}_{x}\in\mathbb{C}^{M\times M}$ contains the coherences between
sensors. Similarly, the covariance matrix of sensor noise is
$\sigma_{w}^{2}\mathbf{I}_{M}$, which reflects independent variables, each
with the same variance $\sigma_{w}^{2}$.

When the physical field is ambient noise111From this point, we strictly
distinguish ambient noise (characterizing the physical field) and sensor noise
(i.e. the additive non-propagative noise which originates from the acquisition
device) in order to avoid any ambiguity between both., a well-established
model consists in considering that the field is the integral of uncorrelated
plane waves of homogeneous variance over the sphere. The latter is known as
diffuse field [27]. If uncorrelated plane waves have variance $\sigma_{x}^{2}$
over the whole sphere, it reads [26]:

$[\bm{\Sigma}_{x}]_{mn}=\operatorname{sinc}\left(\dfrac{2\pi}{\lambda}d_{mn}\right),$
(2)

where $d_{mn}$ is the distance between the $m$-th and the $n$-th sensors.
Also, the covariance matrix from the sensor readings can be derived according
to the snapshot model in Eq. (1):

$\bm{\Sigma}_{y}=\operatorname{\mathrm{E}}\\{\mathbf{y}_{l}\mathbf{y}_{l}^{\mathsf{H}}\\}=\sigma_{x}^{2}\operatorname{diag}(\bm{\alpha})\bm{\Sigma}_{x}\operatorname{diag}(\bm{\alpha})^{\mathsf{H}}+\sigma_{w}^{2}\mathbf{I}_{M}.$
(3)

This covariance can be approximated by the measurements from the zero-mean
sensor readings. The sample covariance matrix, denoted
$\widetilde{\bm{\Sigma}}_{y}$, is estimated from the $L$ snapshots as

$\widetilde{\bm{\Sigma}}_{y}=\dfrac{1}{L}\sum_{l=1}^{L}\mathbf{y}_{l}\mathbf{y}_{l}^{\mathsf{H}}.$
(4)

We emphasize that the use of DFT to build the snapshots $\mathbf{y}_{l}$ is
impacting and should be well-tuned. Indeed, $\bm{\Sigma}_{y}$ describes the
diffuse field for the frequency $f$, whereas $\widetilde{\bm{\Sigma}}_{y}$
results from DFT snapshots, and the DFT acts as a low-pass filter with a
precision frequency denoted $\Delta f$. This results into an integral over a
frequency interval of width $\Delta f$ that is not taken into account in the
derived model (2). Keeping consistency between the model $\bm{\Sigma}_{y}$ and
the measurement $\widetilde{\bm{\Sigma}}_{y}$ requires that DFT is
sufficiently narrowband [35]. Then, respecting

$\Delta f\ll\dfrac{c_{0}}{d_{mn}}\>\forall\>m,n$ (5)

is a necessary condition to preserve the validity of the model (2). In the
acoustical context, this condition should be easily met.

An interesting aspect of using the sample covariance matrix rather than
snapshots is to reduce computation complexity and memory usage: while the
amount of data to process equals $LM$ and scales linearly with the number of
snapshots, it remains constant (and equals $M^{2}$) with the sample covariance
matrix.

##  III Vanilla solution under strong conditions

In this section, we show that a straightforward solution can be derived under
the noiseless assumption. If $\sigma_{w}^{2}\ll\sigma_{x}^{2}$, plugging in
the sample covariance matrix in Eq. (3) rewritten as an element-wise matrix
product, results in the relation
$\widetilde{\bm{\Sigma}}_{y}\approx\sigma_{x}^{2}\left(\bm{\alpha}\bm{\alpha}^{\mathsf{H}}\right)\odot\bm{\Sigma}_{x}$.
Then, a simple vanilla estimation of $\bm{\alpha}$ consists in obtaining the
principal eigenvector of the matrix
$\widetilde{\bm{\Sigma}}_{y}\oslash\bm{\Sigma}_{x}$. This solution provides
relative gains, i.e. up to a complex constant – this point will be discussed
in detail in Sec. IV.3.

The idea is simple, but has a caveat: it is unstable in practice because of
the element-wise division. Instability occurs if the coefficient
$[\bm{\Sigma}_{x}]_{mn}$ is small for some given index pair $(m,n)$, so that
the presence of an additive error term in $[\widetilde{\bm{\Sigma}}_{y}]_{mn}$
is magnified by the division operation. This case exists for example if
$\sigma_{w}^{2}\neq 0$, i.e. even with low sensor noise. The same issue arises
if a small additive error term exists in $[\bm{\Sigma}_{x}]_{mn}$, i.e. in the
model itself. Thus, there is potentially a high sensitivity to sensor noise,
but also to model mismatch.

A simple rule of thumb to safely use the vanilla approach is to keep

$\lambda\gg d_{mn}\;\text{for}\;1\leq m,n\leq M.$ (6)

This condition prevents the $\operatorname{sinc}$ function in (2) from being
close to zero, and leads to the approximation
$\bm{\Sigma}_{x}\approx\mathbf{1}_{M\times M}$. The condition (6) describes a
low frequency regime, and is met when the array aperture is relatively smaller
than the wavelength.

Since, to the best of our knowledge, no state-of-the-art work has addressed
complex gain calibration from covariance readings, the vanilla solution will
be used in the results as the baseline. In the next sections, we focus on the
proposed calibration method of this paper, which does not require a condition
as strict as Eq. (6).

##  IV Complex gain calibration by covariance fitting

###  IV.1 The least square approach

We propose to estimate $\bm{\alpha}$ by solving a least-square problem. In the
general case, the gain values are relative, i.e. up to a complex constant. The
least square minimization can be written by plugging in the sample covariance
matrix in the sensing model (3):

$\operatornamewithlimits{argmin}_{\bm{\alpha}}\dfrac{1}{2}\|\widetilde{\bm{\Sigma}}_{y}-\operatorname{diag}(\bm{\alpha})\bm{\Sigma}_{x}\operatorname{diag}(\bm{\alpha})^{\mathsf{H}}\|^{2}_{\mathsf{F}},$
(7)

the term $\sigma_{x}^{2}$ being discarded without loss of generality, since we
seek relative gains in the first time – see Sec. IV.3 for more details. This
form is non-convex over $\bm{\alpha}$, and to our knowledge, no closed-form
solution can be derived. Let $l$ be the least-square loss function; we remark
that $l$ can also be written under the form

$l(\mathbf{C})=\dfrac{1}{2}\|\widetilde{\bm{\Sigma}}_{y}-\mathbf{C}\odot\bm{\Sigma}_{x}\|^{2}_{\mathsf{F}},$
(8)

turning the input variable $\mathbf{C}=\bm{\alpha}\bm{\alpha}^{\mathsf{H}}$
into a semi-definite positive matrix of rank 1. Thus, the least-square problem
(7) can be recast in the following matrix form:

$\widetilde{\mathbf{C}}=\operatornamewithlimits{argmin}_{\mathbf{C}\in\mathcal{S}_{M}^{+}}l(\mathbf{C})\text{
subject to }\operatorname{rank}(\mathbf{C})=1,$ (9)

with $\mathcal{S}_{M}^{+}$ the set of semi-definite positive matrices. The
relative gains $\tilde{\bm{\alpha}}$ can be calculated straightforwardly from
the principal eigenvector of $\widetilde{\mathbf{C}}$. By applying the
eigenvalue decomposition (EVD)
$\widetilde{\mathbf{C}}=\sum_{m=1}^{M}s_{m}\mathbf{u}_{m}\mathbf{u}_{m}^{\mathsf{H}}$
such that $\mathbf{u}_{m}$ are the eigenvectors and $s_{m}$ are the
eigenvalues ranked in the decreasing order as $s_{1}\geq s_{2}\geq
s_{3}\dots\geq s_{M}$, we choose $\tilde{\bm{\alpha}}=\mathbf{u}_{1}$.

Note that problem (9) is still non-convex, now over $\mathbf{C}$, due to the
rank constraint. Likewise, the loss in Eq. (8) can also be seen as a weighted
case of matrix factorization with factors $\bm{\alpha}$ and
$\bm{\alpha}^{\mathsf{H}}$, and is known to be non-convex as well [13].
However, this reformulation is crucial for the sequel of the paper, because it
unlocks the access to the low-rank matrix approximation, and the possibility
to use efficient solvers accordingly.

###  IV.2 Penalized least-square and convex relaxation

Finding the global minimizer of the problem (9) is generally a NP-hard problem
[13]. However, it is possible to obtain a tractable alternative by applying
convex relaxation to the program. To do so, the rank constraint on
$\mathbf{C}$ is replaced by its convex envelope [16]: the nuclear norm
$\|\mathbf{C}\|_{*}$, defined by the sum of singular values of $\mathbf{C}$.
This modification is enough to complete the convex relaxation, since (i) the
set $\mathcal{S}_{M}^{+}$ is already convex, and (ii) the loss term
$l(\mathbf{C})$ is convex on $\mathbf{C}$.

Rather than choosing a problem defined with constraint as in (9), we favor a
penalized least squares expression. We define the following problem
$\mathcal{P}_{i}~{}$as:

$\widetilde{\mathbf{C}}=\operatornamewithlimits{argmin}_{\mathbf{C}\in\mathcal{S}_{M}^{+}}l(\mathbf{C})+\nu\operatorname{\Omega}_{i}(\mathbf{C}),$
(10)

where $\operatorname{\Omega}_{i}(\mathbf{C})$ is a penalty function (also
known as a regularizer), and $\nu$ is a non-zero positive parameter that
controls the importance of the penalty.

Let us define the convex form of the penalized least square for $i=1$, i.e.
$\mathcal{P}_{1}~{}$, associated to the penalty term
$\operatorname{\Omega}_{1}(\mathbf{C})=\|\mathbf{C}\|_{*}$. The global
minimizer of $\mathcal{P}_{1}~{}$can be easily found via a toolbox like CVX
[24]; however, the merit of using the penalized unconstrained form is to
derive an efficient solver, such as proximal algorithm [36], Alternating
Direction Method of Multipliers [9], or Frank-Wolfe algorithm [28]. A solver
for $\mathcal{P}_{1}~{}$will be derived accordingly in section V.

Finally, for $i=0$ the penalized least-square program $\mathcal{P}_{0}~{}$can
be defined as an equivalent to Eq. (9), if we choose
$\operatorname{\Omega}_{0}(\mathbf{C})=I_{\mathcal{C}}(\mathbf{C})$, such that
$I_{\mathcal{C}}(.)$ is the indicator function defined as [36]

$\displaystyle\begin{split}I_{\mathcal{C}}(\mathbf{M})&=0\text{ if
}\mathbf{M}\in\mathcal{C},\\\ &=\infty\text{ otherwise,}\end{split}$ (11)

and $\mathcal{C}=\\{\mathbf{M}\in\mathbb{C}^{M\times
M}|\operatorname{rank}(\mathbf{M})=1\\}$. Although with only a guarantee to
find a local minimizer, we will see further in Sec. V that an efficient solver
can be derived as well, because the problem $\mathcal{P}_{i}~{}$is a
generalized convenient form to switch between $\mathcal{P}_{0}~{}$and
$\mathcal{P}_{1}~{}$.

###  IV.3 From relative to absolute gains

As seen in Sec. III and IV.1, the provided estimation $\tilde{\bm{\alpha}}$ is
relative up to a constant in the general case. In some cases the relative
solution could be enough (e.g. in source localization). If one wants the
absolute levels of sources, then a scaling correction is needed. We define
$\tilde{\bm{\alpha}}_{\text{sc}}=\eta\tilde{\bm{\alpha}}$ the estimated gain
vector after scaling, and $\eta$ the scaling correction. Below, we present
some case examples of how to treat this issue.

#### Scaling with reference sensor(s)

A first option is to consider the addition of one reference sensor. From
asynchronous measurement of the power spectral density, it is then possible to
estimate the variance of the diffuse field $\sigma_{x}^{2}$, because the
latter is spatially homogeneous. In this case, we choose
$\eta=\sigma_{x}^{-2}$, yet the scaling is on amplitude only. Alternatively,
an amplitude and phase correction is possible if the reference microphone is a
part of the synchronous array (e.g. if it is the 1-st microphone with gain 1,
then $\eta=1/[\tilde{\bm{\alpha}}]_{1}$). If several reference microphones are
included, their gains can be incorporated as anchors by a modification of the
loss term – see [4, sec 3.2] for more details.

#### Scaling with respect to the gain value expectation

A second option is to rely on the expected gain of the whole sensor
collection. This is relevant when the manufacturer provides gain
specifications with uncertainty (e.g. the mean and standard deviation in the
hardware datasheet). Let $\bm{\alpha}_{0}$ be the vector of expected gains.
The optimal correction in the least square sense aligns $\tilde{\bm{\alpha}}$
onto $\bm{\alpha}_{0}$ in the complex plane (like in Procrustes analysis
[23]). Then the scalar $\eta$ minimizes

$\eta=\operatornamewithlimits{argmin}_{a\in\mathbb{C}}\|\bm{\alpha}_{0}-a\tilde{\bm{\alpha}}\|_{2}^{2},$
(12)

so

$\eta=\dfrac{\tilde{\bm{\alpha}}^{\mathsf{H}}\bm{\alpha}_{0}}{\|\tilde{\bm{\alpha}}\|_{2}^{2}}.$
(13)

If all sensors are manufactured through the same process and have the gain
mean value $\overline{\alpha}$, then
$\bm{\alpha}_{0}=\overline{\alpha}\mathbf{1}_{M}$.

#### Scaling with respect to the ground truth

The last option is to use the ground truth as the reference, so
$\bm{\alpha}_{0}$ equals the true gain vector. This is irrelevant in practice,
but is important to assess the performance of the calibration method only.
Indeed, in that case, it is independent of a user-informed or anchor-based
scaling step. This will be used in result sections VI and VII to evaluate the
performance of the calibration methods.

##  V Sensor gain calibration solvers

The problem 10 with the sum of a regression loss and a penalty term is a well-
known general form met in many fields, such as machine learning [13] or
computational imaging ill-posed problems [45]. These applications usually meet
large-scale problems, and require efficient solvers to deal with high
dimensions. As a result, to tackle the non-smoothness of the penalty function,
and to deal with large-scale problems, a plethora of dedicated optimization
methods have been proposed. Below, we propose two solvers for
$\mathcal{P}_{0}~{}$and $\mathcal{P}_{1}~{}$, based on the first-order
iterative proximal algorithm scheme [36, 6].

The proximal gradient descent is an iterative algorithm based on the following
update step [6, ch. 10.2]:

$\mathbf{C}^{[k+1]}=\operatorname{\textbf{prox}}_{\gamma\nu\operatorname{\Omega}_{i}}\left(\mathbf{C}^{[k]}-\gamma\nabla
l(\mathbf{C}^{[k]})\right)$ (14)

where $k$ is the iteration index and $\gamma$ is the gradient step. It
performs a gradient descent, followed by the proximal mapping with the
operator
$\operatorname{\textbf{prox}}_{\gamma\nu\operatorname{\Omega}_{i}}(.)$ defined
by the program

$\displaystyle\operatorname{\textbf{prox}}_{a\operatorname{\Omega}_{i}}(\mathbf{Z})=\operatornamewithlimits{argmin}_{\mathbf{M}}\frac{1}{2}\|\mathbf{M}-\mathbf{Z}\|^{2}_{\mathsf{F}}+a\operatorname{\Omega}_{i}(\mathbf{M})$
(15)

in the case of the mapping of a matrix $\mathbf{Z}$, with $a$ being some
positive real scalar. As a result, the proximal algorithm is flexible, because
solving $\mathcal{P}_{1}~{}$and $\mathcal{P}_{0}~{}$only needs to choose the
proximal operator accordingly. Finally, the loss gradient $\nabla l$ is
obtained with the element-wise Wirtinger derivatives for $1\leq m,n\leq M$:

$\displaystyle\dfrac{\partial l(\mathbf{C})}{\partial[\mathbf{C}]_{mn}^{*}}$
$\displaystyle=\dfrac{\partial}{\partial[\mathbf{C}]_{mn}^{*}}\left(\frac{1}{2}\left|[\widetilde{\bm{\Sigma}}_{y}]_{mn}-[\mathbf{C}]_{mn}[\bm{\Sigma}_{x}]_{mn}\right|^{2}\right)$
(16)
$\displaystyle=-\frac{1}{2}[\bm{\Sigma}_{x}]_{mn}^{*}\left([\widetilde{\bm{\Sigma}}_{y}]_{mn}-[\mathbf{C}]_{mn}[\bm{\Sigma}_{x}]_{mn}\right).$

In the end, the gradient reads

$\nabla
l(\mathbf{C})=-\frac{1}{2}\bm{\Sigma}_{x}^{*}\odot(\widetilde{\bm{\Sigma}}_{y}-\bm{\Sigma}_{x}\odot\mathbf{C}).$
(17)

###  V.1 Resulting solver for $\mathcal{P}_{1}~{}$

With the penalty term $\operatorname{\Omega}_{1}(\mathbf{M})$, the proximal
operator is given by

$\operatorname{\textbf{prox}}_{\gamma\nu\operatorname{\Omega}_{1}}(\mathbf{Z})=\sum_{m=1}^{M}(s_{m}-\gamma\nu)_{+}\mathbf{u}_{m}\mathbf{u}_{m}^{\mathsf{H}}$
(18)

where $(.)_{+}$ is the operator that clips the negative input values to $0$,
and $s_{m}$ (resp. $\mathbf{u}_{m}$) are the eigenvalues (resp. eigenvectors)
of the matrix $\mathbf{Z}\in\mathcal{S}_{M}^{+}$. In the original form, the
operator is defined by singular values and vectors, rather than the eigen-
decomposition. However, the mapping is applied on
$\mathbf{C}^{[k]}\in\mathcal{S}_{M}^{+}$, and in this case singular value and
eigenvalue decompositions are equal. The resulting solver is closely related
to existing low-rank approximation methods for matrix completion – see the
singular value thresholding [10], or [32]. The summary is described in the
algorithm 1 for $i=1$.

The choice of parameter $\nu$ is crucial because it controls the rank sparsity
of the estimation $\widetilde{\mathbf{C}}$. First, it is possible to derive
the upper bound value $\nu_{\text{max}}$ beyond which the rank is null, i.e.
$\widetilde{\mathbf{C}}=\mathbf{0}_{M\times M}$. At convergence – for
$k\rightarrow\infty$ – we have the equality
$\mathbf{C}^{[\infty]}=\operatorname{\textbf{prox}}_{\gamma\nu_{\text{max}}\operatorname{\Omega}_{i}}\left(\mathbf{C}^{[\infty]}-\gamma\nabla
l(\mathbf{C}^{[\infty]})\right)$. Solving the equation for
$\mathbf{C}^{[\infty]}=\mathbf{0}_{M\times M}$ gives

$\nu_{\text{max}}=\frac{1}{2}\|\bm{\Sigma}_{x}^{*}\odot\widetilde{\bm{\Sigma}}_{y}\|_{2}$
(19)

where $\|\mathbf{M}\|_{2}$ denotes the maximal singular value of $\mathbf{M}$.
The authors suggest creating a regularization path, run the solver first with
$\nu=\nu_{\text{max}}-1\cdot 10^{-6}$, then decrease the value of $\nu$
iteratively, and finally choose the smallest one which preserves
$\widetilde{\mathbf{C}}$ with rank one. By doing so, the biasing effect of the
penalty term on the estimation is reduced. We emphasize that this suggestion
is empirical, and that the optimal choice of $\nu$ still remains a general
open problem.

###  V.2 Resulting solver for $\mathcal{P}_{0}~{}$

The penalty term given by the indicator function (11) turns the proximal
mapping in the Euclidean projection onto the set $\mathcal{C}$ of $(M\times
M)$ rank-1 matrices:

$\displaystyle\operatorname{\textbf{prox}}_{a\operatorname{\Omega}_{0}}(\mathbf{Z})$
$\displaystyle=\operatornamewithlimits{argmin}_{\mathbf{M}}\frac{1}{2}\|\mathbf{M}-\mathbf{Z}\|^{2}_{\mathsf{F}}+aI_{\mathcal{C}}(\mathbf{M})$
(20)
$\displaystyle=\operatornamewithlimits{argmin}_{\mathbf{M}\in\mathcal{C}}\frac{1}{2}\|\mathbf{M}-\mathbf{Z}\|^{2}_{\mathsf{F}}.$

If $\mathbf{Z}$ has the EVD
$\mathbf{Z}=\sum_{m=1}^{M}s_{m}\mathbf{u}_{m}\mathbf{u}_{m}^{\mathsf{H}}$, the
Eckhart–Young–Mirsky theorem [22, ch. 2.4] proves that

$\operatorname{\textbf{prox}}_{a\operatorname{\Omega}_{0}}(\mathbf{Z})=s_{1}\mathbf{u}_{1}\mathbf{u}_{1}^{\mathsf{H}}.$
(21)

Note that the operator is independent of $a$, thus from $\nu$, because the
indicator function takes either $0$ or $\infty$ values. For the notation
thought, we keep it in the subscript in order to preserve the generalized
notation from the definition 15. The resulting solver is in algorithm 1, with
$i=0$.

###  V.3 Gradient step and stopping criterion

The step $\gamma$ controls the speed of the gradient descent. The stable
convergence is determined by the interval
$\gamma\in]0,\frac{1}{\mathcal{L}}]$, where $\mathcal{L}$ is the Lipschitz
constant of the gradient $\nabla l$. Here it can be easily shown that
$\displaystyle\mathcal{L}=\frac{1}{2}\max_{m,n}(|[\bm{\Sigma}_{x}]_{mn}|^{2})$.
If the goal is to maximize the convergence speed without falling into
instability, $\gamma$ should be the highest one in this interval. With the
diffuse field model (2) it boils down to $\gamma=2$.

For the convergence of the algorithm, we rely on the evolution of
$\mathbf{C}^{[k]}$ and arbitrarily choose
$\|\mathbf{C}^{[k]}-\mathbf{C}^{[k-1]}\|_{\mathsf{F}}\leq 10^{-6}$. Note that
the convergence rate of low-rank approximation with non-convex penalties
remains an open problem (here the solver for $\mathcal{P}_{0}~{}$); but
interestingly, empirical studies evidence their faster convergence than when
using convex ones [48].

Algorithm 1 sensor calibration in ambient noise by proximal gradient descent

input :

$\bm{\alpha}^{[0]}$, $\gamma$, $\nu$, $i$

initialize
$\mathbf{C}^{[0]}=\bm{\alpha}^{[0]}{\bm{\alpha}^{[0]}}^{\mathsf{H}}$, $k=0$
repeat

$\mathbf{C}^{[k+1]}=\operatorname{\textbf{prox}}_{\gamma\nu\operatorname{\Omega}_{i}}\left(\mathbf{C}^{[k]}-\gamma\nabla
l(\mathbf{C}^{[k]})\right)$ $k\leftarrow k+1$

until  _$\|\mathbf{C}^{[k+1]}-\mathbf{C}^{[k]}\|^{2}_{\mathsf{F}}
<10^{-6}\mathrm{~{}or~{}}k=k_{\text{max}}$_;

obtain the EVD
$\mathbf{C}^{[k]}=\sum_{k=1}^{M}s_{i}\mathbf{u}_{i}\mathbf{u}_{i}^{\mathsf{H}}$
and take $\tilde{\bm{\alpha}}=\mathbf{u}_{1}$ [optional: compute $\eta$ for
relative to absolute scaling]  output : $\tilde{\bm{\alpha}}$ [optional:
$\tilde{\bm{\alpha}}_{\text{sc}}=\eta\tilde{\bm{\alpha}}$]

##  VI Numerical analysis of the performance

The aim of this section is to first study the performance and the robustness
of the proposed gain estimators. Concerning performance, we analyze how the
signal-to-noise ratio (SNR) and the number of snapshots in the sample
covariance matrix affect the calibration results. With linear array
geometries, we will also study how the performance depends on the wavelength
$\lambda$. Concerning robustness, we will evaluate the calibration sensitivity
to mismatch with the exact diffuse noise model. Three methods will be
compared: the proximal gradient descent algorithm 1 for
$\mathcal{P}_{1}~{}$and $\mathcal{P}_{0}~{}$(referred to as
pgd-$\mathcal{P}_{1}~{}$and pgd-$\mathcal{P}_{0}~{}$respectively), and the
vanilla estimator from Sec. III (referred to as vanilla).

The simulation generates snapshot measurements $\mathbf{y}_{l}$, and the field
vectors $\mathbf{x}_{l}$ result in the propagation of uncorrelated plane waves
uniformly distributed on the sphere. The SNR of snapshots respects the ratio

$\text{SNR}=10\log_{10}\left(\dfrac{\operatorname{\mathrm{E}}\\{\|\operatorname{diag}(\bm{\alpha})\mathbf{y}_{l}\|^{2}\\}}{\operatorname{\mathrm{E}}\\{\|{\mathbf{w}}_{l}\|^{2}\\}}\right),$
(22)

the ambient noise variance is fixed at $\sigma_{x}^{2}=1$, and the sensor
noise variance $\sigma_{w}^{2}$ is tuned to match the SNR according to (22).
The gains are scaled with respect to the ground truth (see Sec. IV.3,
paragraph c). Results are provided either by the root-mean-square error
$\text{RMSE}=\|\bm{\alpha}_{0}-\tilde{\bm{\alpha}}_{\text{sc}}\|_{2}$, or by
the relative RMSE:
$\text{rRMSE}=\|\bm{\alpha}_{0}-\tilde{\bm{\alpha}}_{\text{sc}}\|_{2}/\|\bm{\alpha}_{0}\|_{2}$.
In practice, the rRMSE is interesting because it is simple to convert a
relative error to a measured level accuracy or a frequency response
uncertainty in dB, as it is usually provided in a microphone datasheet. For
example, a $10\%$ relative error corresponds to $\pm 0.9$ dB.

Note that in the current particular case where $\bm{\alpha}_{0}$ is the ground
truth, it is easy to show that

$\text{rRMSE}=\sqrt{1-\rho^{2}}$ (23)

where
$\rho=|\bm{\alpha}_{0}^{\mathsf{H}}\tilde{\bm{\alpha}}|/(\|\bm{\alpha}_{0}\|\|\tilde{\bm{\alpha}}\|)$
is the correlation coefficient between the estimated gains and the ground
truth. So here the rRMSE is bounded between $0$ and $1$, and the worst (resp.
best) estimator is when rRMSE = $1$ (resp. rRMSE = $0$) with high probability,
i.e. with a little standard deviation of this rRMSE.

###  VI.1 Effect of SNR and snapshot number

In each simulation, the array is planar, forms a $1\text{\,}\mathrm{m}$
diameter disk in which $M=50$ sensors are randomly positioned. The complex
gains $[\bm{\alpha}]_{m}$ are one-mean, the real and imaginary parts uniformly
distributed in the intervals $[0.5,1.5]$ and $[-0.5,0.5]$ respectively. We
consider an acoustic scenario with
$c_{0}=$343\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$ at frequency
$500\text{\,}\mathrm{Hz}$, so $\lambda=$68.6\text{\,}\mathrm{cm}$$.
Statistical results provide the mean and standard deviations of rRMSEs from
$500$ simulations.

Figure1a.pdf(a) Figure1b.pdf(b)

Figure 1: (color online) Complex gain calibration of $M=50$ sensors in ambient
noise. rRMSE means and standard deviations from 500 simulations, with (a)
$L=500$ snapshots, (b) SNR $=0$ dB. Solvers: $\mathcal{P}_{1}~{}$(blue),
$\mathcal{P}_{0}~{}$(green) and vanilla (red).

The parametric study on SNR if given in Fig. 1, for $L=500$ snapshots. The SNR
ranges from $-30$ to $30\text{\,}\mathrm{dB}$ and includes the asymptotic case
$L=\infty$ – cf. horizontal bars on the right. The trends for
pgd-$\mathcal{P}_{1}~{}$and pgd-$\mathcal{P}_{0}~{}$are close, with a slightly
better rRMSE for pgd-$\mathcal{P}_{0}~{}$in high SNR. The vanilla approach
provides irrelevant estimations for all SNRs in the given configuration. For
the parametric study on snapshots plotted in Fig. 1, $L$ spans from 1 to
10000, and $\text{SNR}=$0\text{\,}\mathrm{dB}$$. Again, the vanilla estimation
fails, with a slight decrease for large $L$. Methods
pgd-$\mathcal{P}_{1}~{}$and pgd-$\mathcal{P}_{0}~{}$also keep close trends,
but with a large number of snapshots pgd-$\mathcal{P}_{1}~{}$keeps decreasing
the rRMSE and performs better.

The two figures reveal the potential instability of the vanilla estimation, as
discussed in Sec. III. The proposed algorithm is stable (i.e. with little
standard deviation). Despite the difficulty of solving a non-convex problem,
pgd-$\mathcal{P}_{0}~{}$shows surprisingly good results, although the two
parametric studies evidence marginal differences between
pgd-$\mathcal{P}_{0}~{}$and pgd-$\mathcal{P}_{1}~{}$except in extreme regimes
– for SNR or $L$ large.

###  VI.2 The importance of spatial over-sampling

The diffuse noise model expression evidences that the array geometry can
affect the least-square loss. Indeed, the term
$\bm{\Sigma}_{x}\odot\mathbf{C}$ weights the unknown matrix by a
$\operatorname{sinc}$ function, taking zero values (when $d_{mn}=k\lambda/2$
for $k\geq 1$) and quickly decaying as $d_{mn}/\lambda$ grows.

Figure2a.pdf(a) Figure2b.jpg(b)

Figure 2: (color online) Analysis of the diffuse field spatial sampling. (a)
Complex gain calibration of uniform and random linear arrays with $M=21$,
aperture $1\text{\,}\mathrm{m}$. ULA inter-element spacing
$d=$5\text{\,}\mathrm{cm}$$. Mean rRMSE on 500 simulations for each point,
with pgd-$\mathcal{P}_{1}~{}$. (b) Illustration in the 1-dimensional case of
the 3 regimes: over-sampling (right), critical sampling (middle), and sub-
sampling (left) of the ambient field covariance with a 8-element ULA (top) and
RLA (bottom).

For illustration purpose only, Fig. 2 depicts a 1-dimensional scenario with a
uniform linear array (ULA) and a random linear array (RLA). It draws the
continuous cross-spectral density on the line between two points of coordinate
$r_{x}$ and $r_{y}$, and has the amplitude
$\operatorname{sinc}(2\pi|r_{x}-r_{y}|/\lambda)$. The black squares represent
the spatial sampling obtained by a 8-element ULA (top) and RLA (bottom), so
that the resulting $8\times 8$ matrix $\bm{\Sigma}_{x}$ consists of the values
indicated in these squares. The uniform sampling (top row) clearly reveals 3
distinct regimes, depending on whether the spatial Shannon-Nyquist criterion
is met or not [33]. Let $d$ be the ULA inter-element spacing, then over-
sampling (resp. sub-sampling) is met when $d/\lambda<1/2$ (resp.
$d/\lambda>1/2$). Critical sampling occurs when the off-diagonal black squares
are located to zeros. In this case, we have $\bm{\Sigma}_{x}=\mathbf{I}_{M}$.
This situation also occurs when $\lambda=2d/k$ for $k\geq 1$.

Figures 2 shows the consequence of spatial sampling as a function of
wavelength, with a $21$-element ULA and RLA of aperture $1\text{\,}\mathrm{m}$
(aperture is defined by $\max_{m,n}(d_{mn})$), SNR = $15\text{\,}\mathrm{dB}$,
and $L=1000$. The averaged rRMSE with the solver for $\mathcal{P}_{1}~{}$are
plotted as a function of the normalized wavelength. When
$\bm{\Sigma}_{x}=\mathbf{I}_{M}$ (indicated with orange vertical lines) the
estimation fails for the ULA. Indeed, the loss as initially defined in Eq. (8)
is reduced to

$\begin{split}l(\bm{\alpha}\bm{\alpha}^{\mathsf{H}})&=\frac{1}{2}\|\widetilde{\bm{\Sigma}}_{y}-\mathbf{I}_{M}\odot\bm{\alpha}\bm{\alpha}^{\mathsf{H}}\|^{2}_{\mathsf{F}}\\\
&=\frac{1}{2}\sum_{m=1}^{M}\left([\widetilde{\bm{\Sigma}}_{y}]_{mm}-|[\bm{\alpha}]_{m}|^{2}\right)^{2},\end{split}$
(24)

and since covariance are in $\mathcal{S}_{M}^{+}$ the diagonal terms are real-
valued, therefore the estimation suffers from phase ambiguities on the gains.
The problem does not arise with RLA, because $\bm{\Sigma}_{x}$ is never
diagonal – see Fig. 2. However, for both ULA and RLA the performance decays in
sub-sampling because $[\bm{\Sigma}_{x}]_{mn}$ ($m\neq n$) fades to $0$ as
$\lambda$ decreases. As a result, the ideal configuration to obtain low rRMSE
remains in this so-called over-sampling regime for any array geometry, which
arises when the main lobe of the spatial sinc function contains more than one
sampling point.

###  VI.3 Model mismatch: a case study

Figure3a.pdf0.4(a) Figure3b.pdf(b)

Figure 3: (color online) Sensor calibration in presence of model mismatch: (a)
retina absorbs plane waves in a solid angle $\theta_{c}$. Illustration with
$M=4$. (b) rRMSE means and standard deviations on 500 simulations, with SNR
$=\infty$ and $L=500$.

We want to evaluate the sensitivity of the calibration due to the model
mismatch from the diffuse field property. To do so, we perturb the uniform
spatial distribution of plane waves by adding an anechoic retina in the far
field, see the geometric configuration in Fig. 3 (4 sensors are drawn only for
the visibility of the illustration). The size of the retina is adjusted by the
solid angle $\theta_{c}\in[0,\pi]$, the case $\theta_{c}=\pi$ depicting the
scenario where the array is nearby an anechoic wall. Such a partially diffuse
field was studied analytically in [8], but we consider this mismatch as
unknown here.

Figure 3 compares the rRMSE means and standard deviations for the 3 solvers,
plus the case in which no calibration is performed in gray. There is no sensor
noise (i.e. $\sigma_{w}^{2}=0$), and we keep the ULA configuration from the
previous section and the wavelength is reset to
$\lambda=$68.6\text{\,}\mathrm{cm}$$ – in the spatial over-sampling regime. As
depicted in Fig. 3, the retina is oriented $\pi/4\text{\,}\mathrm{rad}$ from
the ULA principal axis to discard symmetries. In this configuration, the
vanilla approach works better than in Fig. 1 but remains more sensitive to
model mismatch, with a more significant rRMSE increase than with proximal
solvers. As expected, it remains less stable, with larger standard deviations.
The solver pgd-$\mathcal{P}_{1}~{}$is marginally better than
pgd-$\mathcal{P}_{0}~{}$when $\theta_{c}$ is large, and both do not suffer
from the mismatch when $\theta_{c}\leq\pi/4$. The comparison with uncalibrated
gains (gray) gives an indication of the range of angles $\theta_{c}$ where the
calibration brings an improvement or not; nevertheless it should not be taken
as a generality because the rRMSE level in gray is only representative of the
gain distribution that was chosen – cf. Sec. VI.1.

##  VII Experimental results

This second section of results evaluates the ability of gain calibration on
acoustic real data captured by a “Simcenter Sound Camera” antenna from
Siemens. The scenario is the following: a 43-microphone spiral array is
$30\text{\,}\mathrm{cm}$ wide, and records acoustic pressures in a reverberant
room at the sampling frequency of $51.2\text{\,}\mathrm{kHz}$. We emphasize
that the microphone array is already calibrated off-the-shelf by Siemens, so
we deliberately assign gain coefficients on the measurements in the presented
results below. The picture of the experiment is provided in Fig. 4. First the
quality of the diffuse field is investigated from the raw measurements, then
the calibration step is performed and assessed.

Figure 4: (color online) Experimental setup: a 43-microphone array in a
reverberant room.

###  VII.1 Ambient noise conditions

Figure5a.pdf0.3(a)
Figure5b.jpg0.3(b)
Figure5c.jpg0.3(c)
Figure5d.jpg0.3(d) ,

Figure 5: (color online) Sample covariance $\widetilde{\bm{\Sigma}}_{y}$ off-
diagonal coefficients before calibration in the reverberant room (a)
$\widetilde{\bm{\Sigma}}_{y}$ off-diagonal matrix coefficients at
$f=$2.5\text{\,}\mathrm{k}\mathrm{Hz}$$; black circles indicate values paired
with the $32$-nd microphone. (b) Spatial covariance with respect to the
$32$-nd sensor; analytical value (background) and sampled values (in circles).
Real part of (c) the analytical scaled covariance with
$d_{mn}\in$[0,0.27]\text{\,}\mathrm{m}$$ , and of (d) the
$\widetilde{\bm{\Gamma}}_{y}$ coefficients, in the (frequency, distance)
plane.

The field is known to be diffuse in reverberant rooms above the Schroeder
frequency [38, 35]. Here, the room has a volume of
$411\text{\,}{\mathrm{m}}^{3}$, a reverberation time of about
$9\text{\,}\mathrm{s}$, and a Schroeder frequency of about
$0.3\text{\,}\mathrm{k}\mathrm{Hz}$. For frequency to wavelength conversion,
we measured the sound speed
$c_{0}=$341.5\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$. The diffuse
sound field is excited by two speakers positioned at room corners and emitting
white noise. From a $42\text{\,}\mathrm{s}$-long recording, the sample
covariance matrix $\widetilde{\bm{\Sigma}}_{y}$ is computed with $L=6008$
snapshots. To do so, the DFT is applied to time series blocks of $512$ samples
with frequency resolution $\Delta f=$100\text{\,}\mathrm{Hz}$$, thereby the
condition (5) is respected.

The raw measurements before calibration are visualized in Fig. 5. First, Fig.
5 shows the real and imaginary parts of $[\widetilde{\bm{\Sigma}}_{y}]_{mn}$
coefficients at $f=$2.5\text{\,}\mathrm{k}\mathrm{Hz}$$. Since
$\widetilde{\bm{\Sigma}}_{y}$ is hermitian, and since the diagonal elements
involve variances, the coefficients from the upper triangular part are plotted
only, that is for $n>m$. For comparison with the model, the analytical
expression is plotted in green, with an approximation of $\sigma_{x}^{2}$
based on the average of the measured variances (located in the diagonal of
$\widetilde{\bm{\Sigma}}_{y}$). We recall that the analytical imaginary part
is null. The measured covariance points are well concentrated and close to the
model. The dispersion of these points around their average trend are either
due to the dispersion of the sensor gains, or to the mismatch between the
model and the sensed sound field.

Figure 6: (color online) Evaluating the over-sampling with the experimental
microphone array: number of pairwise distances $d_{mn}$ such that
$d_{mn}/\lambda<1/2$, $\lambda=c_{0}/f$.

In the 2-dimensional case with irregular geometry, determining over-sampling
regime described in Sec. VI.2 is not trivial. Yet,

  * •

a first clue to help answer this question is given in Fig. 5, still for
$f=$2.5\text{\,}\mathrm{k}\mathrm{Hz}$$. The background plots the analytical
covariance between the $32$-nd microphone position and a second coordinate
within a $$30\text{\,}\mathrm{cm}$\times$30\text{\,}\mathrm{cm}$$ square. The
circled points are filled with the real part measured values, that is
$\Re([\widetilde{\bm{\Sigma}}_{y}]_{mn})$ for $m=32$ and $n\neq m$. Note that
these values correspond to the ones circled in black in Fig. 5. The over-
sampling was identified with ULAs when $\lambda/(2d)$. The visual
interpretation of this in Fig. 5 is that nearby sensors are present in the
main lobe of the sinc-shaped spatial covariance. This is true when the
covariance is centered at the $32$-nd sensor, and seems to hold if centered at
other microphone positions. Nevertheless, this visual interpretation is valid
for $f=$2.5\text{\,}\mathrm{k}\mathrm{Hz}$$ only;

  * •

to generalize over other frequencies, a second clue is provided in Fig. 6. It
draws, as a function of $f=c_{0}/\lambda$, the number of pairwise distances
respecting the condition $d_{mn}/\lambda>1/2$. This number quickly decays from
$f=$0.8\text{\,}\mathrm{k}\mathrm{Hz}$$ to $3\text{\,}\mathrm{k}\mathrm{Hz}$,
and reaches 0 at $7.7\text{\,}\mathrm{k}\mathrm{Hz}$. As a rule of thumb, at
least $10\text{\,}\mathrm{\char 37\relax}$ of the distances respect below
$f=$3.4\text{\,}\mathrm{k}\mathrm{Hz}$$.

Finally, Figs. 5 and 5 provide a complete side-by-side comparison of the model
with all the measurements in the $(f,d_{mn})$ plane. Since $M=42$ there are
$M(M-1)/2=861$ pairwise distances; with $79$ frequency bins from $0.1$ to
$7.9\text{\,}\mathrm{k}\mathrm{Hz}$, Fig. 5 consists of $68019$ points
overall. Rather than directly relying on the covariance, we plot the real part
of the coherence. By definition, the sample coherence matrix
$\widetilde{\bm{\Gamma}}_{y}$ is a normalization of the sample covariance
coefficients $[\widetilde{\bm{\Sigma}}_{y}]_{mn}$, that is

$[\widetilde{\bm{\Gamma}}_{y}]_{mn}=\dfrac{[\widetilde{\bm{\Sigma}}_{y}]_{mn}}{\sqrt{[\widetilde{\bm{\Sigma}}_{y}]_{mm}[\widetilde{\bm{\Sigma}}_{y}]_{nn}}}.$
(25)

Although the speakers generates white noise their frequency response is not
flat, thus the source variance $\sigma_{x}^{2}$ varies across frequencies.
Microphone responses are not flat as well. This is why plotting the coherence
is more convenient: it reveals the ambient noise spatial correlation by
mitigating the hardware frequency responses and the variations of
$\sigma_{x}^{2}$ over frequencies. Note that microphone amplitude part of the
gains are discarded, but not the phase part. The side-by-side comparison of
the real part of analytical and measured coherences reveal a good similarity
in the considered frequency range.

In summary, this section validates that the diffuse noise model is a reliable
choice to perform the proposed sensor gain calibration, in the acoustic
scenario with a reverberant room. Although the array geometry is irregular and
planar, a proxy was proposed to identify an approximate frequency range where
spatial over-sampling occurs, that is below $3.5\text{\,}\mathrm{kHz}$.

###  VII.2 Gain calibration results

Figure7a.jpg(a) Figure7b.jpg(b)

Figure 7: (color online) Experimental sensor gain calibration: (a) amplitude
gains $|[\tilde{\bm{\alpha}}_{\text{sc}}]_{m}|$, (b) relative error for each
sensor and frequency.

In this section, we perform the microphone gain calibration from $0.1$ to
$7.9\text{\,}\mathrm{k}\mathrm{Hz}$. Since the sensing model is expressed for
a given frequency, the calibration is launched for each frequency bin,
resulting in $79$ vectors $\tilde{\bm{\alpha}}$, or $79\times 43$ coefficients
overall. The scaling is done by choosing $\bm{\alpha}_{0}$ as the ground
truth, which is known since gains were artificially applied on the snapshots
readings.

The absolute value of these coefficients are drawn in Fig. 7, ground truth
$\bm{\alpha}_{0}$ included. Additionally, Fig. 7 shows the coefficient-wise
relative errors with respect to the ground truth
$|[\tilde{\bm{\alpha}}_{\text{sc}}]_{m}-[\bm{\alpha}_{0}]_{m}|/|[\bm{\alpha}_{0}]_{m}|$.
First, the vanilla approach is the least reliable one, because most estimated
gains are irrelevant, except in low frequencies and around
$2.5\text{\,}\mathrm{k}\mathrm{Hz}$. Concerning low frequencies, the
estimations are correct because the condition (6) is respected (cf. Fig. 6)
and stabilizes the element-wise division by
$\bm{\Sigma}_{x}\approx\mathbf{1}_{M\times M}$. However, the proposed
algorithm performs better on a wider range of frequencies:

  * •

with pgd-$\mathcal{P}_{1}~{}$the estimation works until
$5\text{\,}\mathrm{kHz}$; and fails beyond. Some errors are visible from $3$
to $5\text{\,}\mathrm{kHz}$ on the peripheral microphones, associated to large
sensor indexes;

  * •

the estimations by pgd-$\mathcal{P}_{0}~{}$are good on the whole studied
frequency range, with a homogeneous relative error per sensor and per
frequency. A local increase remains visible in the sector where
$f\in$[3,5]\text{\,}\mathrm{kHz}$$ and where sensor index $m\geq 36$.

Figure 8: (color online) Experimental microphone gain calibration RMSE: total
(top), amplitude only (middle), and phase only (bottom).

Finally, the performance per frequency is given in Fig. 8. In order to
discriminate the nature of errors, the RMSE on raw complex gains (top), their
absolute values (middle) and their phase (bottom) are plotted. The gray graphs
indicate the RMSE if no calibration was performed. The RMSE are the same for
the three compared methods below $600\text{\,}\mathrm{Hz}$, and quickly
diverge beyond. Interestingly, the error in phase are identical for
pgd-$\mathcal{P}_{0}~{}$and pgd-$\mathcal{P}_{1}~{}$, until
pgd-$\mathcal{P}_{1}~{}$fails from $5\text{\,}\mathrm{kHz}$. Thus, the
difference between both is mainly explained by the amplitude component. Again,
pgd-$\mathcal{P}_{0}~{}$reveals to be the most stable: the total RMSE is small
in low frequencies, slightly increases from $2$ to $3.8\text{\,}\mathrm{kHz}$,
and remains stable beyond. Note that this trend is more obvious on the phase
RMSE. The authors suggests explaining this trend by the spatial sampling
density, as it was discussed previously. Indeed, the proxy indicator Fig. 6
reveals that over-sampling is marginal beyond $3.5\text{\,}\mathrm{kHz}$. This
frequency interval matches the one in which RMSE is lower. Thus, this
experimental analysis confirms again that the relative dimension between the
array geometry and the wavelength drives the calibration performance, even
with a relevant model choice.

##  VIII Conclusion

The complex gain calibration of sensors from a synchronous array in ambient
noise is possible by leveraging the spatial structure of the cross-spectral
measurements. Unlike the conventional choice of snapshot readings, it relies
on second-order statistical measurements: the problem is cast in the least-
square form to fit the diffuse noise model with the sample covariance matrix.
The problem is non-convex, and can be modified by convex relaxation via the
low-rank matrix approximation framework.

A proximal algorithm was proposed, in order to deal efficiently with $M$
large, and to solve the problem with or without convex relaxation. The overall
results globally show that the non-convex optimization provides better results
in practice, although it does not guaranty obtaining the global minimizer.
Finally, a general caveat was identified: the spatial sampling should be
sufficiently dense, otherwise the covariance matrix is such that its off-
diagonal elements have too small amplitudes to obtain a proper estimation of
the gain phases.

##  IX Acknowledgements

We thank Johan Clavier and Thomas Rougier for their quality work in the
realization of experiments. We also thank Cédric Herzet for the fruitful
discussions and suggested directions concerning proximal methods. This study
has been produced in the framework of LUG2 supported by Région Auvergne Rhône-
Alpes and BPIFrance (FUI22). It was performed within the framework of the
Labex CeLyA of Université de Lyon, operated by the French National Research
Agency (ANR-10-LABX-0060/ANR-11-IDEX-0007).

## References

  * Akbar _et al._ [2020] Akbar, N., Dickins, G., Thomas, M. R., Samarasinghe, P., and Abhayapala, T. (2020). “A novel method for obtaining diffuse field measurements for microphone calibration,” in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , IEEE, pp. 226–230.
  * Akbar _et al._ [2021] Akbar, N., Dickins, G., Thomas, M. R. P., Samarasinghe, P., and Abhayapala, T. (2021). “Reducing modal error propagation through correcting mismatched microphone gains using rapid,” IEEE, Toronto, ON, Canada, pp. 811–814, 10.1109/ICASSP39728.2021.9415030.
  * Balzano and Nowak [2007] Balzano, L., and Nowak, R. (2007). “Blind calibration of sensor networks,” IEEE, Cambridge, MA, USA, pp. 79–88, 10.1109/IPSN.2007.4379667.
  * Balzano and Nowak [2008] Balzano, L., and Nowak, R. (2008). “Blind Calibration of Networks of Sensors: Theory and Algorithms,” (Springer US, Boston, MA), pp. 9–37, 10.1007/978-0-387-68845-9_1.
  * Barcelo-Ordinas _et al._ [2019] Barcelo-Ordinas, J. M., Doudou, M., Garcia-Vidal, J., and Badache, N. (2019). “Self-calibration methods for uncontrolled environments in sensor networks: A reference survey,” Ad Hoc Networks 88, 142–159, https://www.sciencedirect.com/science/article/pii/S1570870518306115, 10.1016/j.adhoc.2019.01.008.
  * Beck [2017] Beck, A. (2017). _First-order methods in optimization_ , 25 (SIAM).
  * Bilen _et al._ [2014] Bilen, C., Puy, G., Gribonval, R., and Daudet, L. (2014). “Convex optimization approaches for blind sensor calibration using sparsity,” Signal Processing, IEEE Transactions on 62(18), 4847–4856.
  * Blake and Waterhouse [1977] Blake, W. K., and Waterhouse, R. V. (1977). “The use of cross-spectral density measurements in partially reverberant sound fields,” Journal of Sound and Vibration 54(4), 589–599.
  * Boyd _et al._ [2011] Boyd, S., Parikh, N., Chu, E., Peleato, B., Eckstein, J. _et al._ (2011). “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Foundations and Trends® in Machine learning 3(1), 1–122.
  * Cai _et al._ [2010] Cai, J.-F., Candès, E. J., and Shen, Z. (2010). “A Singular Value Thresholding Algorithm for Matrix Completion,” SIAM Journal on Optimization 20(4), 1956–1982, 10.1137/080738970.
  * Chen [2013] Chen, H. (2013). “Robustness analysis of nearfield subband beamformers in the presence of microphone gain and phase errors,” Digital Signal Processing 23(5), 1712–1719.
  * Cui _et al._ [2017] Cui, H., Peng, W., Liu, T., and Wei, X. (2017). “Cramér-rao bounds for multiple near field sources location under unknown gain/phase response,” in _2017 29th Chinese Control And Decision Conference (CCDC)_ , IEEE, pp. 5169–5173.
  * Davenport and Romberg [2016] Davenport, M. A., and Romberg, J. (2016). “An overview of low-rank matrix recovery from incomplete observations,” IEEE Journal of Selected Topics in Signal Processing 10(4), 608–622, 10.1109/JSTSP.2016.2539100.
  * Delaine _et al._ [2019] Delaine, F., Lebental, B., and Rivano, H. (2019). “In Situ Calibration Algorithms for Environmental Sensor Networks: A Review,” IEEE Sensors Journal 19(15), 5968–5978, 10.1109/JSEN.2019.2910317.
  * Druet _et al._ [2019] Druet, T., Recoquillay, A., Chapuis, B., and Moulin, E. (2019). “Passive guided wave tomography for structural health monitoring,” The Journal of the Acoustical Society of America 146(4), 2395–2403, 10.1121/1.5128332.
  * Fazel _et al._ [2001] Fazel, M., Hindi, H., and Boyd, S. P. (2001). “A rank minimization heuristic with application to minimum order system approximation,” in _Proceedings of the 2001 American Control Conference.(Cat. No. 01CH37148)_ , IEEE, Vol. 6, pp. 4734–4739.
  * Frederiksen [2013] Frederiksen, E. (2013). “Acoustic metrology–an overview of calibration methods and their uncertainties,” International journal of metrology and quality engineering 4(2), 97–107.
  * Friedlander and Weiss [1993] Friedlander, B., and Weiss, A. J. (1993). “Performance of direction-finding systems with sensor gain and phase uncertainties,” Circuits Syst. Signal Process. 12(1), 3–35, 10.1007/BF01183145.
  * Garnier and Papanicolaou [2016] Garnier, J., and Papanicolaou, G. (2016). _Passive imaging with ambient noise_ (Cambridge University Press).
  * Gilquin _et al._ [2019] Gilquin, L., Bouley, S., Antoni, J., Le Magueresse, T., and Marteau, C. (2019). “Sensitivity analysis of two inverse methods: Conventional beamforming and bayesian focusing,” Journal of Sound and Vibration 455, 188–202, https://www.sciencedirect.com/science/article/pii/S0022460X19302664, https://doi.org/10.1016/j.jsv.2019.05.002.
  * Godin _et al._ [2010] Godin, O. A., Zabotin, N. A., and Goncharov, V. V. (2010). “Ocean tomography with acoustic daylight,” Geophysical Research Letters 37(13).
  * Golub and Van Loan [2013] Golub, G. H., and Van Loan, C. F. (2013). _Matrix computations_ (JHU press).
  * Gower [2010] Gower, J. C. (2010). “Procrustes methods,” WIREs Comp Stat 2(4), 503–508, 10.1002/wics.107.
  * Grant and Boyd [2014] Grant, M., and Boyd, S. (2014). “CVX: Matlab software for disciplined convex programming, version 2.1” http://cvxr.com/cvx.
  * Gribonval _et al._ [2012] Gribonval, R., Chardon, G., and Daudet, L. (2012). “Blind calibration for compressed sensing by convex optimization,” in _Proc. Speech and Signal Processing (ICASSP) 2012 IEEE Int. Conf. Acoustics_ , pp. 2713–2716, 10.1109/ICASSP.2012.6288477.
  * Hill [1995] Hill, D. A. (1995). “Spatial correlation function for fields in a reverberation chamber,” IEEE Transactions on Electromagnetic Compatibility 37(1), 138, 10.1109/15.350256.
  * Jacobsen [1979] Jacobsen, F. (1979). “The diffuse sound field,” Technical Report Report No. 27.
  * Jaggi [2013] Jaggi, M. (2013). “Revisiting frank-wolfe: Projection-free sparse convex optimization,” in _International Conference on Machine Learning_ , PMLR, pp. 427–435.
  * Li _et al._ [2019] Li, F., Yang, X., Zhang, Y., Luo, W., and Gan, W. (2019). “Passive ocean acoustic tomography in shallow water,” The Journal of the Acoustical Society of America 145(5), 2823–2830.
  * Ling and Strohmer [2015] Ling, S., and Strohmer, T. (2015). “Self-calibration and biconvex compressive sensing,” Inverse Problems 31(11), 115002, 10.1088/0266-5611/31/11/115002.
  * Ling and Strohmer [2018] Ling, S., and Strohmer, T. (2018). “Self-Calibration and Bilinear Inverse Problems via Linear Least Squares,” SIAM Journal on Imaging Sciences 11(1), 252–292, 10.1137/16M1103634.
  * Ma _et al._ [2011] Ma, S., Goldfarb, D., and Chen, L. (2011). “Fixed point and Bregman iterative methods for matrix rank minimization,” Mathematical Programming 128(1), 321–353, 10.1007/s10107-009-0306-5.
  * Menon _et al._ [2012a] Menon, R., Gerstoft, P., and Hodgkiss, W. S. (2012a). “Asymptotic eigenvalue density of noise covariance matrices,” IEEE Transactions on Signal Processing 60(7), 3415–3424, 10.1109/TSP.2012.2193573.
  * Menon _et al._ [2012b] Menon, R., Gerstoft, P., and Hodgkiss, W. S. (2012b). “Cross-correlations of diffuse noise in an ocean environment using eigenvalue based statistical inference,” The Journal of the Acoustical Society of America 132(5), 3213–3224.
  * Nélisse and Nicolas [1997] Nélisse, H., and Nicolas, J. (1997). “Characterization of a diffuse field in a reverberant room,” The Journal of the Acoustical Society of America 101(6), 3517–3524, 10.1121/1.418313.
  * Parikh _et al._ [2014] Parikh, N., Boyd, S. _et al._ (2014). “Proximal algorithms,” Foundations and Trends® in Optimization 1(3), 127–239.
  * Pavlis and Vernon [1994] Pavlis, G. L., and Vernon, F. L. (1994). “Calibration of seismometers using ground noise,” Bulletin of the seismological society of America 84(4), 1243–1255.
  * Pierce [2019] Pierce, A. D. (2019). _Acoustics: an introduction to its physical principles and applications_ (Springer).
  * Rossing and Rossing [2014] Rossing, T. D., and Rossing, T. D. (2014). _Springer handbook of acoustics_ (Springer).
  * Roux _et al._ [2005] Roux, P., Sabra, K. G., Kuperman, W. A., and Roux, A. (2005). “Ambient noise cross correlation in free space: Theoretical approach,” The Journal of the Acoustical Society of America 117(1), 79–84, 10.1121/1.1830673.
  * Seydoux _et al._ [2016] Seydoux, L., Shapiro, N. M., de Rosny, J., Brenguier, F., and Landès, M. (2016). “Detecting seismic activity with a covariance matrix analysis of data recorded on seismic arrays,” Geophysical Journal International 204(3), 1430–1442.
  * Shapiro _et al._ [2005] Shapiro, N. M., Campillo, M., Stehly, L., and Ritzwoller, M. H. (2005). “High-resolution surface-wave tomography from ambient seismic noise,” Science 307(5715), 1615–1618.
  * Siderius _et al._ [2010] Siderius, M., Song, H., Gerstoft, P., Hodgkiss, W. S., Hursky, P., and Harrison, C. (2010). “Adaptive passive fathometer processing,” The Journal of the Acoustical Society of America 127(4), 2193–2200.
  * Weaver and Lobkis [2001] Weaver, R. L., and Lobkis, O. I. (2001). “Ultrasonics without a source: Thermal fluctuation correlations at MHz frequencies,” Physical Review Letters 87(13), 134301\.
  * Wei _et al._ [2020] Wei, K., Aviles-Rivero, A., Liang, J., Fu, Y., Schönlieb, C.-B., and Huang, H. (2020). “Tuning-free plug-and-play proximal algorithm for inverse imaging problems,” in _International Conference on Machine Learning_ , PMLR, pp. 10158–10169.
  * Weiss and Friedlander [1990] Weiss, A. J., and Friedlander, B. (1990). “Eigenstructure methods for direction finding with sensor gain and phase uncertainties,” Circuits, Systems and Signal Processing 9(3), 271–300, 10.1007/BF01201215.
  * Woolfe _et al._ [2015] Woolfe, K. F., Lani, S., Sabra, K. G., and Kuperman, W. A. (2015). “Monitoring deep-ocean temperatures using acoustic ambient noise,” Geophysical Research Letters 42(8), 2878–2884.
  * Yao _et al._ [2019] Yao, Q., Kwok, J. T., Wang, T., and Liu, T.-Y. (2019). “Large-scale low-rank matrix learning with nonconvex regularizers,” IEEE Transactions on Pattern Analysis and Machine Intelligence 41, 2628–2643, 10.1109/TPAMI.2018.2858249.
  * Yardim _et al._ [2014] Yardim, C., Gerstoft, P., Hodgkiss, W. S., and Traer, J. (2014). “Compressive geoacoustic inversion using ambient noise,” The Journal of the Acoustical Society of America 135(3), 1245–1255.


# Evaluating Zero-Shot Scoring for In Vitro Antibody Binding Prediction with
Experimental Validation

Divya Nori
Massachusetts Institute of Technology
Department of Electrical Engineering and Computer Science
77 Massachusetts Ave, Cambridge, MA 02139
<EMAIL_ADDRESS>
Simon V. Mathis11footnotemark: 1
University of Cambridge
Department of Computer Science and Technology
<EMAIL_ADDRESS>
Amir Shanehsazzadeh
Absci Corporation
<EMAIL_ADDRESS>
Work done during internship at Absci.

###### Abstract

The success of therapeutic antibodies relies on their ability to selectively
bind antigens. AI-based antibody design protocols have shown promise in
generating epitope-specific designs. Many of these protocols use an _inverse
folding_ step to generate diverse sequences given a backbone structure. Due to
prohibitive screening costs, it is key to identify candidate sequences likely
to bind in vitro. Here, we compare the efficacy of $8$ common scoring
paradigms based on open-source models to classify antibody designs as binders
or non-binders. We evaluate these approaches on a novel surface plasmon
resonance (SPR) dataset, spanning $5$ antigens. Our results show that existing
methods struggle to detect binders, and performance is highly variable across
antigens. We find that metrics computed on flexibly docked antibody-antigen
complexes are more robust, and ensembles scores are more consistent than
individual metrics. We provide experimental insight to analyze current scoring
techniques, highlighting that the development of robust, zero-shot filters is
an important research gap.

##  1 Introduction

Antibodies are a prominent class of therapeutic agents, primarily because they
can selectively bind to a wide array of disease-causing target antigens. This
binding specificity is governed by interactions between the antibody’s
paratope region and the antigen’s epitope region [1]. Epitope specificity is
primarily facilitated by the structure of hypervariable loops in the paratope
called complimentarity determining regions (CDRs).

Computationally designing the CDRs such that they bind to specific epitopes
with high affinity has become a problem of great interest. Researchers often
start with a known or predicted antibody backbone structure, and the task is
to design sequences that assume this backbone conformation [2]. Recently,
generative machine learning methods have proven effective for this inverse
folding problem on protein design [3]. Given the hypervariability of CDRs,
there are many sequences that can fold into the same 3D backbone structure.
Filtering the diverse outputs of an inverse folding model in silico to select
for sequences with greater chance of in vitro success remains an important
challenge, since experimental resources are often limited.

In the protein space, many in silico metrics are commonly used to evaluate the
quality of designs [4], [5]. Here, we investigate how these metrics perform
for antibodies. We benchmark several scoring paradigms on their ability to
filter antibody designs by predicted success in a surface plasmon resonance
(SPR) assay, an industry standard assay used to measure binding affinity.

##  2 Background

### Problem set-up

As shown in Figure 1, inverse folding models design candidate sequences given
a desired backbone structure. A subset of these sequences will fold into a 3D
structure in vitro, and a smaller subset then assume paratope residue
conformations that interact strongly with the antigen [1].

Figure 1: Antibody-inverse folding model with antibody-antigen complex used as
input predicts a set of candidate CDR sequences. These sequences are then
filtered at the sequence, undocked structure, or docked structure level, which
is the focus of this study. The HCDR3 sequences are from [6].

Here, we assess whether filtering at the sequence level, undocked structure
level, or docked structure level is effective. Particularly, we compare $8$
different scoring paradigms on the binder vs non-binder classification task.
The candidate scoring paradigms we evaluate are all zero-shot methods –
meaning they are not trained on any binding affinity data which is often
available in only limited quantities.

### Related work

In general protein binder design, in silico metrics have been used to filter
libraries and improve binder success rates. Most recently, it was shown that
monomer metrics are less effective in predicting binders than complex metrics
[7]. On antibodies, prior work assesses Rosetta’s predicted free energy of
binding ($\Delta\Delta G$) as a binder scoring function and achieves an AUROC
of $0.55$ on a library designed for HER2 [8]. Recently, a diffusion model was
trained to assess the quality of antibody designs by generating and scoring
docked poses [9]. Language models have also been used to perform efficient
evolution on antibodies, indicating that their likelihoods can be used to
select predicted binders [10]. We assess whether such approaches are
effective, given novel wet-lab insight.

##  3 Methods

### Data

To generate the antibody libraries scored in this work, an antibody-specific
inverse folding model was given an antibody-antigen complex backbone
structure, antibody framework sequences, and antigen sequence as input. The
model was then used to design several candidate antibody sequences,
specifically all three HCDRs of a reference antibody to the antigen of
interest. Libraries were designed for $5$ different antigens, and wet-lab data
was generated by SPR. A future publication will describe additional details on
data generation. See Table 2 for details on library size and % binders.

### Candidate scoring functions

Table 1 describes the $8$ candidate scoring functions we evaluate on the
binder classification task. We hypothesize that antibody candidate sequences
that contain common or predictable motifs are more likely to be binders.
Therefore, we benchmark pseudo-perplexity predicted by the ESM-2 650M
parameter model as a sequence scoring metric [11]. Pseudo-perplexity was
computed for each sequence, scaled by min-max normalization, flipped to range
[-1,0], and translated to range [0,1]. All scoring paradigms going forward are
normalized in a similar manner. Subsequently, we hypothesize that antibody
candidates with biologically plausible structural qualities are more likely to
be successful binders. As an undocked structure scoring metric, we benchmark
the residue-level confidence scores given by ABodyBuilder2 (ABB2), an
antibody-specific folding model [12].

Other discriminative scoring functions used in-practice are self-consistency
metrics, such as similarity between the folded antibody candidate and the
input backbone structure [5]. Here, we benchmark root mean square deviation
(RMSD) as a similarity measure, under three Kabsch alignments [13].

Given that flexible interactions between antibody and antigen is critical for
binding, we benchmark metrics on flexibly docked complexes. To generate these
complexes, we apply a previously proposed architecture, dyMEAN [14], and make
several modifications to the training procedure which are described in the
Appendix. We evaluate three scoring paradigms on the predicted complexes,
selected because they capture established heuristics on antibody-antigen
interactions. The scores we evaluate measure interface proximity, contact
proximity, and interface interactions (details in Appendix).

|  Model/Method |  Metric
---|---|---
Sequence Only |  ESM |  Pseudo-Perplexity
Undocked Structure |  ABodyBuilder2 |  Residue-Level Model Confidence
|  Antibody-Aligned Antibody |  RMSD
|  All HCDR-Aligned Antibody |  RMSD
|  HCDR3-Aligned Antibody |  RMSD
Docked Complex |  Modified dyMEAN |  Interface Proximity
|  Modified dyMEAN |  Contact Proximity
|  Modified dyMEAN |  Interface Interactions
Table 1: Evaluation of $8$ candidate scoring functions on the binder
classification task across $3$ classes of metrics: sequence only, undocked
structure, and docked complex.

##  4 Results

We assess the performance of each scoring paradigm in the classification
setting, evaluating across multiple antigens. Figure 2 shows each scoring
function’s distribution of AUROC values. We observe that antibody-aligned RMSD
performs the best (mean of $0.65$ and standard deviation of $0.15$),
indicating that undocked structure filtering may be effective. However, all
three of the RMSD-based undocked structure metrics display high variance in
performance. Interface proximity is the most consistent metric (mean $0.51$
and standard deviation $0.028$), followed by contact proximity and interface
interactions. This reveals that docked complex metrics exhibit less
variability.

Figure 2: Distribution of AUROC values for various scoring functions across
$5$ antigens.

After evaluating the performance of each metric individually, we assess
whether we can achieve better discrimination with a combined scoring paradigm.
Since the scoring functions are not strongly correlated, as shown in the
Appendix, we hypothesize that a combined method may boost performance.

For any linear combination of metrics, we fit coefficients on $4$ out of $5$
antigens, where our objective is mean AUROC. We do not consider combinations
whose standard deviation exceeds $0.1$. To avoid leakage, we evaluate the
best-performing ensemble score on the held-out antigen. We call this approach
“Best Combination on Hold-Out,” and report retrieval precision on the task of
filtering to $10\%$ of the original antibody design library’s size in Figure
3. Retrieval precision, or Precision @ $K$, measures the proportion of true
binders in the top $K$ highest-scoring candidates. Retrieval precision and
recall curves for the individual scoring metrics are given in the Appendix. We
also evaluate performance on an oracle approach where we find the best
ensemble score on all $5$ antigens, which we refer to “Best Combination on
All.” Our baseline is given by the binder percentage in each of the libraries,
and performance of the best individual metric (antibody-aligned RMSD) is also
shown for comparison.

We observe that ensembling does not improve performance. However, the standard
deviation of the hold-out ensemble scores is $0.074$, compared to $0.22$ for
antibody-aligned RMSD. This indicates that ensembling is more robust to
antigen variation. Additionally, the hold-out ensemble scores are just $2.9$
percentage points below the oracle for retrieval precision, revealing that we
cannot perform better on this classification task using these scoring
criteria. This underscores the need for better prediction methods beyond those
evaluated here. We find that in almost all folds, interface proximity has the
highest weight. The specific weights are given in the Appendix. This reveals
that flexible docking information is critical for this task, and better
antibody-antigen docking models can push the needle.

Figure 3: Distribution of Precision @ $10\%$ for ensemble scores and baselines
across $5$ antigens.

##  5 Conclusions

This study provides experimental insight to show that there is no one-size-
fits-all approach for antibody assessment, and current open-source tools are
insufficient to reliably filter diverse antibody libraries. We learn that
antibody-aligned RMSD outperforms other individual scoring paradigms. However,
sequence and undocked structure-only scoring paradigms struggle to handle
antigen variation, and metrics that use antigen information are more robust.
We find that linear combinations of the individual scoring paradigms exhibit
lower variance in performance but leave much to be desired. Improved flexible
docking models and further experimental benchmark efforts are critical to
identify promising antibody candidates in silico and accelerate their clinical
development.

## Acknowledgments

The authors wish to thank Kieran Didi, Amaro Taylor-Weiner, Christine Lemke,
Byron Olsen, and Kristin Iannacone for critical review of this manuscript.

## Competing interest statement

The authors are current or former employees, contractors, interns, or
executives of Absci Corporation and may hold shares in Absci Corporation.

## References

  * Peng et al. [2014] Hung-Pin Peng, Kuo Hao Lee, Jhih-Wei Jian, and An-Suei Yang.  Origins of specificity and affinity in antibody–protein interactions.  _Proceedings of the National Academy of Sciences_ , 111(26):E2656–E2665, 2014.
  * Hsu et al. [2022] Chloe Hsu, Robert Verkuil, Jason Liu, Zeming Lin, Brian Hie, Tom Sercu, Adam Lerer, and Alexander Rives.  Learning inverse folding from millions of predicted structures.  pages 8946–8970, 2022.
  * Dauparas et al. [2022] J. Dauparas, I. Anishchenko, N. Bennett, H. Bai, R. J. Ragotte, L. F. Milles, B. I. M. Wicky, A. Courbet, R. J. de Haas, N. Bethel, P. J. Y. Leung, T. F. Huddy, S. Pellock, D. Tischer, F. Chan, B. Koepnick, H. Nguyen, A. Kang, B. Sankaran, A. K. Bera, N. P. King, and D. Baker.  Robust deep learning–based protein sequence design using ProteinMPNN.  _Science_ , 378(6615):49–56, October 2022.  doi: 10.1126/science.add2187.  URL https://doi.org/10.1126/science.add2187.
  * Johnson et al. [2023] Sean R Johnson, Xiaozhi Fu, Sandra Viknander, Clara Goldin, Sarah Monaco, Aleksej Zelezniak, and Kevin K Yang.  Computational scoring and experimental evaluation of enzymes generated by neural networks.  _bioRxiv_ , pages 2023–03, 2023.
  * Wu et al. [2022] Kevin Eric Wu, Kevin K Yang, Rianne van den Berg, James Zou, Alex Xijie Lu, and Ava P Amini.  Protein structure generation via folding diffusion.  2022\.
  * Shanehsazzadeh et al. [2023] Amir Shanehsazzadeh, Sharrol Bachas, Matt McPartlon, George Kasun, John M Sutton, Andrea K Steiger, Richard Shuai, Christa Kohnert, Goran Rakocevic, Jahir M Gutierrez, et al.  Unlocking de novo antibody design with generative artificial intelligence.  _bioRxiv_ , pages 2023–01, 2023.
  * Bennett et al. [2023] Nathaniel R Bennett, Brian Coventry, Inna Goreshnik, Buwei Huang, Aza Allen, Dionne Vafeados, Ying Po Peng, Justas Dauparas, Minkyung Baek, Lance Stewart, et al.  Improving de novo protein binder design with deep learning.  _Nature Communications_ , 14(1):2625, 2023.
  * Mason et al. [2019] Derek M Mason, Simon Friedensohn, Cédric R Weber, Christian Jordi, Bastian Wagner, Simon Meng, Pablo Gainza, Bruno E Correia, and Sai T Reddy.  Deep learning enables therapeutic antibody optimization in mammalian cells by deciphering high-dimensional protein sequence space.  _BioRxiv_ , page 617860, 2019.
  * Peng et al. [2023] Zhangzhi Peng, Chenchen Han, Xiaohan Wang, Dapeng Li, and Fajiie Yuan.  Generative diffusion models for antibody design, docking, and optimization.  _bioRxiv_ , 2023.  doi: 10.1101/2023.09.25.559190.
  * Hie et al. [2023] Brian L Hie, Varun R Shanker, Duo Xu, Theodora UJ Bruun, Payton A Weidenbacher, Shaogeng Tang, Wesley Wu, John E Pak, and Peter S Kim.  Efficient evolution of human antibodies from general protein language models.  _Nature Biotechnology_ , 2023.
  * Lin et al. [2022] Zeming Lin, Halil Akin, Roshan Rao, Brian Hie, Zhongkai Zhu, Wenting Lu, Allan dos Santos Costa, Maryam Fazel-Zarandi, Tom Sercu, Sal Candido, et al.  Language models of protein sequences at the scale of evolution enable accurate structure prediction.  _BioRxiv_ , 2022:500902, 2022.
  * Abanades et al. [2023] Brennan Abanades, Wing Ki Wong, Fergus Boyles, Guy Georges, Alexander Bujotzek, and Charlotte M Deane.  Immunebuilder: Deep-learning models for predicting the structures of immune proteins.  _Communications Biology_ , 6(1):575, 2023.
  * Kabsch [1976] Wolfgang Kabsch.  A solution for the best rotation to relate two sets of vectors.  _Acta Crystallographica Section A: Crystal Physics, Diffraction, Theoretical and General Crystallography_ , 32(5):922–923, 1976.
  * Kong et al. [2023] Xiangzhe Kong, Wenbing Huang, and Yang Liu.  End-to-end full-atom antibody design.  _ICML Conference Poster_ , 2023.
  * Dunbar et al. [2014] James Dunbar, Konrad Krawczyk, Jinwoo Leem, Terry Baker, Angelika Fuchs, Guy Georges, Jiye Shi, and Charlotte M Deane.  Sabdab: the structural antibody database.  _Nucleic acids research_ , 42(D1):D1140–D1146, 2014.
  * Eastman et al. [2017] Peter Eastman, Jason Swails, John D Chodera, Robert T McGibbon, Yutong Zhao, Kyle A Beauchamp, Lee-Ping Wang, Andrew C Simmonett, Matthew P Harrigan, Chaya D Stern, et al.  Openmm 7: Rapid development of high performance algorithms for molecular dynamics.  _PLoS computational biology_ , 13(7):e1005659, 2017.

##  Appendix A Dataset details

Table 2 shows the number of screened antibodies and the percentage of binding
antibodies for each library.

| Antigen #1 | Antigen #2 | Antigen #3 | Antigen #4 | Antigen #5
---|---|---|---|---|---
# of Sequences | $64$ | $56$ | $68$ | $68$ | $58$
# of Binders | $9$ | $15$ | $10$ | $6$ | $14$
% Binders | $14.1$ | $26.8$ | $14.7$ | $8.8$ | $24.1$
Table 2: Five libraries, each designed and screened against one of five
antigens using surface plasmon resonance (SPR), are used as experimental data
in this study. Each library consists of between 56 and 68 unique sequences
with between 8.8% and 24.1% of the sequences binding to their target antigen.

##  Appendix B dyMEAN training

We make the following modifications to the training procedure of dyMEAN. While
we train on the Structural Antibody Database [15] as in the original dyMEAN
paper, we split SAbDab into training and validation sets according to clusters
formed by complexes sharing $\geq 40\%$ antigen sequence similarity. The
clusters containing any of the $5$ evaluation antigens were dropped from both
training and validation. Additionally, we use OpenMM [16] to relax the bound
antigen structures before model input to ensure that we are not leaking
information about docked position. Finally, we provide the model with relaxed
ABB2-predicted antibody structures to focus on learning flexible docking
rather than folding.

##  Appendix C Docked complex scoring paradigms

We benchmark three scoring functions computed on docked complexes predicted by
dyMEAN. The first is an interface proximity score which gives the average
distance between the residues on the antibody interface and their nearest
corresponding antigen residues. Here, we define interface as the contiguous
region between the first residue of HCDR1 to the last residue of HCDR3,
including HFWR2 and HFWR3. Residue-residue distances are computed as minimum
distance between the side chain atoms of both residues. We also compute a
contact proximity score which computes the same score on just the antibody
contact residues which are defined as having at least one side chain atom
within $5$ Angstroms of an antigen atom.

Lastly, we compute an interface interactions score which evaluates the types
of interactions at the antibody-antigen interface. We begin by identifying all
antibody-antigen residue pairs where an antibody residue atom is within at
most 3 Angstroms of an antigen residue atom. We then classify each residue-
residue pair, awarding higher scores to polar-polar interactions, nonpolar-
nonpolar Van der Waals interactions, positive-negative ionic interactions, and
aromatic stacking interactions.

##  Appendix D Correlation analysis

Figure 4 shows correlations between the $8$ scoring paradigms. Correlation
coefficients were computed on all pairs of each antigen’s set of scores, and
these coefficients were averaged across all $5$ targets. As expected, the
RMSD-based undocked structure metrics exhibit some correlation, and the
metrics computed on flexibly docked complexes are also correlated. With the
exception of ESM pseudo-perplexity and ABB2 pRMSD, which are weakly
correlated, we do not observe much correlation across the sequence, undocked
structure, and docked complex metric groups. Therefore, we hypothesize that
combining individual metrics may reveal new information useful for binder
classification.

Figure 4: Correlations between the $8$ scoring paradigms, averaged across
antigens, show that each metric is capturing fairly different signal.
Combining metrics could yield favorable performance.

##  Appendix E Retrieval precision and recall

Retrieval precision and recall metrics are particularly relevant because we
can measure the number of true binders amongst the $K$ we may select to test
in the lab, while ensuring that a high proportion of binders are not missed.
Precision @ $K$ measures the proportion of true binders in the top $K$
highest-scoring candidates, and Recall @ $K$ measures the proportion of
relevant items in the top $K$. We evaluate retrieval precision and recall
values for each scoring function averaged across the $5$ antigens, as shown in
Figure 5. We focus on the task of filtering to $10\%$ of the original antibody
design library’s size. We observe that antibody-aligned RMSD has the highest
Precision @ $10\%$ with a value of $0.26$, averaged across antigens, compared
to a baseline mean binder rate of $0.13$. Antibody-aligned RMSD also has the
highest retrieval recall of $0.16$ at this threshold.

Figure 5: The retrieval precision and recall curves, averaged across antigens,
show how effectively the different scoring paradigms can select $K$ designs.
At a threshold of $10\%$, which is approximately $K=6$, antibody-aligned RMSD
displays the best retrieval precision and recall.

##  Appendix F Combined scoring function weights

In Table 3, we report the coefficients found when each group of $4$ antigens
was used for optimization.

Held-out Antigen | Metric 1 | Metric 2 | Metric 3
---|---|---|---
Antigen #1 | Interface Proximity, $0.7$ | Contact Proximity, $0.2$ | Interface Interactions, $0.1$
Antigen #2 | Antibody-Aligned RMSD, $0.5$ | Interface Interactions, $0.4$ | ESM Pseudo-Perplexity, $0.1$
Antigen #3 | Interface Proximity, $0.7$ | Contact Proximity, $0.2$ | Interface Interactions, $0.1$
Antigen #4 | Interface Proximity, $0.7$ | Contact Proximity, $0.2$ | Interface Interactions, $0.1$
Antigen #5 | Interface Proximity, $0.7$ | ESM Pseudo-Perplexity, $0.2$ | Contact Proximity, $0.1$
Oracle | Interface Proximity, $0.6$ | Contact Proximity, $0.3$ | HCDR3-Aligned RMSD, $0.1$
Table 3: Best weighted ensemble scores for each of the $5$ folds, as measured
by mean AUROC.