# Supervised Metric Learning for Retrieval
via Contextual Similarity Optimization

Christopher Liao
Boston University
<EMAIL_ADDRESS>&Theodoros Tsiligkaridis
MIT Lincoln Laboratory
<EMAIL_ADDRESS>&Brian Kulis
Boston University
<EMAIL_ADDRESS>
###### Abstract

Existing deep metric learning approaches fall into three general categories:
contrastive learning, average precision (AP) maximization, and classification.
We propose a novel alternative approach, _contextual similarity optimization_
, inspired by work in unsupervised metric learning. Contextual similarity is a
discrete similarity measure based on relationships between neighborhood sets,
and is widely used in the unsupervised setting as pseudo-supervision. Inspired
by this success, we propose a framework which optimizes _a combination of
contextual and cosine similarities_. Contextual similarity calculation
involves several non-differentiable operations, including the heaviside
function and intersection of sets. We show how to circumvent non-
differentiability to explicitly optimize contextual similarity, and we further
incorporate appropriate similarity regularization to yield our novel metric
learning loss. The resulting loss function achieves state-of-the-art Recall @
1 accuracy on standard supervised image retrieval benchmarks when combined
with the standard contrastive loss. Code is released here:
https://github.com/Chris210634/metric-learning-using-contextual-similarity

##  1 Introduction

Supervised metric learning aims to learn a transformation from data to an
embedding space, where similar samples are close together and dis-similar
samples are far apart. Three types of metric learning approaches currently
exist: contrastive learning, average precision (AP) maximization, and
classification. Contrastive losses roughly follow the idea proposed by Hadsell
et al. (2006): minimize the distance between positive pairs and maximize the
distance between negative pairs. Alternatively, approaches such as ProxyNCA
(Movshovitz-Attias et al. (2017)) use a classification loss to optimize a
linear layer on top of the embedding. The linear layer and the embedding are
trained jointly, and the former is discarded. More recently, there is strong
interest in average precision (AP) maximization. AP is a metric which takes
into account the rank of samples relative to each other and is non-
differentiable. Several studies tackle this non-differentiability issue: Fast-
AP (Cakir et al. (2019)) uses a soft-binning approach; Smooth-AP (Brown et al.
(2020)) and ROADMAP (Ramzi et al. (2021)) use sigmoid-like functions to
approximate the heaviside function used for ranking. In the current work, we
take a different approach to metric learning, using the notion of contextual
similarity from unsupervised metric learning (Zhong et al. (2017) and Kim et
al. (2022)).

Unsupervised metric learning aims to learn a good representation without using
labels. Most existing work in this area focus on a semi-supervised pseudo-
labeling approach. In particular, a teacher model is used to guess the true
similarity between unlabeled samples; the guess is then used as artificial
supervision for training a student model. Instead of directly using the cosine
similarity between teacher features as the pseudo-similarity label, most
state-of-the-art methods rely on a weighted combination of the cosine
similarity and contextual similarity (Ge et al. (2020a), Ge et al. (2020b),
Kim et al. (2022)). Inspired by this success, we endeavor to explicitly
optimize contextual similarity.

The definition of contextual similarity varies. In the person re-ID
literature, most papers use the $k$-reciprocal re-rank distance proposed in
Zhong et al. (2017). Contextual similarity is based on the following
intuition: (1) neighbors are likely to have the same label regardless of
absolute similarity, (2) two samples are likely to have the same label if they
share the same neighbors, (3) reciprocal relationships are important: sample
$i$ is likely to be similar to sample $j$ if $i$ is a neighbor of $j$ and $j$
is a neighbor of $i$, (4) neighbors of neighbors can be treated as neighbors.
These observations lead directly to a rigorous definition of contextual
similarity presented in Section 3.1. Existing supervised metric learning
approaches already implicitly minimize some notion of contextual similarity –
this is why contextual similarity works as a pseudo-label. We hypothesize that
_explicitly_ optimizing contextual similarity boosts recall performance when
combined with standard contrastive learning.

We use the term “context” to refer to the set of neighbors around a sample.
Roughly speaking, the contextual similarity between samples $i$ and $j$ is a
fraction between 0 and 1 indicating the degree of overlap between the contexts
of $i$ and $j$. We seek to minimize the contextual similarity between samples
with different labels, and vice versa. Clearly, this requires the contextual
similarity to be differentiable with respect to the embeddings. Contextual
similarity is non-differentiable because it involves manipulation of indicator
functions, e.g. to calculate the intersection between contexts. We find that
existing tools to circumvent non-differentiability, such as using sigmoid
approximations, are insufficient for our problem. Therefore, we develop a
novel optimization approach for contextual similarity in Section 3.2.

We make two technical and one experimental contribution to the metric learning
literature:

  1. 1.

We derive a highly non-trivial differentiable contextual loss function.

  2. 2.

We propose a simple but novel similarity regularizer which significantly
improves the standard contrastive loss and may be of independent interest.

  3. 3.

Our framework improves the Recall @ 1 performance significantly over the
current state-of-the-art across four diverse image retrieval benchmarks.

##  2 Literature Review

Supervised Metric Learning  There are three main approaches: (1) contrastive
(2) classification and (3) AP maximization. The contrastive loss was first
proposed by Hadsell et al. (2006). The triplet loss (Weinberger et al. (2005))
operates on triplets (a pair of positive samples with a negative sample)
instead of pairs and is competitive with the contrastive loss when combined
with an effective triplet sampling strategy (Wu et al. (2017)). While
contrastive and triplet losses are well-known, they remain the go-to method
for metric learning and Musgrave et al. (2020a) show that they are comparable
in performance to many recent methods. Variations of the contrastive loss
include the lifted structure loss (Oh Song et al. (2016)) and multi-similarity
loss (MS) (Wang et al. (2019)). Standard classification methods can achieve
state-of-the-art metric learning results when appropriate tricks are used
(Boudiaf et al. (2020), Zhai & Wu (2018), Teh et al. (2020)). AP maximization
methods learn to rank samples within a mini-batch correctly. The AP objective
is non-differentiable because of the heaviside used for ranking. As a
workaround, Fast-AP uses soft-binning; Smooth-AP uses a low-temperature
sigmoid; ROADMAP uses an upper bound on the heaviside instead of an
approximation. Our work is distinct from all three categories while having
similarities to learning to rank (see Figure 2).

Unsupervised Metric Learning  There is extensive interest in unsupervised and
domain adaptive metric learning, mainly in the context of person re-ID (see
survey Ye et al. (2021)). Contextual similarity is a cornerstone for most
recent work in this area. Most unsupervised person re-ID papers use the
$k$-reciprocal re-rank distance, which is a weighted combination of Euclidean
distance and Jaccard distance between reciprocal-neighbor sets. Their
formulation has three disadvantages: (1) it operates on the entire dataset
instead of a batch; (2) it involves a series of set expansion operations which
are hard to vectorize, much less make differentiable; (3) it is optimized for
the person re-ID application. More recently, STML (Kim et al. (2022)) propose
an unsupervised metric learning framework for image retrieval which uses a
simpler batch-wise contextual similarity measure. We mostly follow STML’s
contextual similarity definition, with some minor tweaks. We emphasize that
all work cited above use contextual similarity as a pseudo-label for
unsupervised learning, while our work optimizes the contextual similarity in
the supervised setting.

##  3 Method

This section is organized as follows: (1) We state the mathematical definition
of contextual similarity, which loosely follows Kim et al. (2022), then (2) we
explain how to optimize the objective.

Figure 1: Method Overview. Our final framework consists of three loss terms.
$\mathcal{L}_{\text{context}}$ pulls apart the contexts of dis-similar samples
and pushes together the contexts of similar samples. This is illustrated on
the left, where black triangles indicate samples with any label. See also
analysis in Section 4.2. $\mathcal{L}_{\text{reg}}$ is a regularizer that
encourages utilization of the entire embedding space by regularizing average
similarities to a fixed value. $\mathcal{L}_{\text{contrast}}$ is the standard
contrastive loss.

Notation Denote the normalized output of the embedding network as
$f_{i}\in\mathbb{R}^{d}$. $s_{ij}\in[-1,1]$ denotes the cosine similarity
between the samples $i$ and $j$ ($s_{ij}=\langle f_{i},f_{j}\rangle$).
$y_{ij}\in\\{0,1\\}$ denotes the true similarity between $i$ and $j$.
$y_{ij}=1$ if samples $i$ and $j$ share the same label; otherwise, $y_{ij}=0$.
We use uppercase letters to denote matrices, math script to denote sets, and
lowercase letters to denote scalars. $i$, $j$ and $p$ are reserved for sample
indices.

###  3.1 Contextual Similarity Definition

We compute the contextual similarity on a single batch, not the entire
dataset. Two samples are contextually similar if they share the same
neighbors, i.e. the intersection of their $k$-neighbor sets is large.
Following this intuition, we calculate $\tilde{w}_{ij}$, the preliminary
contextual similarity between samples $i$ and $j$:

$\mathcal{N}_{k+\epsilon}(i)=\\{j\>|\>s_{ij}\leq s_{ip}+\epsilon\text{ where
$p$ denotes the $k$-th closest neighbor of $i$}\\}$ (1)
$\tilde{w}_{ij}=\begin{cases}|\mathcal{N}_{k+\epsilon}(i)\cap\mathcal{N}_{k+\epsilon}(j)|\>/\>|\mathcal{N}_{k+\epsilon}(i)|&,\>\text{if
}j\in\mathcal{N}_{k+\epsilon}(i)\\\ 0&,\text{ otherwise.}\end{cases}$ (2)

We include $i$ in $\mathcal{N}_{k+\epsilon}(i)$ (so if $k=2$, then
$\mathcal{N}_{k+0}(i)$ includes two elements: $i$ and its closest neighbor).
Then, we use query expansion and symmeterize to obtain the final contextual
similarity. Query expansion (Arandjelović & Zisserman (2012)) is a standard
evaluation-time trick in metric learning. It boosts retrieval performance by
retrieving neighbors of a sample’s neighbors. Analogously, we adjust the
contextual similarity by averaging $\tilde{w}_{ij}$ over the set of close
reciprocal neighbors $\mathcal{R}_{k/2}(i)$.

$\mathcal{N}_{k/2+\epsilon}(i)=\\{j\>|\>s_{ij}\leq s_{ip}+\epsilon\text{ where
$p$ denotes the $k/2$-th closest neighbor of $i$}\\}$ (3)
$\mathcal{R}_{k/2+\epsilon}(i)=\\{j\>|\>j\in\mathcal{N}_{k/2+\epsilon}(i)\text{
and }i\in\mathcal{N}_{k/2+\epsilon}(j)\\}$ (4)
$\hat{w}_{ij}=\frac{1}{|\mathcal{R}_{k/2+\epsilon}(i)|}\sum_{p\in\mathcal{R}_{k/2+\epsilon}(i)}\tilde{w}_{pj}\>\>\>,\>\>\>w_{ij}=\frac{1}{2}(\hat{w}_{ij}+\hat{w}_{ji}).$
(5)

Note that $w_{ij}\in[0,1]$ and only depend on the cosine similarities $s_{ij}$
between embeddings. We want to optimize the embeddings such that $w_{ij}$
converges to $y_{ij}$. The value of $k$ is not arbitrary; it must be set to
the number of samples per label in the mini-batch. For example, a standard
metric learning setup is to randomly sample 32 labels and then sample 4 images
per label. In this case, we set $k=4$.

###  3.2 Optimization

The definition in the previous sub-section clearly contains three discrete
operations: (1) greater-than, (2) logical-and, and (3) intersection. For
optimization, we will deal extensively with indicator matrices: we denote as
$\mathbbm{1}_{\mathcal{N}}\in\mathbb{R}^{n\times n}$ the indicator matrix
where $\mathbbm{1}_{\mathcal{N}}(i,j)=1\>\text{if}\>j\in\mathcal{N}(i)$ for
some set $\mathcal{N}$. We use $\odot$ to denote element-wise multiplication.

Greater-than This is used in Eq. 1 and 3 and is equivalent to the non-
differentiable heaviside function $\theta$. Our approach is to use the exact
value of $\theta(\cdot)$ in the forward pass and a constant positive gradient
$\alpha$ in the backward pass. A constant positive gradient is reasonable
since $\theta$ is a (non-strictly) increasing function.

$\text{Forward: }\theta(x)=1\text{ if }x\geq 0\>;\>0\text{ otherwise
}\>\>\>\>\text{Backward: }\frac{\partial\theta(x)}{\partial x}=\alpha.$ (6)

In Section 4.1 we show with a toy experiment that this approach is robust
despite being somewhat heuristic. In contrast, $\theta$ is traditionally
approximated by a sigmoid (e.g. for AP approximation and in Gated Recurrent
Networks Cho et al. (2014) ):

$\theta_{\sigma}(x)=\frac{1}{1+\exp\left(x/\tau\right)}.$ (7)

$\theta_{\sigma}$ trades off the quality of the approximation with the domain
where gradients are non-zero. As temperature $\tau$ decreases,
$\theta_{\sigma}$ approaches $\theta$, but gradient vanishes everywhere except
in a small region around the boundary. This behavior is not intuitive: it is
undesirable to only have a large gradient at the boundary. Some prior work
(e.g. ROADMAP) side-step this issue by using an upper-bound to the heaviside
function (where the right side of the heaviside function increases linearly),
which solves the gradient issue at the expense of grossly over-estimating the
true objective. In our case, this is especially concerning since we will be
multiplying together indicator functions. In Section 5.1, we show that our
approach in Eq. 6 achieves better empirical results than a sigmoid
approximation.

Logical-and  A logical-and is used explicitly in Eq. 4 and implicitly in Eq. 2
and 5. There are two differentiable substitutes for logical-and: $\min$ and
multiplication. Multiplication is smooth, but has no gradient at the origin.
$\min$ is logically consistent on the continuous domain [0,1], but the
gradient is not continuous when inputs are equal. We found experimentally that
multiplication is the best option, see row 5 of Table 3. A gradient of zero at
the origin is desireable in the case of Eq. 4 and 5 (query expansion step). A
work-around for the zero-gradient issue for Eq. 2 is presented in the next
sub-section.

Intersection  The number of elements in the intersection in Eq. 2 is
calculated using matrix multiplication. A sample
$p\in\mathcal{N}_{k+\epsilon}(i)\cap\mathcal{N}_{k+\epsilon}(j)$ if
$p\in\mathcal{N}_{k+\epsilon}(i)\text{ and }p\in\mathcal{N}_{k+\epsilon}(j)$.
Using multiplication for the logical-and, the number of elements in the
intersection can be expressed compactly as a function of the indicator
matrices:

$M_{+}=\mathbbm{1}_{\mathcal{N}_{k+\epsilon}}\mathbbm{1}_{\mathcal{N}_{k+\epsilon}}^{\intercal}\>,\>\text{where}\>M_{+}(i,j)=\left|\mathcal{N}_{k+\epsilon}(i)\cap\mathcal{N}_{k+\epsilon}(j)\right|.$
(8)

We are now ready to state the loss function by combining the definition of
contextual similarity with the optimization method.

Figure 2: Illustration of $\mathcal{N}_{k+\epsilon}$ optimization using the
simplified loss function $\mathcal{L}_{1}$ in Eq. 18. Negative samples which
are in the set $\mathcal{N}_{k+\epsilon}(i)$ are pushed away from $i$ by a
constant gradient. Positive samples which are not in the set
$\mathcal{N}_{k+\epsilon}(i)$ are pulled closer. This behavior is analogous to
the standard contrastive loss: $\mathcal{L}_{\text{contrast}}$ applies a
constant gradient to pairs violating a fixed margin, while $\mathcal{L}_{1}$
applies a constant gradient to pairs violating the flexible margin defined by
the $k$-th neighbor. $k=4$ in the illustration.

###  3.3 Loss

Our loss function straightforwardly combines the previous two sub-sections,
with one exception for the intersection calculation in Eq. 2. Optimizing the
intersection with $M_{+}$ tends to focus on pairs of neighboring samples
because using multiplication for logical-and zeros out the gradient when both
inputs are 0. This is problematic since the resulting embedding becomes
clustered regardless of true similarity (see Figure 4). We mitigate this
problem by optimizing the intersection of the complements
$\mathcal{N}_{k+\epsilon}^{c}(i)\cap\mathcal{N}_{k+\epsilon}^{c}(j)$ in
addition to optimizing the original intersection in Eq. 2. c denotes the
complement of a set. Maximizing the size of the intersection between two sets
is equivalent to maximizing the size of the intersection between their
complements, and vice versa. This can be trivially proven under the
assumptions that the universal set is constant and that the size of both sets
is constant.

Figure 3: This figure justifies the non-standard approach of optimizing
$\mathcal{L}_{1}$ (Eq. 18) by overriding the gradient of $\theta$. We generate
synthetic 2-D data consisting of five concentric circles; each color
represents a label. We use a standard three layer MLP with ReLU activations to
output a 2D embedding. We train on increasingly noisy data and plot the
resulting embeddings on clean data. We only sample 4 points per label per
batch ($k=4$). $\mathcal{L}_{1}$ minimization results in reasonable embeddings
despite its simplicity.

Let $D\in\mathcal{R}^{n\times n}$ denote the matrix where $D(i,j)\in[0,4]$ is
the squared Euclidean distance between the normalized features of sample $i$
and $j$. $D(i,j)=2-2s_{ij}$. sg denotes stop gradient.

Step 1 (Neighborhood Optimization):

$\mathbbm{1}_{\mathcal{N}_{k+\epsilon}}(i,j)=\theta(-D(i,j)+\text{sg}(D(i,p))+\epsilon)\text{
where $p$ denotes the $k$-th closest neighbor of $i$.}$ (9)

Step 2 (Optimizing Intersection of Neighborhoods):

$\begin{split}\tilde{W}=\frac{1}{2}\left(\frac{M_{+}}{\text{sg}(|\mathcal{N}_{k+\epsilon}(i)|)}+\frac{M_{-}}{\text{sg}(|\mathcal{N}_{k+\epsilon}^{c}(i)|)}\right)\odot\mathbbm{1}_{\mathcal{N}_{k+\epsilon}}\>\\\
\text{where}\>\>M_{+}=\mathbbm{1}_{\mathcal{N}_{k+\epsilon}}\mathbbm{1}_{\mathcal{N}_{k+\epsilon}}^{\intercal}\text{
and
}M_{-}=\mathbbm{1}_{\mathcal{N}_{k+\epsilon}^{c}}\mathbbm{1}_{\mathcal{N}_{k+\epsilon}^{c}}^{\intercal}.\end{split}$
(10)

Step 3 (Query Expansion):

$\mathbbm{1}_{\mathcal{N}_{k/2+\epsilon}}(i,j)=\theta(-D(i,j)+\text{sg}(D(i,p))+\epsilon)\text{
where $p$ denotes the $k/2$-th closest neighbor of $i$}$ (11)
$\mathbbm{1}_{\mathcal{R}_{k/2+\epsilon}}=\mathbbm{1}_{\mathcal{N}_{k/2+\epsilon}}\odot\mathbbm{1}_{\mathcal{N}_{k/2+\epsilon}}^{\intercal}$
(12)
$\hat{W}=\frac{\mathbbm{1}_{\mathcal{R}_{k/2+\epsilon}}\tilde{W}}{|\mathcal{R}_{k/2+\epsilon}(i)|}\>,\>W=\frac{1}{2}\left(\hat{W}+\hat{W}^{\intercal}\right).$
(13)

Finally, we use the MSE loss to optimize $w_{ij}:=W(i,j)$ against the true
similarity labels $y_{ij}$:

$\mathcal{L}_{\text{context}}=\frac{1}{n^{2}}\sum_{i,j|i\neq
j}(y_{ij}-w_{ij})^{2}.$ (14)

$\mathcal{L}_{\text{context}}$ is the key to our framework and works
reasonably on its own. However, it has two vulnerabilities: (1) similar to
learning to rank losses, $\mathcal{L}_{\text{context}}$ can be small even if
the distance between positive pairs is large, so long as there are no closer
negative pairs (see Figure 2); (2) $\mathcal{L}_{\text{context}}$ tends to
converge to a solution where the average cosine similarity between all pairs
is large, suggesting that only a small portion of the available embedding
space is utilized. In response to problem (1), we add the standard contrastive
loss (Hadsell et al. (2006)), which explicitly optimizes the cosine similarity
toward fixed margins $\delta_{+}$ and $\delta_{-}$. In response to problem
(2), we add a non-standard but straightforward similarity regularizer, which
regularizes the average cosine similarity between all pairs toward a fixed
value $\tilde{s}$. Our final framework (Eq. 17) minimizes a weighted
combination of the three losses.

$\mathcal{L}_{\text{contrast}}=\frac{\sum_{i,j|y_{ij}=1}(\delta_{+}-s_{ij})_{+}}{\big{|}\\{i,j|y_{ij}=1\text{
and
}\delta_{+}-s_{ij}>0\\}\big{|}}+\frac{\sum_{i,j|y_{ij}=0}(s_{ij}-\delta_{-})_{+}}{\big{|}\\{i,j|y_{ij}=0\text{
and }s_{ij}-\delta_{-}>0\\}\big{|}}$ (15)
$\mathcal{L}_{\text{reg}}=\left(\tilde{s}-\frac{1}{n^{2}}\sum_{i,j}^{n^{2}}s_{ij}\right)^{2}$
(16)
$\mathcal{L}_{\text{ours}}=\lambda\mathcal{L}_{\text{context}}+\gamma\mathcal{L}_{\text{reg}}+(1-\lambda-\gamma)\mathcal{L}_{\text{contrast}}.$
(17)

Explanation of Hyperparameters  $\alpha$ controls the magnitude of the
heaviside gradient. $\epsilon$ is a margin-like parameter for the neighborhood
set. $\delta_{+}$ and $\delta_{-}$ are the positive and negative margins resp.
for the contrastive loss. $\tilde{s}$ is the desired average cosine similarity
between all pairs. $\lambda$ and $\gamma$ control the relative weighting
between the three losses.

##  4 Analysis

The contextual similarity loss function $\mathcal{L}_{\text{context}}$ is
highly non-trivial. In this section we demystify each of the three steps using
a toy experiment and gradient analysis.

###  4.1 Step 1: Neighborhood optimization

Consider the following simplified version of $\mathcal{L}_{\text{context}}$ :

$\mathcal{L}_{1}=\mathcal{L}_{\text{MSE}}(y_{ij},\mathbbm{1}_{\mathcal{N}_{k+\epsilon}}(i,j)).$
(18)

This is a valid loss function that works reasonably on its own (see Figure 3).
During each iteration, we sample a batch with exactly $k$ samples from each
label. Intuitively, $\mathcal{L}_{1}$ reaches a minimum of zero when
$\mathcal{N}_{k+\epsilon}(i)$ includes only $i$ and the $k-1$ other samples
with the same label, i.e. all samples are correctly ranked. Using the gradient
of $\theta(\cdot)$ defined in Eq. 6, we see that negative samples which are in
$\mathcal{N}_{k+\epsilon}(i)$ receive a gradient of magnitude $\alpha$ pushing
it away from $i$, while positive samples which are not in
$\mathcal{N}_{k+\epsilon}(i)$ receive a gradient of magnitude $\alpha$ pushing
it toward $i$. See Figure 2.

Figure 4: This figure justifies the $M_{-}$ term in Eq. 10. The left plot
shows the distribution of distances between pairs in normalized embedding
space when we only optimize $M_{+}$; the right plot shows the same when we
optimize both $M_{-}$ and $M_{+}$. A uniformly distributed embedding space has
an average pair-wise distance of $\sqrt{2}$ because most directions are
orthogonal in high dimensions. This is indicated by the green distribution.
Clearly, more of the embedding space is utilized when we optimize both $M_{-}$
and $M_{+}$.

###  4.2 Step 2: Optimizing Intersection of neighborhoods

While $\lambda\mathcal{L}_{1}+(1-\lambda)\mathcal{L}_{\text{contrastive}}$ is
a reasonable loss function already, row 1 of Table 3 shows that
$\mathcal{L}_{1}$ collapses without the contrastive loss ($\lambda=1$). We
propose that this is because $\mathcal{L}_{1}$ provides very sparse gradients:
$\partial\mathcal{L}_{1}/\partial s_{ij}=0$ for most $(i,j)$ pairs. Consider
the following more complicated loss:

$\mathcal{L}_{2}=\mathcal{L}_{\text{MSE}}(y_{ij},\tilde{W}(i,j)).$ (19)

$\tilde{W}$ is defined in Eq. 10. $\tilde{W}(i,j)$ has two terms multiplied
together element-wise. The $\mathbbm{1}_{\mathcal{N}_{k+\epsilon}}$ term was
already addressed in the previous sub-section, so we focus on the first term,
which optimizes the intersection between neighborhood sets $M_{+}$. We offer a
straight-forward intuition: _Maximizing the intersection between the
neighborhood sets of two samples is equivalent to pushing one sample towards
the context of the other sample and vice versa_ ; _minimizing the intersection
is equivalent to pulling apart the contexts of the two samples_. We show this
intuition by analyzing the gradient.

For simplicity, consider $\epsilon=0$, such that the normalization factors in
Eq. 10 are equal to $k$. Further consider a negative pair of samples $i$ and
$j$ ($y_{ij}=0$), where $j$ is wrongly ranked w.r.t. $i$ (i.e.
$j\in\mathcal{N}_{k}(i)$). Following the loss function equations:

$\begin{split}M_{+}(i,j)&=\langle\mathbbm{1}_{\mathcal{N}_{k}}(i),\mathbbm{1}_{\mathcal{N}_{k}}(j)\rangle\>\>\>\>\>\>\>\>M_{-}(i,j)=\langle\mathbbm{1}_{\mathcal{N}_{k}}^{c}(i),\mathbbm{1}_{\mathcal{N}_{k}}^{c}(j)\rangle\\\
\tilde{W}(i,j)&=\frac{1}{2}\left(\frac{1}{k}M_{+}(i,j)+\frac{1}{k}M_{-}(i,j)\right)\cdot\underbrace{\mathbbm{1}_{\mathcal{N}_{k}}(i,j)}_{=1}.\end{split}$
(20)

$\tilde{W}(i,j)\in(0,1]$ must be a non-zero positive number. Using the
equation for $\mathcal{L}_{2}$, we see that the gradient w.r.t.
$\tilde{W}(i,j)$ must be positive because $y_{ij}=0$. We are now ready for the
backward pass. For simplicity of notation, $\partial\tilde{W}(i,j):=g_{ij}>0$
denotes the gradient of the loss w.r.t. $\tilde{W}(i,j)$.

$\begin{split}\partial M_{+}(i,j)&=\partial M_{-}(i,j)=\frac{1}{2k}g_{ij}\\\
\partial\mathbbm{1}_{\mathcal{N}_{k}}(i)=\frac{1}{2k}g_{ij}\mathbbm{1}_{\mathcal{N}_{k}}(j)\>&,\>\partial\mathbbm{1}_{\mathcal{N}_{k}}(j)=\frac{1}{2k}g_{ij}\mathbbm{1}_{\mathcal{N}_{k}}(i)\\\
\partial\mathbbm{1}_{\mathcal{N}_{k}}^{c}(i)=\frac{1}{2k}g_{ij}\mathbbm{1}_{\mathcal{N}_{k}}^{c}(j)\>&,\>\partial\mathbbm{1}_{\mathcal{N}_{k}}^{c}(j)=\frac{1}{2k}g_{ij}\mathbbm{1}_{\mathcal{N}_{k}}^{c}(i).\end{split}$
(21)

Use the fact that
$\mathbbm{1}_{\mathcal{N}_{k}^{c}}=1-\mathbbm{1}_{\mathcal{N}_{k}}$ and
$\partial\mathbbm{1}_{\mathcal{N}_{k}}=-\partial\mathbbm{1}_{\mathcal{N}_{k}^{c}}$
and vice versa:

$\begin{cases}\partial\mathbbm{1}_{\mathcal{N}_{k}}(i)&=\frac{g_{ij}}{2k}\left(\mathbbm{1}_{\mathcal{N}_{k}}(j)-\mathbbm{1}_{\mathcal{N}_{k}^{c}}(j)\right)=\frac{g_{ij}}{2k}\left(2\mathbbm{1}_{\mathcal{N}_{k}}(j)-1\right)\\\
\partial\mathbbm{1}_{\mathcal{N}_{k}}(j)&=\frac{g_{ij}}{2k}\left(\mathbbm{1}_{\mathcal{N}_{k}}(i)-\mathbbm{1}_{\mathcal{N}_{k}^{c}}(i)\right)=\frac{g_{ij}}{2k}\left(2\mathbbm{1}_{\mathcal{N}_{k}}(i)-1\right).\\\
\end{cases}$ (22)

Going backward through $\theta(\cdot)$ multiplies the gradient by $-\alpha$
(negative sign accounts for optimizing distance instead of similarity):

$\partial D(i)=-\frac{\alpha
g_{ij}}{2k}\left(2\mathbbm{1}_{\mathcal{N}_{k}}(j)-1\right)\>\text{and vice
versa for }\partial D(j).$ (23)

From the above equation, $\partial D(i,p)<0$ when
$\mathbbm{1}_{\mathcal{N}_{k}}(j,p)=1$ and $\partial D(i,p)>0$ when
$\mathbbm{1}_{\mathcal{N}_{k}}(j,p)=0$, for all samples $p$ in the batch. In
words, we increase the distance between $i$ and all samples in the
neighborhood of $j$, and we decrease the distance between $i$ and all points
outside the neighborhood of $j$. Recall that we assumed $i$ and $j$ to be a
negative pair, so it makes sense to pull $i$ away from the context of $j$ in
this fashion.

Table 1: R@$k$ Results. We use ResNet-50 with an embedding size of 512 for all
experiments. $\dagger$ indicates results reported by the original authors; we
re-run all other baselines using the implementation by Musgrave et al.
(2020b). Standard deviations are based on three trials with the same train-
test split. We emphasize that our R@1 performance is at least two standard
deviations better than the next best baseline on all datasets. We are at least
comparable to baselines on other R@$k$ metrics.

| CUB | Cars
---|---|---
Method | R@1 | R@2 | R@4 | R@8 | R@1 | R@2 | R@4 | R@8
Contrastive | 68.5 $\pm$ 0.3 | 78.3 $\pm$ 0.1 | 86.0 $\pm$ 0.2 | 91.3 $\pm$ 0.1 | 85.4 $\pm$ 0.2 | 91.1 $\pm$ 0.3 | 94.6 $\pm$ 0.3 | 96.8 $\pm$ 0.1
Triplet | 67.3 $\pm$ 0.2 | 77.9 $\pm$ 0.1 | 85.6 $\pm$ 0.2 | 91.2 $\pm$ 0.1 | 77.6 $\pm$ 1.3 | 85.4 $\pm$ 0.8 | 90.8 $\pm$ 0.7 | 94.1 $\pm$ 0.4
NtXent | 65.7 $\pm$ 0.4 | 76.3 $\pm$ 0.2 | 84.3 $\pm$ 0.4 | 90.0 $\pm$ 0.4 | 79.0 $\pm$ 0.6 | 86.0 $\pm$ 0.3 | 91.0 $\pm$ 0.2 | 94.4 $\pm$ 0.3
MS | 68.9 $\pm$ 0.5 | 78.5 $\pm$ 0.4 | 86.0 $\pm$ 0.6 | 91.4 $\pm$ 0.5 | 88.7 $\pm$ 0.4 | 93.0 $\pm$ 0.2 | 95.7 $\pm$ 0.1 | 97.3 $\pm$ 0.1
N-Softmax$\dagger$ | 61.3 | 73.9 | 83.5 | 90.0 | 84.2 | 90.4 | 94.4 | 96.9
ProxyNCA++$\dagger$ | 69.0 $\pm$ 0.8 |  79.8 $\pm$ 0.7 |  87.3 $\pm$ 0.7 |  92.7 $\pm$ 0.4 | 86.5 $\pm$ 0.4 | 92.5 $\pm$ 0.3 | 95.7 $\pm$ 0.2 |  97.7 $\pm$ 0.1
Fast-AP | 63.3 $\pm$ 0.1 | 73.7 $\pm$ 0.4 | 82.2 $\pm$ 0.3 | 88.5 $\pm$ 0.2 | 74.7 $\pm$ 0.4 | 82.5 $\pm$ 0.7 | 88.0 $\pm$ 0.6 | 92.2 $\pm$ 0.2
Smooth-AP | 66.5 $\pm$ 0.9 | 76.6 $\pm$ 0.5 | 84.8 $\pm$ 0.6 | 90.8 $\pm$ 0.4 | 81.1 $\pm$ 0.2 | 87.8 $\pm$ 0.4 | 92.2 $\pm$ 0.3 | 95.1 $\pm$ 0.3
ROADMAP | 68.7 $\pm$ 0.5 | 78.3 $\pm$ 0.3 | 86.1 $\pm$ 0.3 | 91.1 $\pm$ 0.1 | 84.5 $\pm$ 0.5 | 90.3 $\pm$ 0.0 | 93.9 $\pm$ 0.0 | 96.2 $\pm$ 0.1
Ours | 69.8 $\pm$ 0.2 |  79.8 $\pm$ 0.1 |  87.1 $\pm$ 0.1 | 92.3 $\pm$ 0.2 | 89.3 $\pm$ 0.0 |  93.7 $\pm$ 0.2 |  96.3 $\pm$ 0.1 |  97.8 $\pm$ 0.2

| SOP | mini-iNaturalist
---|---|---
Method | R@1 | R@10 | R@100 | R@1 | R@4 | R@16 | R@32
Contrastive | 82.4 $\pm$ 0.0 | 91.9 $\pm$ 0.0 | 96.0 $\pm$ 0.0 | 43.5 $\pm$ 0.1 | 62.7 $\pm$ 0.1 | 77.6 $\pm$ 0.1 | 83.2 $\pm$ 0.1
Triplet | 82.0 $\pm$ 0.0 | 92.5 $\pm$ 0.1 | 96.7 $\pm$ 0.0 | 35.4 $\pm$ 0.1 | 56.5 $\pm$ 0.1 | 74.7 $\pm$ 0.1 | 81.7 $\pm$ 0.1
NtXent | 79.7 $\pm$ 0.2 | 90.8 $\pm$ 0.0 | 96.1 $\pm$ 0.0 | 40.8 $\pm$ 0.1 | 61.6 $\pm$ 0.1 | 78.0 $\pm$ 0.0 | 83.9 $\pm$ 0.0
MS | 81.4 $\pm$ 0.0 | 91.4 $\pm$ 0.0 | 96.1 $\pm$ 0.1 | 44.9 $\pm$ 0.1 | 63.9 $\pm$ 0.1 | 78.4 $\pm$ 0.1 | 83.9 $\pm$ 0.1
N-Softmax$\dagger$ | 78.2 | 90.6 | 96.2 | – | – | – | –
ProxyNCA++$\dagger$ | 80.7 $\pm$ 0.5 | 92.0 $\pm$ 0.3 | 96.7 $\pm$ 0.1 | – | – | – | –
Fast-AP | 80.3 $\pm$ 0.1 | 91.0 $\pm$ 0.1 | 96.0 $\pm$ 0.0 | 35.6 $\pm$ 0.2 | 55.8 $\pm$ 0.1 | 72.8 $\pm$ 0.0 | 79.3 $\pm$ 0.0
Smooth-AP | 82.0 $\pm$ 0.0 | 92.6 $\pm$ 0.0 |  96.9 $\pm$ 0.0 | 42.7 $\pm$ 0.0 | 63.3 $\pm$ 0.0 | 79.0 $\pm$ 0.0 | 84.7 $\pm$ 0.0
ROADMAP | 83.1 $\pm$ 0.1 | 92.6 $\pm$ 0.0 | 96.6 $\pm$ 0.0 | 45.9 $\pm$ 0.1 |  65.8 $\pm$ 0.0 |  80.4 $\pm$ 0.1 |  85.7 $\pm$ 0.0
Ours | 83.3 $\pm$ 0.0 |  92.9 $\pm$ 0.1 | 96.7 $\pm$ 0.0 | 46.2 $\pm$ 0.0 |  65.8 $\pm$ 0.1 | 80.2 $\pm$ 0.1 | 85.4 $\pm$ 0.1

###  4.3 Step 3: Query Expansion

Query expansion (QE) is an established trick for image retrieval Arandjelović
& Zisserman (2012). QE expands the neighborhood set by additionally retrieving
the neighbors of very close neighbors. In the unsupervised setting, QE refers
to averaging the contextual similarity scores for samples in the
$\mathcal{R}_{k/2+\epsilon}$ neighborhood (see Eq. 13); this leads to more
robust pseudo-supervision. In our setting, we have two reasons for using QE
(Eq. 11 12 and 13): (1) averaging $\tilde{w}_{ij}$ with very close neighbors
could have the same effect as label smoothing and (2) the
$\mathcal{N}_{k/2+\epsilon}$ neighborhood is optimized by Eq. 12 and 13.
Empirically, QE is necessary to achieve the optimal performance of our
framework, since $\mathcal{L}_{\text{context}}$ achieves higher R@1 than
$\mathcal{L}_{2}$ in Table 3.

##  5 Experiments

Datasets  We experiment on two small-scale and two large-scale datasets:
Caltech-UCSD Birds (CUB-200) (Wah et al. (2011)), Stanford Cars-196 (Krause et
al. (2013)), Stanford Online Products (SOP) (Oh Song et al. (2016)), and mini-
iNaturalist-2021 (Van Horn et al. (2018)). CUB200 and Cars196 are smaller
fine-grain classification datasets with 200 and 196 unique labels,
respectively. SOP is a large-scale dataset with 120,053 product images from
22,634 classes. mini-iNaturalist-2021 is a subset of the iNaturalist 2021
species classification competition dataset, with 50 images each from 10,000
species. We are the first metric learning study to use this version of the
iNaturalist dataset, so we randomly withhold 5,000 species for testing and use
the other 5,000 for training. We will release our train-test split for
reproducibility.

Baselines  We choose a representative set of baselines for comparison in Table
1. Refer to Sections 1 and 2 for citations. Contrastive and triplet losses are
the accepted standard in the field. Multi-similarity (MS) is a popular pair-
wise similarity loss. From classification-inspired methods, we compare against
Normalized-Softmax (N-Softmax) (Zhai & Wu (2018)) and ProxyNCA++ (Teh et al.
(2020)). From AP maximization methods, we compare against Fast-AP, Smooth-AP,
and ROADMAP. For fair comparison, we re-run baselines using the same setup,
with the same learning rates and schedules. We use the published optimal
values for any internal hyperparameters. We were unable to match the published
performances for ProxyNCA++ and N-Softmax, so we use the published R@$k$
values for these two baselines.

Hyperparameters and Setup  We use hyperparameter values that are approximately
optimal across all datasets. $\lambda=0.4$, $\gamma=0.1$, $\alpha=10.0$,
$\epsilon=0.05$, $k=4$, $\delta_{+}=0.75$, $\delta_{-}=0.6$, $\tilde{s}=0.25$.
We follow the same training procedure as ROADMAP, but with a slightly faster
learning-rate schedule for time efficiency. We use Adam with a learning-rate
schedule that multiplies the learning rate by 0.3 at Epochs 15, 30, and 45. We
train for 80 epochs. We report results on the model with the best test R@1
metric, as is standard in the literature. We use an initial learning rate of
$8\times 10^{-5}$ on CUB, $0.00016$ on Cars, $4\times 10^{-5}$ on SOP, and
$8\times 10^{-5}$ on iNaturalist. We use a batch size of 256 for iNaturalist
and 128 for other datasets; the larger batch size is necessary to achieve
reasonable performance on iNaturalist. We use a 4 per class sampler. For CUB
and Cars, we use random sampling (sample 32 classes at random, then sample 4
images per class). For SOP and iNaturalist, we use hierarchical sampling
(Cakir et al. (2019)), following prior work. We use an embedding size of 512
and ResNet-50 with a linear embedding layer. Following prior work, we use
layer-norm and max-pooling on the smaller CUB and Cars datasets. We always
freeze batch-norm.

Performance Metrics  We report Recall @ $k$ for select $k$ in Table 1. R@$k$
is the percentage of test samples where at least one of the $k$ closest
neighbors have the same label. We report the average and standard deviations
across three trials with the same train-test split.

Discussion  Our R@1 results are at least two standard deviations better than
the best baseline across all datasets. Our results are especially good on CUB
and Cars, where we achieve R@1 gains of 0.8% and 0.6 %, resp. We achieve more
modest R@1 gains of 0.2% and 0.3% on SOP and iNaturalist, resp. We note that
these gains are significant because the standard deviation is relatively
small. We also note that while ProxyNCA++ and MS are the best baselines on CUB
and Cars, ROADMAP is the best baseline on SOP and iNaturalist. Our method is
the best across _all_ datasets, suggesting that it is more versatile.

Table 2: Ablation Results 1. We experiment with removing one or two of the three loss terms. Overall, all three losses are necessary to achieve the best R@1 results. Reasonable results are also achieved with just $\mathcal{L}_{\text{reg}}$ and $\mathcal{L}_{\text{contrast}}$. We emphasize that both $\mathcal{L}_{\text{context}}$ and $\mathcal{L}_{\text{reg}}$ are our contributions. $\mathcal{L}_{\text{context}}$ | $\mathcal{L}_{\text{reg}}$ | $\mathcal{L}_{\text{contrast}}$ | $\lambda$ | $\gamma$ | CUB R@1 | Cars R@1 | SOP R@1 | iNat. R@1
---|---|---|---|---|---|---|---|---
✓ |  | ✓ | 0.45 | 0 | 66.2 $\pm$ 0.5 | 87.9 $\pm$ 0.5 | 82.88 | 20.54
✓ | ✓ |  | 0.8 | 0.2 |  70.4 $\pm$ 0.6 | 82.8 $\pm$ 0.9 | 81.65 | 31.78
| ✓ | ✓ | 0 | 0.17 | 69.6 $\pm$ 0.4 | 88.5 $\pm$ 0.4 | 83.34 | 44.19
✓ |  |  | 1 | 0 |  70.3 $\pm$ 0.4 | 76.8 $\pm$ 1.0 | 81.91 | 21.63
✓ | ✓ | ✓ | 0.4 | 0.1 | 69.8 $\pm$ 0.3 |  89.2 $\pm$ 0.2 | 83.28 | 46.20
Table 3: Ablation Results 2. Here, we experiment with variations of $\mathcal{L}_{\text{context}}$, both by itself ($\lambda=1$) and regularized by a small amount of contrastive loss ($\lambda=0.8$). We exhaustively test various ways to simplify or modify the contextual loss presented in Eq. 9 \- 14. On the SOP dataset, we show that all of the modifications to $\mathcal{L}_{\text{context}}$ decrease R@1.  |  | SOP R@1 |
---|---|---|---
Ablation ($\gamma=0$) | $\lambda=0.8$ | $\lambda=1.0$ | Explanation
$\lambda\mathcal{L}_{1}$ | $+(1-\lambda)\mathcal{L}_{\text{contrast}}$ | 83.13 | 41.99 | Eq. 18 (step 1 only)
$\lambda\mathcal{L}_{1,\sigma}$ | $+(1-\lambda)\mathcal{L}_{\text{contrast}}$ | 79.20 | 75.61 | Eq. 18 but using $\theta_{\sigma}$ in Eq. 7
$\lambda\mathcal{L}_{2}$ | $+(1-\lambda)\mathcal{L}_{\text{contrast}}$ | 82.39 | 81.05 | Eq. 19 (skip step 3)
$\tilde{W}=\mathbbm{1}_{\mathcal{N}_{k+\epsilon}}$ (skip step 2) | 82.74 | 81.74 |
$\lambda\mathcal{L}_{\text{context},\min}$ | $+(1-\lambda)\mathcal{L}_{\text{contrast}}$ | 66.57 | – | Use $\min$ instead of $\odot$ for logical-and
$\lambda\mathcal{L}_{\text{context},\sigma}$ | $+(1-\lambda)\mathcal{L}_{\text{contrast}}$ | 82.39 | 77.84 | Use $\theta_{\sigma}$ for all steps
$\lambda\mathcal{L}_{\text{context},M_{+}}$ | $+(1-\lambda)\mathcal{L}_{\text{contrast}}$ | 82.31 | 80.72 | Remove $M_{-}$ term from Eq. 10
No stop gradient in Eq. 10 | 73.03 | – |
$\lambda\mathcal{L}_{\text{context}}$ | $+(1-\lambda)\mathcal{L}_{\text{contrast}}$ | 83.20 | 82.04 | Final results without $\mathcal{L}_{\text{reg}}$

###  5.1 Ablation Results

Table 2 shows R@1 on all datasets with one or two of the three losses removed.
This table shows that in general, the best R@1 performance is achieved by a
combination of the three losses in our framework. A combination of only
$\mathcal{L}_{\text{reg}}$ and $\mathcal{L}_{\text{contrast}}$ achieves better
performance than $\mathcal{L}_{\text{contrast}}$ alone across all datasets
(see first row of Table 1). The last row of Table 2 shows that adding
$\mathcal{L}_{\text{context}}$ achieves 0.7% and 2% R@1 gains on Cars and
iNaturalist, resp. These observations show that both of our contributions
($\mathcal{L}_{\text{reg}}$ and $\mathcal{L}_{\text{context}}$) are necessary
to achieve sate-of-the-art performance.

In Table 3, we evaluate the contribution of each step in the calculation of
$\mathcal{L}_{\text{context}}$. These experiments focus on dissecting
$\mathcal{L}_{\text{context}}$, so we set $\gamma=0$ (no similarity
regularizer) and $\lambda=0.8$ or $1.0$. Overall, all of the modifications
tested in Table 3 decrease the performance of $\mathcal{L}_{\text{context}}$
in terms of R@1. In particular: row 1 shows that including only step 1 leads
to a collapsed representation when $\lambda=1.0$; rows 3 and 4 show that steps
3 and 2 of the loss calculation are necessary; row 7 shows that the $M_{-}$
term in Eq. 10 is necessary; rows 2 and 6 show that our approach to optimizing
$\theta$ in Eq. 6 is better than using a sigmoid. We defer hyperparameter
tuning results to the Appendix.

##  6 Conclusion

In this work we propose a novel supervised metric learning loss based on
contextual similarity optimization. Our contextual loss improves the R@1
performance significantly over the current state-of-the-art on four diverse
benchmarks, when regularized by the contrastive loss and a similarity
regularizer. As our core technical contribution, we propose a highly non-
trivial framework to optimize the contextual similarity, overcoming several
non-differentiability issues. Namely, we show that the heaviside function can
be effectively optimized by setting a positive gradient in the backward pass,
and the intersection between sets can be optimized by indicator matrix
multiplication. We carefully justify our framework by analyzing the gradient
of the loss and with an exhaustive ablation study.

#### Ethics Statement

We note that metric learning can be applied to controversial problems such as
person re-identification and face re-identification. Our work is mainly
foundational, so does not contribute directly to these applications. We also
limit our experimentation to the image retrieval aspect of metric learning.

#### Reproducibility Statement

We will release code publicly on GitHub. Our main results are reproducible by
running the code. Instructions on how to run the code will be provided in a
README file. We include details on hardware requirements in the Appendix. Code
is released here: https://github.com/Chris210634/metric-learning-using-
contextual-similarity

#### Acknowledgments

DISTRIBUTION STATEMENT A. Approved for public release. Distribution is
unlimited.

This material is based upon work supported by the Department of the Air Force
under Air Force Contract No. FA8702-15-D-0001. Any opinions, findings,
conclusions or recommendations expressed in this material are those of the
author(s) and do not necessarily reflect the views of the Department of the
Air Force.

## References

  * Arandjelović & Zisserman (2012) Relja Arandjelović and Andrew Zisserman.  Three things everyone should know to improve object retrieval.  In _2012 IEEE conference on computer vision and pattern recognition_ , pp. 2911–2918. IEEE, 2012.
  * Boudiaf et al. (2020) Malik Boudiaf, Jérôme Rony, Imtiaz Masud Ziko, Eric Granger, Marco Pedersoli, Pablo Piantanida, and Ismail Ben Ayed.  A unifying mutual information view of metric learning: cross-entropy vs. pairwise losses.  In _European conference on computer vision_ , pp. 548–564. Springer, 2020.
  * Brown et al. (2020) Andrew Brown, Weidi Xie, Vicky Kalogeiton, and Andrew Zisserman.  Smooth-ap: Smoothing the path towards large-scale image retrieval.  In _European Conference on Computer Vision_ , pp. 677–694. Springer, 2020.
  * Cakir et al. (2019) Fatih Cakir, Kun He, Xide Xia, Brian Kulis, and Stan Sclaroff.  Deep metric learning to rank.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pp. 1861–1870, 2019.
  * Cho et al. (2014) Kyunghyun Cho, Bart Van Merriënboer, Dzmitry Bahdanau, and Yoshua Bengio.  On the properties of neural machine translation: Encoder-decoder approaches.  _arXiv preprint arXiv:1409.1259_ , 2014.
  * Ge et al. (2020a) Yixiao Ge, Dapeng Chen, and Hongsheng Li.  Mutual mean-teaching: Pseudo label refinery for unsupervised domain adaptation on person re-identification.  _arXiv preprint arXiv:2001.01526_ , 2020a.
  * Ge et al. (2020b) Yixiao Ge, Feng Zhu, Dapeng Chen, Rui Zhao, and Hongsheng Li.  Self-paced contrastive learning with hybrid memory for domain adaptive object re-id.  In _Advances in Neural Information Processing Systems_ , 2020b.
  * Hadsell et al. (2006) Raia Hadsell, Sumit Chopra, and Yann LeCun.  Dimensionality reduction by learning an invariant mapping.  In _2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’06)_ , volume 2, pp. 1735–1742. IEEE, 2006.
  * Kim et al. (2022) Sungyeon Kim, Dongwon Kim, Minsu Cho, and Suha Kwak.  Self-taught metric learning without labels.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 7431–7441, 2022.
  * Krause et al. (2013) Jonathan Krause, Michael Stark, Jia Deng, and Li Fei-Fei.  3d object representations for fine-grained categorization.  In _Proceedings of the IEEE international conference on computer vision workshops_ , pp. 554–561, 2013.
  * Liu et al. (2016) Ziwei Liu, Ping Luo, Shi Qiu, Xiaogang Wang, and Xiaoou Tang.  Deepfashion: Powering robust clothes recognition and retrieval with rich annotations.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 1096–1104, 2016.
  * McFee & Lanckriet (2010) Brian McFee and Gert RG Lanckriet.  Metric learning to rank.  In _ICML_ , 2010.
  * Movshovitz-Attias et al. (2017) Yair Movshovitz-Attias, Alexander Toshev, Thomas K Leung, Sergey Ioffe, and Saurabh Singh.  No fuss distance metric learning using proxies.  In _Proceedings of the IEEE International Conference on Computer Vision_ , pp. 360–368, 2017.
  * Musgrave et al. (2020a) Kevin Musgrave, Serge Belongie, and Ser-Nam Lim.  A metric learning reality check.  In _European Conference on Computer Vision_ , pp. 681–699. Springer, 2020a.
  * Musgrave et al. (2020b) Kevin Musgrave, Serge Belongie, and Ser-Nam Lim.  Pytorch metric learning, 2020b.
  * Oh Song et al. (2016) Hyun Oh Song, Yu Xiang, Stefanie Jegelka, and Silvio Savarese.  Deep metric learning via lifted structured feature embedding.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 4004–4012, 2016.
  * Ramzi et al. (2021) Elias Ramzi, Nicolas Thome, Clément Rambour, Nicolas Audebert, and Xavier Bitot.  Robust and decomposable average precision for image retrieval.  _Advances in Neural Information Processing Systems_ , 34:23569–23581, 2021.
  * Teh et al. (2020) Eu Wern Teh, Terrance DeVries, and Graham W Taylor.  Proxynca++: Revisiting and revitalizing proxy neighborhood component analysis.  In _European Conference on Computer Vision_ , pp. 448–464. Springer, 2020.
  * Van Horn et al. (2018) Grant Van Horn, Oisin Mac Aodha, Yang Song, Yin Cui, Chen Sun, Alex Shepard, Hartwig Adam, Pietro Perona, and Serge Belongie.  The inaturalist species classification and detection dataset.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 8769–8778, 2018.
  * Wah et al. (2011) Catherine Wah, Steve Branson, Peter Welinder, Pietro Perona, and Serge Belongie.  The caltech-ucsd birds-200-2011 dataset.  2011\.
  * Wang et al. (2019) Xun Wang, Xintong Han, Weilin Huang, Dengke Dong, and Matthew R Scott.  Multi-similarity loss with general pair weighting for deep metric learning.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 5022–5030, 2019.
  * Weinberger et al. (2005) Kilian Q Weinberger, John Blitzer, and Lawrence Saul.  Distance metric learning for large margin nearest neighbor classification.  _Advances in neural information processing systems_ , 18, 2005.
  * Wu et al. (2017) Chao-Yuan Wu, R Manmatha, Alexander J Smola, and Philipp Krahenbuhl.  Sampling matters in deep embedding learning.  In _Proceedings of the IEEE international conference on computer vision_ , pp. 2840–2848, 2017.
  * Ye et al. (2021) Mang Ye, Jianbing Shen, Gaojie Lin, Tao Xiang, Ling Shao, and Steven CH Hoi.  Deep learning for person re-identification: A survey and outlook.  _IEEE transactions on pattern analysis and machine intelligence_ , 44(6):2872–2893, 2021.
  * Zhai & Wu (2018) Andrew Zhai and Hao-Yu Wu.  Classification is a strong baseline for deep metric learning.  _arXiv preprint arXiv:1811.12649_ , 2018.
  * Zhong et al. (2017) Zhun Zhong, Liang Zheng, Donglin Cao, and Shaozi Li.  Re-ranking person re-identification with k-reciprocal encoding.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 1318–1327, 2017.

##  Appendix A Appendix

###  A.1 Code and Environment

We include code with our submission. We run experiments on 1 V100 GPU with 16
GB of memory. The CUB and Cars experiments take under one hour. The SOP and
iNaturalist experiments take 4 hours and 6 hours, respectively. Some of our
code is borrowed from ROADMAP. For faster experimentation, we use mixed
precision floating point. Experiments take more than 12 hours on P100 GPUs,
partially because mixed precision arithmetic does not appear to speed-up
experiments as much on P100 GPUs compared to on V100 GPUs.

###  A.2 Minor Experimental Details

Augmentation  On CUB, Cars and SOP, we use random resized crop with default
parameters to crop the image to 256 $\times$ 256 pixels. We horizontally flip
the image with 50% probability. For testing, we resize the image to 288
pixels, then center crop to 256 pixels. On iNaturalist, we random resize crop
to 224 $\times$ 224 pixels, then flip horizontally with 50% probability. We
use a smaller image size on iNaturalist so that the larger batch size can fit
in the 16 GB of GPU memory. For testing, we resize the image to 256 pixels
then center crop to 224 pixels.

Sampling  We use a batch size of 128 on CUB, Cars and SOP. We use a batch size
of 256 on iNaturalist. On all datasets, we train for 80 epochs. An epoch is
defined as 15, 30, 655, and 576 batches for CUB, Cars, SOP, and iNaturalist
respectively. On CUB and Cars, we use a random 4 per class sampler. On SOP and
iNaturalist, we use a balanced hierarchical sampling strategy, since there are
super-labels. In particular, half of each batch is randomly sampled from one
super-label, and the other half is randomly sampled from another super-label.
In order to sample in a balanced manner (to prevent over-emphasis of uncommon
super-labels), we arrange all samples in the training dataset into half-
batches with same super-labels at the beginning of each epoch. We then form
batches by pairing up half-batches in a round-robin fashion.

Hyperparameters  We use the same set of hyperparameters for all experiments.
We present SOP R@1 results for various hyperparameter values in Figure 6. In
Figure 5, we show the R@1 result on CUB and Cars for varying $\lambda$ when
$\gamma=0$. We compare this result against ROADMAP with varying $\lambda$.
ROADMAP is a recent state-of-the-art metric learning method, which uses a
$\lambda$ hyperparameter to balance between an AP loss and the contrastive
loss. We show that even if $\lambda$ and learning rates are tuned for ROADMAP
on CUB and Cars, our method still achieves higher R@1 when the optimal
$\lambda$ and learning rates are chosen.

Loss Plot  Figure 7 plots the contextual loss and test R@1 over the course of
training on SOP, with varying $\lambda$. $\gamma=0$. These plots show that the
loss decreases over the course of training; in general, a lower contextual
loss corresponds to a higher test R@1. This result suggests that the
contextual loss is a reasonable objective for image retrieval.

Note on Contrastive Loss  The contrastive loss has two margin parameters
$\delta_{-}$ and $\delta_{+}$. The optimal values for $\delta_{+}$, the
positive margin, is different depending on the value of $\gamma$ (how much
similarity regularization is used). For the “contrastive” results in row 1 of
Table 1, we use $\delta_{-}=0.6$ and $\delta_{+}=0.9$. For our results where
$\gamma=0.1$, we set $\delta_{-}=0.6$ and $\delta_{+}=0.75$. We find this
tighter positive margin to be optimal under similarity regularization.
However, we note that $\delta_{+}=0.75$ is likely to be too small when only
the contrastive loss is used, since positive pairs are not encouraged to be
more similar than a cosine similarity of 0.75.

Note on Stop Gradient  Note that we detach the normalization factors in Eq.
10. This is necessary because optimizing the size of the neighborhood is
undesirable. On the contrary, we do not detach the normalization factor in Eq.
13. This is intentional even if somewhat un-intuitive. We show in Figure 9
that detaching $|\mathcal{R}_{k/2+\epsilon}(i)|$ in Eq. 13 increases R@1 on
SOP when $\lambda$ is high, but decreases R@1 in all other scenarios.

In-Shop Results  We offer R@$k$ results for standard $k$ values on the In-Shop
dataset in Table 5. In-Shop (Liu et al. (2016)) is a popular image retrieval
benchmark. We use the standard train-test split: we use 25,882 images from
3,997 classes for training; for testing, we use a query set with 14,218 images
and a gallery set (aka. reference set) with 12,612 images, both from 3,985
classes. We use the same augmentation and batch size as SOP. An epoch is
defined as 600 batches.

Average Precision Results  In the main paper, we focused on the R@$k$ metric,
as is standard in the image retrieval literature. We offer Average Precision
(AP) results in Table 4. We present results for two flavors of AP (McFee &
Lanckriet (2010)) : mAP and mAP@R. mAP is the precision at $k$, averaged
across correctly retrieved samples, averaged across the query set. mAP@R
(Musgrave et al. (2020a)) only averages across retrievals up to the number of
positive pairs in the gallery set. Hence, mAP@R is always smaller than mAP.
The two versions of AP are equivalent bases for comparison. In the following
equations, $\mathcal{X}^{+}_{i}$ denotes the set of samples in the gallery set
with the same label as query $i$ and $\mathcal{X}^{-}_{i}$ denotes the set of
samples in the gallery set with a different label than query $i$. Prec@$k$
denotes the precision at $k$, which is the percentage of samples ranked less
than or equal to $k$ with the same label as the query.

$\text{mAP}_{i}=\frac{1}{|\mathcal{X}^{+}_{i}|}\sum_{k=1}^{|\mathcal{X}^{+}_{i}|+|\mathcal{X}^{-}_{i}|}\text{Prec@$k$}\>\mathbbm{1}[k\in\mathcal{X}^{+}_{i}]$
$\text{mAP@R}_{i}=\frac{1}{|\mathcal{X}^{+}_{i}|}\sum_{k=1}^{|\mathcal{X}^{+}_{i}|}\text{Prec@$k$}\>\mathbbm{1}[k\in\mathcal{X}^{+}_{i}]$

Experimental Results on Similarity Regularizer  In the main paper, we
presented the similarity regularizer $\mathcal{L}_{\text{reg}}$ as
complementary to the contextual similarity loss
$\mathcal{L}_{\text{context}}$. However, there is no theoretical reason for
limiting the regularizer to our contextual similarity framework. We show in
Table 6 that the similarity regularizer offers an improvement in retrieval
performance when combined with some baselines. In particular, the regularizer
consistently improves R@1 performance of Triplet, Multi-Similarity, Smooth-AP,
and ROADMAP losses across Cars, SOP and iNaturalist benchmarks. This is a
promising result.

Figure 5: Comparison of our method to ROADMAP with varying $\lambda$.
$\gamma=0$. The orange lines show R@1 of our method with varying $\lambda$ and
learning rates. The purple lines show R@1 of ROADMAP with varying $\lambda$
and learning rates. $\lambda=0$ corresponds to the contrastive loss. Optimal
performance of our method is always better than ROADMAP on CUB and Cars.

Figure 6: Hyperparameter tuning on SOP. We experiment with various values for
$\epsilon$, $\gamma$ and $\tilde{s}$. The R@1 values are similar for the
different hyperparameter choices.

Figure 7: Plot of contextual similarity loss $\mathcal{L}_{\text{context}}$
with varying $\lambda$. Each color is a different $\lambda$. Results are on
the SOP dataset. Observe that the contrastive loss implicitly minimizes the
contextual loss (blue line). Also observe that in general, a lower contextual
loss corresponds to a higher test R@1. This shows that our loss function is a
reasonable objective. Note that the plot on the left plots
$\mathcal{L}_{\text{context}}$, not the complete training loss. $\gamma=0$.

Figure 8: Ablation results on CUB and Cars. The blue line represents test R@1
of our framework for varying $\lambda$ with $\gamma=0$. The other colors
represent various modifications or simplifications of the loss function
(reference Table 3 and Section 4 for full description). We perform each
ablation experiment with varying $\lambda$ since the optimal $\lambda$ is not
constant across all experiments. Clearly, the modified versions of
$\mathcal{L}_{\text{context}}$ are sub-optimal.

Figure 9: Investigating the effect of detaching $|\mathcal{R}_{k/2+\epsilon}(i)|$ in Eq. 13. The blue line shows R@1 of our framework with varying $\lambda$. $\gamma=0$. The orange line shows R@1 after detaching $|\mathcal{R}_{k/2+\epsilon}(i)|$. Overall, detaching this normalization factor is undesirable. Table 4: mAP and mAP@R Results. We use ResNet-50 with an embedding size of 512 for all experiments. We re-run all baselines in this Table using the implementation by Musgrave et al. (2020b). Standard deviations are based on three trials with the same train-test split.  | CUB | Cars | SOP | mini-iNaturalist
---|---|---|---|---
Method | mAP | mAP@R | mAP | mAP@R | mAP | mAP@R | mAP | mAP@R
Contrastive | 36.6 $\pm$ 0.4 | 26.6 $\pm$ 0.5 | 38.2 $\pm$ 0.4 | 28.8 $\pm$ 0.5 | 62.9 $\pm$ 0.1 | 57.0 $\pm$ 0.1 | 16.0 $\pm$ 0.1 | 11.6 $\pm$ 0.0
Triplet | 36.9 $\pm$ 0.4 | 26.5 $\pm$ 0.2 | 34.9 $\pm$ 0.8 | 24.7 $\pm$ 0.7 | 63.1 $\pm$ 0.0 | 56.4 $\pm$ 0.0 | 12.1 $\pm$ 0.0 | 7.9 $\pm$ 0.0
NtXent | 36.3 $\pm$ 0.3 | 26.0 $\pm$ 0.3 | 36.2 $\pm$ 0.7 | 26.3 $\pm$ 0.7 | 60.1 $\pm$ 0.1 | 53.4 $\pm$ 0.1 | 15.9 $\pm$ 0.0 | 10.6 $\pm$ 0.0
MS | 38.0 $\pm$ 0.1 | 27.7 $\pm$ 0.1 |  43.1 $\pm$ 0.3 |  33.7 $\pm$ 0.3 | 62.5 $\pm$ 0.0 | 56.2 $\pm$ 0.1 | 16.6 $\pm$ 0.0 | 11.9 $\pm$ 0.0
Fast-AP | 34.4 $\pm$ 0.8 | 24.4 $\pm$ 0.7 | 32.7 $\pm$ 0.0 | 23.5 $\pm$ 0.0 | 60.7 $\pm$ 0.1 | 54.2 $\pm$ 0.1 | 13.9 $\pm$ 0.1 | 9.1 $\pm$ 0.1
Smooth-AP | 36.8 $\pm$ 0.4 | 26.4 $\pm$ 0.4 | 37.5 $\pm$ 0.3 | 27.5 $\pm$ 0.2 | 63.1 $\pm$ 0.2 | 56.4 $\pm$ 0.2 | 16.5 $\pm$ 0.0 | 11.3 $\pm$ 0.0
ROADMAP | 37.7 $\pm$ 0.1 | 27.5 $\pm$ 0.2 | 38.3 $\pm$ 0.4 | 28.7 $\pm$ 0.4 | 64.4 $\pm$ 0.1 | 58.2 $\pm$ 0.1 |  17.8 $\pm$ 0.0 |  12.7 $\pm$ 0.0
Ours |  38.3 $\pm$ 0.2 |  28.0 $\pm$ 0.2 | 42.4 $\pm$ 0.4 | 33.0 $\pm$ 0.3 |  65.0 $\pm$ 0.1 |  58.6 $\pm$ 0.1 | 17.1 $\pm$ 0.0 | 12.4 $\pm$ 0.0
Table 5: In-Shop Results. We use ResNet-50 with an embedding size of 512 for all experiments. $\dagger$ indicates results reported by the original authors; we re-run all other baselines using the implementation by Musgrave et al. (2020b). Our R@1 performance is better than all baselines. We are at least comparable to baselines on other R@$k$ metrics. | In-Shop
---|---
Method | R@1 | R@10 | R@20 | R@30 | R@40 | R@50
Contrastive | 90.1 | 97.4 | 98.3 | 98.6 | 98.8 | 98.9
Triplet | 90.2 | 98.0 | 98.7 | 99.0 | 99.2 | 99.3
NtXent | 89.3 | 97.6 | 98.3 | 98.7 | 98.9 | 99.0
MS | 86.9 | 96.1 | 97.3 | 97.9 | 98.1 | 98.3
N-Softmax$\dagger$ | 88.6 | 97.5 | 98.4 | 98.8 | – | –
ProxyNCA++$\dagger$ | 90.4 | 98.1 | 98.8 | 99.0 | 99.2 | –
Fast-AP | 89.7 | 97.5 | 98.3 | 98.6 | 98.8 | 98.9
Smooth-AP | 90.2 | 97.9 | 98.7 | 99.0 | 99.2 | 99.2
ROADMAP | 90.4 | 97.6 | 98.3 | 98.6 | 98.8 | 99.0
Ours | 90.7 | 97.8 | 98.5 | 98.9 | 99.1 | 99.2
Table 6: Similarity Regularizer Results. Here, we experiment with naively
adding our similarity regularizer to a subset of the baseline methods. The
overall loss is $\gamma\mathcal{L}_{\text{reg}}+(1-\gamma)(\text{baseline
loss})$ where $\gamma=0.1$. We try two different values for $\tilde{s}$ and
compare the difference in various performance metrics with the original
baseline. The colored subscript denotes the increase or decrease in each
metric when our similarity regularizer is added. Observe that in most cases,
the regularizer improves the performance of the baseline or has negligible
effect, with an appropriate $\tilde{s}$ value.

| Cars
---|---
| mAP | mAP@R | R@1
Method | No $\mathcal{L}_{\text{reg}}$ | $\tilde{s}=0.25$ | $\tilde{s}=0.3$ | No $\mathcal{L}_{\text{reg}}$ | $\tilde{s}=0.25$ | $\tilde{s}=0.3$ | No $\mathcal{L}_{\text{reg}}$ | $\tilde{s}=0.25$ | $\tilde{s}=0.3$
Contrastive | 38.2 | 40.6 +2.4 | 39.4 +1.2 | 28.8 | 31.2 +2.4 | 30.0 +1.2 | 85.4 | 88.0 +2.6 | 87.1 +1.7
Triplet | 34.9 | 39.2 +4.3 | 38.5 +3.6 | 24.7 | 29.9 +5.2 | 29.2 +4.5 | 77.6 | 87.8 +10.2 | 86.7 +9.1
NtXent | 36.2 | 37.0 +0.8 | 35.8 -0.4 | 26.3 | 27.0 +0.7 | 26.0 -0.3 | 79.0 | 79.9 +0.9 | 78.7 -0.3
MS | 43.1 | 42.7 -0.4 |  43.6 +0.5 | 33.7 | 33.3 -0.4 | 34.1 +0.4 | 88.7 | 88.8 +0.1 | 89.3 +0.6
Smooth-AP | 37.5 | 38.9 +1.4 | 38.1 +0.6 | 27.5 | 28.8 +1.3 | 28.2 +0.7 | 81.1 | 82.1 +1.0 | 82.4 +1.3
ROADMAP | 38.3 | 39.2 +0.9 | 39.6 +1.3 | 28.7 | 29.6 +0.9 | 30.0 +1.3 | 84.5 | 85.7 +1.2 | 85.7 +1.2
Ours | – | 42.4 | – | – | 33.0 | – | – | 89.3 | –

| SOP
---|---
| mAP | mAP@R | R@1
Method | No $\mathcal{L}_{\text{reg}}$ | $\tilde{s}=0.25$ | $\tilde{s}=0.3$ | No $\mathcal{L}_{\text{reg}}$ | $\tilde{s}=0.25$ | $\tilde{s}=0.3$ | No $\mathcal{L}_{\text{reg}}$ | $\tilde{s}=0.25$ | $\tilde{s}=0.3$
Contrastive | 62.9 | 62.8 -0.1 | 62.7 -0.2 | 57.0 | 56.9 -0.1 | 56.8 -0.2 | 82.4 | 82.3 -0.1 | 82.3 -0.1
Triplet | 63.1 | 64.4 +1.3 | 64.6 +1.5 | 56.4 | 57.8 +1.4 | 58.0 +1.6 | 82.0 | 83.1 +1.1 | 83.3 +1.3
NtXent | 60.1 | 60.3 +0.2 | 60.0 -0.1 | 53.4 | 53.7 +0.3 | 53.4 +0.0 | 79.7 | 80.0 +0.3 | 79.7 +0.0
MS | 62.5 | 62.4 -0.1 | 62.6 +0.1 | 56.2 | 56.1 -0.1 | 56.3 +0.1 | 81.4 | 81.4 +0.0 | 81.8 +0.4
Smooth-AP | 63.1 | 63.2 +0.1 | 63.5 +0.4 | 56.4 | 56.5 +0.1 | 56.8 +0.4 | 82.0 | 82.1 +0.1 | 82.3 +0.3
ROADMAP | 64.4 | 64.4 +0.0 | 64.4 +0.0 | 58.2 | 58.2 +0.0 | 58.2 +0.0 | 83.1 | 83.1 +0.0 | 83.1 +0.0
Ours | – | 65.0 | – | – | 58.6 | – | – | 83.3 | –

| iNaturalist
---|---
| mAP | mAP@R | R@1
Method | No $\mathcal{L}_{\text{reg}}$ | $\tilde{s}=0.25$ | $\tilde{s}=0.3$ | No $\mathcal{L}_{\text{reg}}$ | $\tilde{s}=0.25$ | $\tilde{s}=0.3$ | No $\mathcal{L}_{\text{reg}}$ | $\tilde{s}=0.25$ | $\tilde{s}=0.3$
Contrastive | 16.0 | 15.1 -0.9 | 15.5 -0.5 | 11.6 | 11.2 -0.4 | 11.5 -0.1 | 43.5 | 43.4 -0.1 | 43.7 +0.2
Triplet | 12.1 | 13.3 +1.2 | 13.6 +1.5 | 7.9 | 9.8 +1.9 | 10.0 +2.1 | 35.4 | 41.4 +6.0 | 41.7 +6.3
NtXent | 15.9 | 15.9 +0.0 | 15.9 +0.0 | 10.6 | 10.7 +0.1 | 10.6 +0.0 | 40.8 | 40.7 -0.1 | 40.7 -0.1
MS | 16.6 | 16.6 +0.0 | 16.6 +0.0 | 11.9 | 11.9 +0.0 | 11.9 +0.0 | 44.9 | 44.9 +0.0 | 45.0 +0.1
Smooth-AP | 16.5 | 16.6 +0.1 | 16.7 +0.2 | 11.3 | 11.4 +0.1 | 11.4 +0.1 | 42.7 | 43.1 +0.4 | 43.2 +0.5
ROADMAP | 17.8 |  18.1 +0.3 | 18.0 +0.2 | 12.7 |  13.1 +0.4 | 13.0 +0.3 | 45.9 | 47.1 +1.2 | 46.8 +0.9
Ours | – | 17.1 | – | – | 12.4 | – | – | 46.2 | –


††thanks: Also at Acoustic Science and Technology Laboratory, Harbin
Engineering University, Harbin 150001,China.††thanks:

# Subspace-based compressive sensing algorithm for raypath separation in a
shallow-water waveguide

Longyu Jiang  The Laboratory of Image Science and Technology, Southeast
University, Nanjing 210096, China.  Zhe Zhang  The Laboratory of Image Science
and Technology, Southeast University, Nanjing 210096, China.  Rui Jin  The
Laboratory of Image Science and Technology, Southeast University, Nanjing
210096, China.  Xiao Zhou  The Laboratory of Image Science and Technology,
Southeast University, Nanjing 210096, China<EMAIL_ADDRESS>Philippe Roux
Institut des Sciences de la Terre, Université Joseph Fourier, Centre National
de la Recherche Scientifique, 1381 Rue de la Piscine, Saint-Martin d’Hères,
France.

###### Abstract

Compressive sensing (CS) has been applied to estimate the direction of arrival
(DOA) in underwater acoustics. However, the key problem needed to be resolved
in a multipath propagation environment is to suppress the interferences
between the raypaths. Thus, in this paper, a subspace-based compressive
sensing algorithm that formulates the statistic information of the signal
subspace in a CS framework is proposed. The experiment results show that (1)
the proposed algorithm enables the separation of raypaths that arrive closely
at the receiver array and (2) the existing algorithms fail, especially in a
low signal-to-noise ratio (SNR) environment.

##  I Introduction

Multipath acoustic waves propagate in shallow water due to the reflection and
(or) refraction of acoustic eigenrays. In ocean acoustic
tomography(OAT)Munk1995 , we use the variation in the arrival time of rays to
inversely derive fluctuations of sound speed profile in the water column.
Recent works have emphasized the use of angle variations to image sound speed
perturbations using the sensitivity kernel approach AulanierShallow2013 ;
AulanierTime2013 . Multiple paths cover the different parts of the ocean, so
they can provide spatially-extended information for the inversion process.
However, they also produce interferences between them. Thus, successful
separation of these raypaths based on direction-of-arrival (DOA) estimation is
the first step of OAT. Researchers have paid great attention to solving this
problem, and many algorithms have been proposedschimidt1986 ; jiang2017raypath
; jiang2017Active . Conventional beamforming (CBF) and Capon BF are two such
conventional algorithms, but they fail in separating raypaths that arrive
close to each other at the receiver array. Moreover, Iturbe and co-authors
proposed a double-beamforming algorithm in a double-array configuration to
improve the performance of the conventional beamforming algorithm Iturbe2008 .
On the other hand, high-resolution algorithms are also introduced to improve
the angular resolution of beamforming. The MUltiple SIgnal Classification
(MUSIC) algorithm schimidt1986  is a representative of the high-resolution
algorithms, which is based on the orthogonality of the signal and noise
subspace. Following the strategy of high-resolution algorithms, exploring the
spectral information of the emitted signal, Jiang et al proposed a smoothing-
MUSICAL algorithm jiang2017raypath  to separate the correlated or coherent
raypaths and a higher-order algorithmjiang2017Active  in the frequency domain
to enable the separation the raypaths in an ocean environment in the presence
of colored noise. Experimental results illustrate that both algorithms achieve
resolution improvement compared to beamforming-like algorithms. However, the
noise and signal subspaces are not always totally orthogonal in actual
operationthomas1995 . Thus, the DOA estimation based on subspace decomposition
still faces great challenges,especially when multiple raypaths propagate in a
strongly noisy ocean environment. Recently, the DOA estimation based on CS has
been very attractive to many researchers in this area. CS theory as a new
information technology makes use of the fact that received signals with
inherent sparsity can be reconstructed perfectly from very few measurements
through the solution of a convex minimization problem. In underwater
acoustics, Edelmann and Gaumond Edelmann2011Beamforming  first applied a CS
algorithm to estimate DOAs using ship-towed horizontal arrays. Compared with
conventional beamforming, the algorithm achieved a resolution improvement and
interference suppression. Xenaki and Gerstoft Xenaki2014Compressive  developed
a compressive beamforming algorithm and further discussed its performance,
especially with coherent arrivals, single-snapshot data and random array
geometries. In addition, the authors indicated that the limitations of CS-
based algorithms are related to the beampattern. Gerstoft et al.
gerstoft2015multiple  extended the algorithm to the case of multiple snapshots
combining a maximum a posteriori method with the least absolute shrinkage and
selection operator (LASSO) method. The experiment shows that this method also
outperforms conventional high-resolution methods even under challenging
scenarios such as coherent arrivals and single-snapshot data. Recently, a
compressive sensing beamforming method (CSB-II) zhong2013compressive  is
developed based on sampling covariance matrix in the application context of
aeroacoustic experiments containing strong background noise at broadband
frequencies. The experiment results illustrate that the CSB-II method achieves
a higher resolution and is robust to sensing noise for presumably spatially
sparse and incoherent signals. However, in the context of ocean acoustic
tomography, the multiple raypaths are produced by the emitted signal through
the refraction or the reflection of the ocean surface and the continental
shelf. Thus, they are always fully correlated or coherent. The scattering at
different points on the sea surface or the eigenpaths crossing different areas
of the water column may attenuate these correlations; however, there is still
the possibility that the case of fully correlated or coherent raypaths exists.
This is a very difficult separation case. In addition, even in the case of
partly correlated raypaths, the correlations between the raypaths impede the
successful separation or decrease the separation resolution of raypaths.
However, the existing CS-based algorithms do not focus on solving this
problem.

In order to perform the separation between eigenrays, the present paper
proposes a subspace-based compressive sensing algorithm that formulates the
statistic information of the signal subspace in a CS framework. The signal
subspace is first obtained through an eigenvalues decomposition to the
spectral matrix of the received signal. Then, to suppress the interferences
produced by the correlations between the raypaths, the sparse signal
reconstruction is achieved by designing a subspace-based $L_{1}$
regularization optimization.

The paper is organized as follows. In Section II, the signal model is
established. In Section III, the principle of the comparative method–reweighed
compressive sensing is described. In Section IV, we introduce our proposed
algorithm. In Section V, we illustrate the performance of the proposed
algorithm based on simulations, a tank experiment and an ocean experiment.

##  II Signal model

In a shallow-water waveguide, we assume that $P$ raypaths emitted by a far-
field source $e(t)$ are received by each element of a vertical array composed
of $M$ sensors. $s_{p}(t)=a_{p}e(t)$ denotes the $p^{th}$ raypath, where
$a_{p}$ is the random amplitude of the $p^{th}$ path. The received signal on
the $m^{th}~{}(m=1,2,...,M)$ sensor, denoted by $y_{m}(t)$, can be modeled in
the time domain as follows.

${y_{m}}(t)=\sum\limits_{p=1}^{P}{{s_{p}}(t-{\tau_{mp}})+{n_{m}}(t)}$ (1)

where $s_{p}$ denotes the amplitude of the $p^{th}$ source on the $m^{th}$
sensor; and $n_{m}(t)$ is the white Gaussian noise received at the $m^{th}$
sensor. $\tau_{mp}$ represents the propagation delay of the $p^{th}$ raypath
between the $m^{th}$ sensor and the reference sensor.

In the frequency domain, Eq. 1 can be rewritten as

${Y_{m}}(\nu)=\sum\limits_{p=1}^{P}{{S_{p}}(\nu)G_{mp}+{N_{m}}(\nu)}$ (2)

with $G_{mp}=exp(-j\nu(m-1)2\pi\frac{dsin\theta_{p}}{c})$, $d$ is the distance
between two adjacent sensors and $c$ is the propagation velocity of the
acoustic signal, which is assumed to be constant in this paper. The term
$S_{p}(\nu)=a_{p}e_{p}(\nu)$, where $a_{p}$ is the random amplitude of the
$p^{th}$ path and $e_{p}(\nu)$ is the deterministic Fourier Transform of the
emitted signal $e(t)$. ${N_{m}}(\nu)$ is the received white Gaussian noise at
frequency $\nu$. and $\theta_{p}$ is the direction of arrival of the $p^{th}$
raypath at the reference sensor.

Considering the received signals of all the $M$ sensors and constructing the
sparse signal, Eq. 2 can be rewritten in matrix form as follows.

${\bf Y}=\bf{G}{\mathbf{S}}+{\bf{N}}$ (3)

where $\mathbf{Y}(\nu)=[Y_{1}(\nu),Y_{2}(\nu),\ldots,Y_{M}(\nu)]^{T}$ and
$\mathbf{Y}(\nu)$ represents the received signal on all the $M$ sensors, whose
size is $M\times 1$.
$\mathbf{N}(\nu)=[n_{1}(\nu),n_{2}(\nu),\ldots,n_{M}(\nu)]^{T}$ is a vector of
dimension $M\times 1$, which is obtained by the concatenation of the noise at
each sensor.
$\mathbf{G}=[\mathbf{G}_{1},\mathbf{G}_{2},\ldots,\mathbf{G}_{M}]^{T}$, and
$\mathbf{G}_{m}=[e^{-2i\pi\nu\tau_{m1}},\ldots,e^{-2i\pi\nu\tau_{mQ}}]$, which
is an over-complete dictionary. $\mathbf{G}$ contains the terms
$e^{-2i\pi\nu\tau_{mp}},p=1,2,\ldots,P$, which describe the transfer functions
between the source and the sensors.
$\mathbf{S}=[S_{1},S_{2},\cdots,S_{P},\cdots,S_{Q}]^{T}$ denotes the sparse-
signal vector containing the angle and the amplitude information of the
raypaths, which are produced from the emitted signal. $Q$, (${P<M},P<<Q$ and
$M<Q$) is the dimension of the sparse-signal vector, and $T$ means
“transposed”. The so-called sparsity means that the number of nonzero entries
in $\bf S$ is pretty small, i.e., ${\parallel\bf S\parallel}_{0}<Q$.

##  III Reweighed compressive sensing

Based on the signal model in frequency domain described above, we first review
the compressive sensing beamforming algorithm. The signal $\bf S$ is sparse,
as it has only $P$ nonzero components, with $P<<Q$. Thus, the seemingly
underdetermined problem is considered. The sparse vector estimate $\mathbf{S}$
contains the angle and amplitude information of the raypaths, and it can be
estimated by the following $l_{1}$ minimization problem:

$\centering\min_{\mathbf{S}\in\mathbb{C}^{M}}\|\mathbf{S}\|_{1}\textstyle{\
subject\ to}\ \|\mathbf{GS}-\mathbf{Y}\|_{2}\leq\epsilon\@add@centering$ (4)

where the noise is assumed to be spatially white and $\varepsilon$ is an upper
bound for the noise norm. , such that $\|\mathbf{N}\|_{2}\leq\epsilon$.

The reweighed CS Xenaki2014Compressive  improves the performance of the
conventional CS algorithm through the use of a diagonal weight matrix ${\bf
W}_{nn}\in{\mathbb{R}}$ to enhance the sparsity of the solution:

$\centering\min_{\mathbf{S}\in\mathbb{C}^{M}}\|\mathbf{WS}\|_{1}\textstyle{\
subject\ to}\ \|\mathbf{GS}-\mathbf{Y}\|_{2}\leq\epsilon\@add@centering$ (5)

The weight matrix, $\mathbf{W}=\textstyle{diag}[w_{1},w_{2},\cdots,w_{Q}]$, is
initialized with the identity matrix $\mathbf{I}_{Q}$, and the elements
$w_{n}^{k+1}$ at the $(k+1)^{th}$ iteration are updated as

$\centering
w_{n}^{k+1}=\frac{1}{|\hat{s}_{n}^{k}|+{\xi}^{{}^{\prime}}}\@add@centering$
(6)

where
$\widehat{\mathbf{S}}_{n}^{k}=\textstyle{diag}[\hat{s}_{n}^{1},\hat{s}_{n}^{2},\cdots,\hat{s}_{n}^{k}]$
is the solution in the previous iteration. Setting parameter $\xi>0$ ensures
that a null element in the current estimate, $\widehat{\mathbf{S}}_{n}^{k}$,
prevents instability in the iteration. The algorithm will proceed until
solutions in adjacent iterations finally reach approximately equal values;
i.e., $\widehat{\mathbf{S}}_{n}^{k+1}=\widehat{\mathbf{S}}_{n}^{k}$, where

$\centering|w_{n}^{k+1}{s}_{n}^{k+1}|=\begin{cases}\frac{s_{n}^{k}}{|s_{n}^{k}|+\xi}\approx
1\ \ |s_{n}^{k}|>0,\\\ 0,\ \ |s_{n}^{k}|=0\\\ \end{cases}\@add@centering$ (7)

which results in its convergence.

In addition, for stationary sources, which do not move in space, and for the
whole $L$ snapshots, which is a single observation vector of an array data, in
the time domain, due to the consistency of the sparsity the problem is solved
by minimizing the $l_{1}$-norm of the product of the weight $\bf W$ and the
vector ${\bf S}_{l_{1}}$ combined with the $l_{1}$-norm of the row vectors in
$\bf S$,

${\min}{\left\|{\bf WS}_{l_{1}}\right\|_{1}}\text{ subject to }\left\|\bf
GS-Y\right\|_{2}\leq\epsilon$ (8)

where the sparse signal $\bf S$ is a $Q\times L$ matrices and the received
data $\bf Y$ has dimensions $M\times L$. The equation is transformed into the
equation (9),(10).

$w_{q}^{k+1}=\frac{1}{\left|(\hat{s}_{q}^{k})_{l_{1}}\right|+\xi}$ (9)
$\left|w_{q}^{k+1}(s_{q}^{k+1})_{l_{1}}\right|=\left\\{\begin{array}[]{rcl}\frac{(s_{q}^{k})_{l_{1}}}{(s_{q}^{k})_{l_{1}}+\xi}\approx
1,&&\left|(s_{q}^{k})_{l_{1}}\right|>0\\\
0,&&\left|(s_{q}^{k})_{l_{1}}\right|=0\end{array}\right.$ (10)

where $(s_{q}^{k})_{l_{1}}$ is the $q^{th}$ element of the vectors ${\bf
S}_{l_{1}}$ at the $k^{th}$ iteration.

##  IV Subspace-based compressive sensing algorithm

In this paper, the broadband acoustic pulse emission combined to a wavelength
that is small compared to water depth allows us to describe the waveguide
propagation with the ray model. In the case of tank data, the propagation is
dominated by reflections on the top and bottom waveguide boundaries while both
refraction and reflections in shallow waters are observed with the ocean data
used in the present study (Fig. 3c). This means that the $P$ largest
eigenvalues of the spectral covariance matrix can be associated to the $P$
main eigenrays, which may not be the case in oceanic waveguides where modal
dispersion is dominating. We first compute the spectral matrix of the received
signals as follows.

$\hat{\mathbf{R}}=E\\{\mathbf{Y}\mathbf{Y^{H}}\\}$ (11)

where $\hat{\mathbf{R}}\in\mathbb{R}^{M\times M}$. $E\\{\cdot\\}$ denotes the
expectation. The expectation is computed by a frequency-smoothing operation
jiang2017raypath . $\cdot^{\textbf{H}}$ denotes the conjugate transpose. In
addition, the multiple raypaths are produced by the same source, and there
exist correlations between the raypaths, which inevitably lead to rank
deficiency, so we use frequency smoothing for preprocessing. For the wideband
signal, two-sided correlation (TSC) matrices are designed at each frequency.
Each transforms the spectral matrix at a certain frequency into the focusing
spectral matrix at central frequency Valaee1995Wideband . Then, we apply an
eigenvalue decomposition to the spectral matrix $\hat{\mathbf{R}}$. That is,

$\hat{\mathbf{R}}=\hat{\mathbf{R}}_{s}+\hat{\mathbf{R}}_{n}=\sum_{k=1}^{P}\lambda_{k}\mathbf{\mu}_{k}\mathbf{\mu}_{k}^{H}+\sum_{k=P+1}^{M}\lambda_{k}\mathbf{\mu}_{k}\mathbf{\mu}_{k}^{H}$
(12)

where $\lambda_{1},\lambda_{2},\cdots,\lambda_{M}$ are the eigenvalues of
$\hat{\mathbf{R}}$. The eigenvalues are sequentially arranged as
$\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{M}$, and
$\mathbf{\mu}_{1},\cdots,\mathbf{\mu}_{M}$ denote the corresponding
eigenvectors; Therefore, the signal subspace is projected by the largest $P$
eigenvalues (We assume that the number of raypaths $P$ is known in this
paperJiang2014Automatic .) and the corresponding eigenvectors.

$\hat{\mathbf{R}}_{s}=\sum_{k=1}^{P}\lambda_{k}\mathbf{\mu}_{k}\mathbf{\mu}_{k}^{H}$
(13)

$\hat{\mathbf{R}_{s}}$ is rearranged into a long vector
$\hat{\mathbf{r}_{V}}\in\mathbb{R}^{M^{2}\times 1}$ by orderly stacking each
row of $\hat{\mathbf{R}_{s}}$ behind the previous ones. That is,
$\hat{\mathbf{r}}_{V}=[R_{s,1,1},R_{s,1,M},R_{s,M,1}\cdots,R_{s,M,M}]^{T}$,
where $R_{s,i,j}$ denotes the $i^{th}$ row $j^{th}$ column element of
$\hat{\mathbf{R}_{s}}$.

According to Eq. 3, we have

$E\\{\mathbf{YY}^{H}\\}=E\\{\mathbf{GSS}^{H}\mathbf{G}^{H}\\}+E\\{\mathbf{GSN}^{H}\\}+E\\{\mathbf{NS}^{H}\mathbf{G}^{H}\\}+E\\{\mathbf{NN}^{H}\\}$
(14)

Assuming that the signals and the noise are uncorrelated, we obtain

$E\\{\mathbf{GSN}^{H}\\}=E\\{\mathbf{NS}^{H}\mathbf{G}^{H}\\}=\mathbf{0}$ (15)

and the signal subspace is denoted as follows:

$\hat{\mathbf{R}}_{s}=E\\{\mathbf{GSS}^{H}\mathbf{G}^{H}\\}$ (16)

In the context of OAT, the raypaths are produced by the same source. They may
be de-correlated when they scatter at different points on the sea surface or
cross different areas of the water column; however, there is still the
possibility that the case of correlated or coherent raypaths exists. To
further improve separation resolution, it is necessary to consider and
suppress the interferences produced by the correlations between the raypaths.
Thus, the above equation can be rewritten as follows:

$\begin{split}\mathbf{GSS}^{H}\mathbf{G}^{H}==&\mathbf{D}_{1}+\mathbf{D}_{2}\end{split}$
(17)

where

$\begin{split}\mathbf{D}_{1}=\begin{bmatrix}G_{11}&\cdots&G_{1N}\\\
\vdots&\ddots&\vdots\\\
G_{M1}&\cdots&G_{MN}\end{bmatrix}\begin{bmatrix}S_{1}S_{1}^{H}&0&\cdots&0\\\
0&S_{2}S_{2}^{H}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&S_{N}S_{N}^{H}\end{bmatrix}\begin{bmatrix}G_{11}^{H}&\cdots&G_{M1}^{H}\\\
\vdots&\ddots&\vdots\\\ G_{1N}^{H}&\cdots&G_{MN}^{H}\end{bmatrix}\\\
\end{split}$ (18)

and

$\begin{split}\mathbf{D}_{2}=\begin{bmatrix}G_{11}&\cdots&G_{1N}\\\
\vdots&\ddots&\vdots\\\
G_{M1}&\cdots&G_{MN}\end{bmatrix}\begin{bmatrix}0&S_{1}S_{2}^{H}&\cdots&S_{1}S_{N}^{H}\\\
S_{2}S_{1}^{H}&0&\cdots&S_{2}S_{N}^{H}\\\ \vdots&\vdots&\ddots&\vdots\\\
S_{N}S_{1}^{H}&S_{N}S_{2}^{H}&\cdots&0\end{bmatrix}\begin{bmatrix}G_{11}^{H}&\cdots&G_{M1}^{H}\\\
\vdots&\ddots&\vdots\\\ G_{1N}^{H}&\cdots&G_{MN}^{H}\end{bmatrix}\\\
\end{split}$ (19)

Stacking each row of $\mathbf{D}_{1}$ below the previous ones, we rearrange
$\mathbf{D}_{1}$ into a long vector $\mathbf{d_{1V}}$ and further reshape it
as follows:

$\begin{split}\mathbf{d_{1V}}=={\mathbf{G}}^{\prime}\mathbf{S^{{}^{\prime}}}\end{split}$
(20)

where
${\mathbf{G}}^{\prime}=\begin{bmatrix}G_{11}G_{11}^{H}&G_{12}G_{12}^{H}&\cdots&G_{1N}G_{1N}^{H}\\\
G_{11}G_{21}^{H}&G_{12}G_{22}^{H}&\cdots&G_{1N}G_{2N}^{H}\\\ &&\vdots&\\\
G_{M1}G_{M1}^{H}&G_{M2}G_{M2}^{H}&\cdots&G_{MN}G_{MN}^{H}\end{bmatrix}\in\mathbb{C}^{M^{2}\times
N}$ and
$\mathbf{S^{{}^{\prime}}}=\begin{bmatrix}S_{1}S_{1}^{H}&S_{2}S_{2}^{H}&\cdots&S_{N}S_{N}^{H}\end{bmatrix}^{T}\in\mathbb{R}^{N\times
1}$. $\cdot^{T}$ denotes the matrix transpose.

Similarly, we rearrange $\mathbf{D}_{2}$ into a long vector
$\mathbf{{d}_{2V}}$ by stacking each row of $\mathbf{D}_{2}$ below the
previous ones, and $\mathbf{{d}_{2V}}$ will be set as a variable in the
calculation with the convex optimization tools.

Finally, we obtain

$\mathbf{r}_{V}={\mathbf{G}}^{\prime}\mathbf{S^{{}^{\prime}}}+\mathbf{{d}_{2V}}$
(21)

Based on CS theory, an accurate reconstruction is assured under two
conditions: sparsity of the underlying signal and sufficient incoherence of
the process that maps the underlying signals to the observations. In our case,
$Q>>P$ assures sparsity, as there are only $P$ nonzero values. The incoherence
of the process is assured, since the first column to the $M^{th}$ column of
the sensing matrix ${\mathbf{G}}^{\prime}$ constitute a Vandermonde matrix
(The element $a_{ij}$ at $i^{th}$ row and the $j^{th}$ column of a Vandermonde
matrix can be expressed as $a_{ij}=(a_{i})^{j-1}$.). If any two $P$ columns
extracted from the sensing matrix are orthogonal, then the Restricted Isometry
Property (RIP) condition is satisfied. At this time, the rank of the sub-
matrix formed by the extracted arbitrary two $P$ columns is $2P$. As a
Vandermonde matrix satisfies the orthogonality of any two $P$ columns
mentioned above under certain $P$Jafarpour2011Deterministic ;
Cohen2009Compressed , Vandermonde matrixusually satisfies the RIP condition of
compressed sensing.

Thus, the sparse signal can then be reconstructed by solving the following
$L_{1}$ regularization optimization.

$\left\\{\begin{matrix}\hat{\mathbf{S^{{}^{\prime}}}}=\arg\min\left\|\mathbf{S^{{}^{\prime}}}\right\|_{1},\text{
subject to }\mathbf{S^{{}^{\prime}}}>0\\\
\left\|\mathbf{r}_{V}-{\mathbf{G}}^{\prime}\mathbf{S^{{}^{\prime}}}\right\|_{2}\leq\delta,\delta>0\end{matrix}\right.$
(22)

We solved the above programming using available convex optimization tools,
such as CVX Chen2001Atomic . The peak positions of
$\hat{\mathbf{S^{{}^{\prime}}}}$ represent the DOA estimations of the
raypaths.

##  V Performance test

In this section, we illustrate the performance of the subspace-based CS
algorithm with simulation data in a multipath environment, a set of real data
obtained in an ultrasonic waveguideRoux2011Travel  and, finally, a set of
ocean shallow water dataRoux2008The . These tests are all performed in the
same configuration, which is composed of a point source and a vertical array
of receivers. We used $150$, $300$ and $60$ snapshots in the simulations, tank
experiments and at-sea experiments, respectively. The sound speed is uniform
in the experiments, $1500$ m/s for the simulation, $1473$ m/s for tank data,
$1509$ m/s for the real data. We denote the wavelength as $\lambda$, the
height of the water column as $H$, the interval between two adjacent sensors
as $d$ and the distance between the point source and the reference sensor as
$D$. Their specific values corresponding to each test are shown in Table 1. In
order to test the robustness of the subspace-based CS algorithm to the
variation of sensor number and the interval between the two adjacent sensors,
we chose vastly different values of $M$, and $d/\lambda$ for the three cases.

MUSIC and reweighed CS beamforming Xenaki2014Compressive  are provided as
comparative algorithms. Fig. 1 shows the received signals in simulations, tank
experiments and at-sea experiments. Figs. 2, 3 and 4 illustrate the separation
contrast results with simulation data in a multipath environment, a set of
real data obtained in an ultrasonic waveguide and, finally, a set of ocean
shallow water data, respectively. The SNR is defined as the ratio of signal
power to the noise power in the frequency band of the signal for the array. It
is denoted by the equation $SNR=10*lg(P_{s}/P_{n})$, where $P_{s}$ and $P_{n}$
represent the power of signal and noise in the frequency band of the signal,
respectively.

Using different line types, Fig. 2(a) shows the separation results with
simulation data when the SNR is equal to -5 dB, and the black crosses denote
the theoretical positions. The theoretical positions indicate the correct
arrival angles (or arrival times), which are computed using ray theory
Cornuelle1987Compressive ; Skarsoulis2004Travel . It is clear that the
subspace-based CS algorithm can successfully separate all the five raypaths,
while the MUSIC algorithm totally fails in separating them. In addition, the
reweighed CS algorithm detects four peaks, but it is difficult to identify
them by their expected positions. When the SNR is equal to 0 dB, five raypaths
are successfully separated with the subspace-based CS algorithm, and the
corresponding results are shown in Fig. 2(b). In addition, the MUSIC algorithm
does not resolve any raypath. The reweighed CS algorithm provides four peaks,
and each of them is located between two theoretical values. Thus, it is also
difficult to know the corresponding theoretical values of these peaks and
whether these peaks are pseudo-peaks. Moreover, we test the performance of the
subspace-based CS algorithm at different SNRs and compare its separation root-
mean-square errors (RMSE) with those of other algorithms. We define the
average RMSE for the direction of arrival ${\hat{\theta}_{p}}$ for the raypath
$p$ by the quantity:

$RMSE_{\theta_{p}}=\sqrt{\frac{1}{K_{ip}}\sum_{k=1}^{K_{ip}}|{\hat{\theta}_{p}}-\theta_{p}|^{2}},$
(23)

Where $K_{i}$ is the number of the trials, which is equal to 20 for the
simulation. That is, the simulations are performed 20 times for different
SNRs. Generally, the subspace-based CS algorithm is more robust to noise, as
shown in Figs. 2(a) and (b). In Fig. 2(c), even for these raypaths, which can
be separated by all three algorithms, the subspace-based CS algorithm achieves
the minimum RMSEs.

Fig. 4 shows the separation contrast with the data obtained in an ultrasonic
waveguide immersed in a water tank. If the frequency of the signals is
multiplied by a factor and the spatial distances, including both the one
between the source and the receivers and the one between the adjacent
receivers, are divided by the same factor, then the physical phenomena
occurring in the environment remain the same. Thus, the small-scale experiment
can reproduce the actual physical phenomena occurring in nature at a smaller
scale inside the laboratory, which achieves a reduced cost and a totally
controlled experiment. In the present work, a steel bar acts as the bottom,
for which the boundary conditions are nearly perfect at the water-bottom
interface, and a 1.10 m long, 5.4 cm deep acoustic waveguide is constructed
Roux2014Inverting . Two coplanar 64-element vertical line arrays (VLAs) are
placed, and A broadband 1-$\mu$s-long ultrasonic pulse is emitted by each
source successively at the 3 MHz central frequency of the transducers
Roux2011Travel . The transducer dimensions $0.75~{}mm\times 12~{}mm$ are used
to make the linear arrays omnidirectional in the plane defined by the source-
receiver arrays, and the beams are collimated in the plane perpendicular to
the waveguide axis. We chose five point-to-array configurations (one emitted
source and 31 receivers ) and we described their specific positions for the
benchmark configuration in Table I. Five eigenrays are expected between the
source-receiver arrays with this experimental configuration. To further test
the performance of the subspace-based algorithm, we study statistically the
outputs of the beamforming algorithms using a set of sources around the actual
source [the two above and the two below in the source array (SA)] and sets of
subarrays among the vertical received array that are close to the center of
the receiver array (RA). The parameters of the benchmark configuration are
listed in the second line of Table I. In addition, the noise (SNR=-3 dB) has
been added to the received signal for the ultrasound scenario. The same
processing using the subspace-based CS algorithm, the CS beamforming algorithm
and the MUSIC algorithm was performed to produce the beamformer outputs. In
Fig. 3(a), the average beamformer outputs are shown. The purple line
corresponds to the average output of the CS beamforming algorithm; the blue
line with circles corresponds to that of the subspace-based CS algorithm; and
the yellow line with very short vertical lines indicates that of the MUSIC
algorithm. Clearly, the subspace-based CS algorithm correctly separates all
the five raypaths with minimal bias compared to the theoretical position. At
the same time, the MUSIC algorithm fails in separating three of the five
raypaths, and the CS beamforming algorithm gives additional pseudo-peaks. In
Fig. 3(b), only the peaks of these spectra are retained. The purple triangles
denote the peaks of the CS beamforming algorithm; the blue circles represent
those of the subspace-based CS algorithm; and the yellow squares correspond to
those of the MUSIC algorithm. It can be seen that the subspace-based algorithm
shows the smallest dispersion to the theoretical values of the DOAs, which
means that it is more robust to slight changes in the experimental
configuration than the other two algorithms.

Fig. 3 shows the contrasts of separation results with five data sets obtained
in an at-sea experiment Roux2013Analyzing . The FAF05 (Focused Acoustic Field
2005) experiment was conducted in July 2005, to the north of Elba Island,
Italy with repetitive data collection over more than eight consecutive hours
between two source-receiver vertical arrays that were separated by a distance
4.071 km in a 123 m deep waveguide. The FAF05 at-sea experimental setup is
similar to that of the small-scale ultrasonic experiment, although at a much
larger scale. The configuration is composed of two equally spaced VLAs. The
source array (SA) has 29 transducers spanning 78 m in 120 m water, and the
receiver array (RA) has 32 hydrophones covering 62 m. The central frequency of
the transducers is 3.2 kHz, with 1 kHz bandwidth. The source signals were
200-ms linear frequency modulated chirps that were compressed after reception
to their pulse equivalent by cross-correlation. This provided 40 dB of signal-
to-noise ratio for reception with power-limited transmission. Similar to the
case of the tests with tank data, we use five different subarrays that will
provide different experimental configurations with arrival angles that are
close to each other. Fig. 4(c) shows the five eigenrays (including one
surface-reflected ray and four refracted rays) propagating between two
elements of the SA and the RARoux2013Analyzing . The eigenrays are either
refracted at thermoclines or reflected at the air-water interface. All of the
eigenrays that interact with the bottom and the surface-reflected eigenray
have lower amplitudes than the refracted eigenrays because of reflection loss.
The FAF05 data set has been used to develop a coordinated source-receiver
array processing procedure and illustrate its effectiveness in an example of
tracking raylike arrivals in a fluctuating ocean environment Roux2008The . In
addition, using these data, Roux et el. Roux2013Analyzing  analyzed the time-
dependent sound speed fluctuations in shallow water from group-velocity versus
phase-velocity data representation. The subspace-based CS algorithm, the CS
beamforming algorithm and the MUSIC algorithm are applied to the five sets of
data. We calculate the expected DOAs based on the propagation of the raypaths.
The average spectra for the five sets of data, shown in Fig. 4(a), are used to
statistically study the performance of each algorithm. The subspace-based CS
algorithm provides five peaks around the expected theoretical DOAs, while the
CS beamforming algorithm gives multiple peaks with relatively large deviation.
In addition, we only keep the peaks for each algorithm and plot them in Fig.
4(b). It can be seen that there is no dispersion in the results of the
subspace-based CS algorithms and that there is much more dispersion with those
of the CS beamforming algorithm. Additionally, although the results of the
MUSIC algorithm gather together around the expected values, it fails in
detecting some of the raypaths for some subarray configurations. Based on
these results, we can conclude that the subspace-based CS algorithm provides
more reliable results with higher angular resolution.

| $M$ | $d~{}(m)$ | $D~{}(m)$ |  $H$ (m) | $\lambda~{}(m)$ | $d/\lambda$ | $D/H$ | D/$\lambda$
---|---|---|---|---|---|---|---|---
Simulation | 11 | $2.5$ | $2\times 10^{3}$ | 100 | 1 | 2.5 | 20 | $2\times 10^{3}$
| Small-scale
---
experiment
31 |  | $1.5\times$
---
$10^{-3}$
$1$ |  | $5.5\times$
---
$10^{-2}$
| $1.473$
---
$\times 10^{-3}$
1.02 | 18.18 | 678.89
Ocean data | 5 | 2 | $4.701\times 10^{3}$ | 120 | $0.43$ | 4.64 | 39.175 |  | $1.0914$
---
$\times 10^{4}$
Table 1: The configuration parameters of the tests, where $M$ is the number of
sensors, $d$ is the interval between two adjacent sensors, $D$ is the distance
between the point source and the reference sensors, $H$ is the height of the
water column and $\lambda$ is the wavelength.

Figure 1: The received signals: (a) simulations; (b) tank experiments; (c) at-
sea experiments. Figure 2: Separation results comparison of different
algorithms using simulation data. (a) SNR= -5 dB; (b) SNR= 0 dB; (c) RMSE
comparison  Figure 3: (color online) Separation results comparison of
different algorithms in the case of data obtained in a tank (nine sets in
total). A line type (or a mark) corresponds to a result of an algorithm. The
black crosses denote the theoretical values, and the benchmarks are plotted
with dashed lines. (a) Average; (b) Angular peaks. Figure 4:  Separation
results comparison of the different algorithms in the case of ocean data
obtained during the FAF05 at-sea experiment (five sets in total). A linetype
(or a mark) corresponds to a result of an algorithm. The black crosses denote
the theoretical values, and the benchmarks are plotted with dashed lines. (a)
Average; (b) Spectrum peaks; (c) Set of eigenrays propagating between two
elements of the SA and the RA (at a depth of 90 m).

##  VI Conclusions

In this paper, a subspace-based CS algorithm is proposed for separating
raypaths, which introduces the subspace-based algorithm into a CS framework.
Compared to MUSIC and CS beamforming, the proposed algorithm achieved a better
performance in terms of resolution and robustness to ocean fluctuations,
especially when in a strongly noisy environment.

###### Acknowledgements.

The Project Supported by the National Natural Science Foundation of China (No:
61871124 and 61876037), The national defense Pre-Research foundation of China,
by the fund of Science and Technology on Sonar Laboratory (No:
6142109KF201806), by the Stable Supporting Fund of Acoustic Science and
Technology Laboratory (No: JCKYS2019604SSJSSO12). The Focused Acoustic
Forecasting experiment (FAF05) was performed in collaborative experiments with
the NATO Underwater Research Centre (NURC), La Spezia, Italy, with Mark
Stevenson as Chief Scientist. Scientists who contributed to these experiments
include Tuncay Akal, W. A. Kuperman, W. H. Hodgkiss, H. C. Song, B. D.
Cornuelle, Piero Boni, Piero Guerrini, other NURC staff, and the officers and
crew of the RV Alliance.

## References

  * (1) W. Munk, P. Worcester, and C. Wunsch, ” Ocean Acoustic Tomography,” Cambridge, U.K.: Cambridge Univ. Press, 1995\.
  * (2) F. Aulanier, P. Roux, B. Nicolas, R. Brossier, and J. Mars. ”Shallow water acoustic tomography from angle measurements instead of travel-time measurements,” Journal of the Acoustic Society of America, vol. 134, pp. EL373-EL379, 2013.
  * (3) F. Aulanier, B. Nicolas, P. Roux, and J. Mars. ”Time-angle sensitivity kernels for sound-speed perturbations in a shallow ocean,” Journal of the Acoustic Society of America, vol. 134, no. 1, pp. 88-96, 2013.
  * (4) R. Schmidt, ” Multiple emitter location and signal parameter estimation,” IEEE Transactions on Antennas and Propagation, vol. 34, no. 3, pp. 276-280, 1986.
  * (5) L. Jiang, P. Roux, and J. I. Mars, ” Raypath separation with a high-resolution algorithm in a shallow-water waveguide,” IEEE Journal of Oceanic Engineering, vol. 43, no. 1, pp. 119-130, 2018.
  * (6) L. Jiang, Y. Hong, P. Roux, J. Wu, and H. Shu, ” Active wideband higher-order raypath separation in multiple environment,” The Journal of the Acoustical Society of America, vol. 141, no. 1, pp. EL38-EL44, 2017.
  * (7) I. Iturbe, P. Roux, B. Nicolas, and J. Mars. ”Ocean acoustic tomography using a double-beamforming algorithm.” The Journal of the Acoustical Society of America, vol.123, no. 2, pp. 3912, 2008.
  * (8) J. K. Thomas, L. L. Scharf, and D. W. Tufts. ”Probability of a subspace swap in the SVD.” IEEE Transactions on Signal Processing, vol. 43, no. 3, pp. 730-736, 1995.
  * (9) G. F. Edelmann and C. F. Gaumond, ” Beamforming using compressive sensing,” The Journal of the Acoustical Society of America, vol. 130, no. 4, pp. EL232-EL237, 2011.
  * (10) A. Xenaki, P. Gerstoft, and K. Mosegaard, ” Compressive beamforming,” The Journal of the Acoustical Society of America, vol. 136, no. 1, pp. 260-271, 2014.
  * (11) P. Gerstoft, A. Xenaki, and C. F. Mecklenbräuker, ” Multiple and single snapshot compressive beamforming,” The Journal of the Acoustical Society of America, vol. 138, no. 4, pp. 2003-2014, 2015.
  * (12) S. Zhong, Q. Wei, and X. Huang, ” Compressive sensing beamforming based on covariance for acoustic imaging with noisy measurements,” The Journal of the Acoustical Society of America, vol. 134, no. 5, pp. EL445-EL451, 2013.
  * (13) S. Valaee and P. Kabal, ”Wideband array processing using a two-sided correlation transformation,” IEEE Trans. Sign. Process. vol. 43, no. 1, pp.160?172, 1995.
  * (14) L. Jiang, J. I. Mars, ”Automatic detection of the number of raypaths in a shallow-water waveguide”, IEEE Journal of Oceanic Engineering, vol. 39, no. 4, 2014, pp. 713-723.
  * (15) S. Jafarpour. ”Deterministic compressed sensing”. Princeton university, 2011.
  * (16) A. Cohen, W. Dahmen, R. DeVore. ”Compressed sensing and best $k$-term approximation”. Journal of the American mathematical society, vol. 22, no. 1, pp. 211-231, 2009.
  * (17) S. S. Chen, D. L. Donoho, M. A. Saunders, ”Atomic decomposition by basis pursuit”. SIAM review, vol. 43, no. 1, pp. 129-159, 2001.
  * (18) Roux P, Iturbe I, Nicolas B, et al. ”Travel-time tomography in shallow water: experimental demonstration at an ultrasonic scale”. The Journal of the Acoustical Society of America, vol. 130, no. 3, pp. 1232-1241, 2011.
  * (19) Roux P., Cornuelle B. D., Kuperman W. A, WS Hodgkiss ”The structure of raylike arrivals in a shallow-water waveguide”, The Journal of the Acoustical Society of America, vol.124, no. 6, pp. 3430-3439, 2008.
  * (20) B. Cornuelle and B. M. Howe, ”High spatial resolution in vertical slice ocean acoustic tomography,” Journal of Geophysical Research. Atmospheres,vol. 921, no. C11, pp. 11 680-11 692, 1987.
  * (21) E. K. Skarsoulis and B. D. Cornuelle, ”Travel-time sensitivity kernels in ocean acoustic tomography,” The Journal of the Acoustical Society of America, vol. 116, no. 1, pp. 227-238, 2004.
  * (22) Roux P and Nicolas B, ”Inverting for a deterministic surface gravity wave using the sensitivity-kernel approach”, The Journal of the Acoustic Society of America, vol.135, no. 4, pp. 1789-1799, 2014.
  * (23) Roux P., Kuperman W. A, Cornuelle B. D., Aulanier F, Hodgkiss W. S. and Song H. C. ”Analyzing sound speed fluctuations in shallow water from group-velocity versus phase-velocity data representation”, The Journal of the Acoustical Society of America vol.133, no. 3, pp. 1945-1952, 2013.