# NASiam: Efficient Representation Learning using Neural Architecture Search
for Siamese Networks ††thanks: This work was performed using HPC resources
from GENCI-IDRIS (Grant 20XX-AD011012644).

Alexandre Heuillet  IBISC
Université Paris-Saclay, Univ Evry
Evry-Courcouronnes, France
MRL
Massachusetts Institute of Technology
Cambridge, MA, USA
<EMAIL_ADDRESS>Hedi Tabia  IBISC
Université Paris-Saclay, Univ Evry
Evry-Courcouronnes, France
<EMAIL_ADDRESS>Hichem Arioui  IBISC
Université Paris-Saclay, Univ Evry
Evry-Courcouronnes, France
<EMAIL_ADDRESS>
###### Abstract

Siamese networks are one of the most trending methods to achieve self-
supervised visual representation learning (SSL). Since hand labeling is
costly, SSL can play a crucial part by allowing deep learning to train on
large unlabeled datasets. Meanwhile, Neural Architecture Search (NAS) is
becoming increasingly important as a technique to discover novel deep learning
architectures. However, early NAS methods based on reinforcement learning or
evolutionary algorithms suffered from ludicrous computational and memory
costs. In contrast, differentiable NAS, a gradient-based approach, has the
advantage of being much more efficient and has thus retained most of the
attention in the past few years. In this article, we present NASiam, a novel
approach that uses for the first time differentiable NAS to improve the
multilayer perceptron projector and predictor (encoder/predictor pair)
architectures inside siamese-networks-based contrastive learning frameworks
(e.g., SimCLR, SimSiam, and MoCo) while preserving the simplicity of previous
baselines. We crafted a search space designed explicitly for multilayer
perceptrons, inside which we explored several alternatives to the standard
ReLU activation function. We show that these new architectures allow ResNet
backbone convolutional models to learn strong representations efficiently.
NASiam reaches competitive performance in both small-scale (i.e.,
CIFAR-10/CIFAR-100) and large-scale (i.e., ImageNet) image classification
datasets while costing only a few GPU hours. We discuss the composition of the
NAS-discovered architectures and emit hypotheses on why they manage to prevent
collapsing behavior. Our code is available on GitHub.

###### Index Terms:

Deep Learning, NAS, Self-Supervised Learning, Siamese Networks

##  I Introduction

Deep Learning (DL) has experienced rapid growth in the past few years. Two DL
subfields have received much attention: Unsupervised Representation Learning
and Automated Deep Learning (AutoDL).

Unsupervised representation learning aims to make DL models learn strong
representations from unlabeled data. This is especially useful when
considering that data labeling is often a costly and laborious human-made
process. One of the most common approaches to unsupervised visual
representation learning is siamese networks [1]. Siamese networks consist of
two weight-sharing branches (i.e., ”twins”) applied to two or more inputs. The
output feature vectors of the two branches are compared to compute a loss
(e.g., a contrastive loss). In the case of representation learning, the inputs
are usually data augmentations of the same image, and the siamese networks
seek to maximize the similarity between the output feature vectors of the two
branches [2, 3, 4].

On the other hand, AutoDL tries to remove the human factor from the DL
pipeline. Architecture design is one part of this pipeline that has proven
particularly relevant to automate. Most DL architectures are handcrafted and
lack the certainty of an optimal solution [5, 6, 7]. Neural Architecture
Search (NAS) aims to solve this issue by using a meta-learner to search for
neural network architectures relevant to a given task (e.g., image
classification, semantic segmentation, or object detection). NAS algorithms
efficiently browse large search spaces that would prove challenging to
navigate manually. The first NAS works used reinforcement learning [8, 9] or
evolutionary methods [10, 11] but proved particularly inefficient, with
thousands of GPU days needed to obtain a competitive architecture (e.g., 2000
GPU days for NASNet [8]). Nowadays, most approaches use differentiable NAS
[12, 13, 14] (i.e., gradient-based search process) as it requires far less
computational resources and often yields better results.

This article leverages differentiable NAS to discover encoder (projector) and
predictor architectures (i.e., Multilayer Perceptrons) that enable backbone
Convolutional Neural Networks (CNNs) to efficiently learn strong
representations from unlabeled data. To the extent of our knowledge, this is
the first time that NAS has been applied to enhance the architecture of
Siamese networks. Thus, we improved the performance of several siamese network
frameworks such as SimSiam [15], SimCLR [2], or MoCo [4], with an encoder-
predictor pair discovered by a meta-learner inspired by DARTS [16], a popular
differentiable NAS method. We dubbed our approach NASiam (Neural Architecture
Search for Siamese Networks). We show that NASiam reaches competitive results
on small-scale (CIFAR-10, CIFAR-100 [17]) and large-scale (ImageNet [18])
datasets. Thus, Section III highlights our main following contributions:

  * •

A novel way to design encoder/predictor pairs for siamese networks using
differentiable neural architecture search.

  * •

A novel search space specifically designed for the Multi-Layer Perceptron
(MLP) heads of encoder/predictor pairs.

The rest of the article is structured as follows: Section II features a short
survey on related differentiable NAS and siamese networks works. In Section
III, the proposed method is presented. Section IV presents the results of
different image classification experiments and showcases a discussion on the
composition of the discovered encoder/predictor pair architectures, and
Section V brings a conclusion to our work while giving some insights about
future work.

##  II Related Work

This section briefly recalls related work in differentiable NAS, siamese
networks for representation learning, and neural architecture search for
contrastive learning.

###  II-A Differentiable Neural Architecture Search

One of the most trending differentiable NAS family of methods is derived from
Differentiable ARchiTecture Search (DARTS) [16]. This method uses Stochastic
Gradient Descent (SGD) to optimize a set of architectural parameters (denoted
$\alpha$) that represent operations inside building blocks (i.e., elementary
components of the network) called ”cells”. A cell $C$ can be considered as a
direct acyclic graph whose nodes represent states. Its edges are a mix of $K$
different operations $O=\\{o_{1},...,o_{K}\\}$ that define the search space
$S$. Multiple cells can be stacked up to form a global ”supernet” that
encompasses all candidate architectures. Furthermore, as part of a weight-
sharing mechanism, DARTS only searches for two types of cells: normal cells
(i.e., cells that make up most of the network) and reduction cells (i.e.,
cells that perform dimension reduction). The reduction cells are typically
positioned at the 1/3 and 2/3 of the supernet. As the SGD on $\alpha$ occurs
while the supernet containing all cells is being trained on a given dataset,
DARTS is practically solving a bi-level optimization problem. Moreover, the
$\alpha$ weights are relaxed into a continuous form by discretization through
a softmax [19] operation as follows:

$\overline{o}_{i,j}(x)=\sum^{K}_{k=1}\frac{exp(\alpha^{k}_{i,j})}{\sum^{K}_{k^{\prime}=1}exp(\alpha^{k^{\prime}}_{i,j})}o_{k}(x)$
(1)

where $\overline{o}_{i,j}(x)$ is the mixed output of edge $e_{i,j}$ for input
feature $x$ and $\alpha^{k}_{i,j}\in\alpha_{i,j}$ is the weight associated
with operation $o_{k}\in O$ for $e_{i,j}$.

Several works attempted to improve on DARTS. P-DARTS [3] significantly reduced
the search time by progressively deepening the architecture when searching,
leading to a better search space approximation and regularization. PC-DARTS
[20] attempted to reduce DARTS’ memory cost by sampling only a portion of the
supernet to avoid redundancy in the search space exploration. FairDARTS [12]
tried to solve two critical problems that occurred in DARTS, the over-
representation of skip connections and the uncertainty in the probability
distribution of operations. To this end, the authors used the sigmoid function
rather than softmax (see Eq. 1) and crafted a novel loss function that can
push $\alpha$ values towards 0 or 1. DARTS- [21] introduced auxiliary skip
connections that are less prone to become dominant, thus ensuring fairer
competition with the other operations. $\beta$-DARTS [13] introduced a new
regularization method, called Beta-Decay, that prevents the architectural
parameters from saturating. Beta-Decay led to increased robustness and better
generalization ability. Finally, D-DARTS [22] proposed a mechanism to
distribute the search process to the cell level. This approach led to the
individualization of each cell, thus increasing the diversity among the
candidate architectures and expanding the search space.

###  II-B Siamese Neural Networks

Bromley et al. [1] first proposed the Siamese Neural Networks (SNNs)
architecture as twin (i.e., identical and sharing the same weights) models
that process two or more inputs and compare their outputs. The central
intuition behind this concept is that comparing the output feature vectors
will highlight the discrepancies between the inputs. Hence, this approach is
advantageous in signature [1] or face [23] recognition applications.

Another application of SNNs is unsupervised representation learning, also
designated as Self-Supervised Learning (SSL). In particular, it is possible to
learn representations from unlabeled data by feeding variations of the same
input to twin Convolutional Neural Network (CNN) [24] models and computing the
similarity between the output feature vectors. This similarity metric is used
as a loss function, leading the SNNs to learn robust representations (i.e.,
resisting disturbance in the input data). This process is denoted as
contrastive unsupervised learning. Momentum Contrast (MoCo) [25] pre-trains a
CNN using unsupervised learning with a momentum encoder and fine-tunes its
classifier head on standard supervised linear classification. The authors of
MoCo show that the unsupervised pre-trained approach can surpass standard CNN
on multiple ImageNet [18] computer vision tasks. SimCLR [2] added a Multi-
Layer Perceptron (MLP) head as a predictor and highlighted the critical role
of strong data augmentation and large batches (e.g., 8000) in contrastive
learning. Following up on this, [4] proposed an improved version of MoCo
(dubbed MoCo V2) that added a two-layer MLP head in the encoder and modified
the data transforms according to those of SimCLR. Bootstrap Your Own Latent
(BYOL) [26] proposed an SNN framework centered around an online network and a
target network. The output of the target network is iteratively bootstrapped
to serve as input to the online network. The authors showed that BYOL could
learn more robust representations than previous approaches. Finally, SimSiam
[15] introduced a simpler SNN architecture that removes the need for negative
sample pairs, momentum encoders, and large batches. More specifically, SimSiam
implements a stop-grad mechanism that stops gradient backpropagation in one of
the two branches of the twin model. Despite being a more straightforward
approach than previous baselines, SimSiam reaches a competitive score on
ImageNet classification.

###  II-C Neural Architecture Search for Contrastive Self-Supervised Learning

A handful of previous works have already explored using NAS for contrastive
Self-Supervised Learning (SSL). [27] first introduced a method to leverage NAS
to improve existing SSL frameworks such as SimCLR [2]. Their approach, dubbed
SSNAS, is derived from DARTS [16] and reached competitive performance compared
with supervised models. [28] proposed CSNAS, a novel way to search for SSL-
focused CNN architectures using Sequential Model-Based Optimization. CSNAS
leverages a cell-based search space similar to DARTS [16] and performs
contrastive SSL using PIRL [29]. The authors showed that CSNAS managed to
overperform or match both handcrafted architectures and supervised NAS models
on image classification tasks. Another work of note is SSWP-NAS [30], a proxy-
free weight-preserving NAS method for SSL. Similarly to CSNAS, SSWP-NAS is
based on DARTS and navigates through a cell-based search space to discover new
CNN architectures. SSWP-NAS overperformed previous SSL NAS methods and reached
competitive results compared to supervised NAS approaches. In a drastically
different approach, Contrastive Neural Architecture Search (CTNAS) [31]
refactors NAS with Contrastive Learning. A Neural Architecture Comparator is
designed to drive the search process by comparing candidate architectures with
a baseline architecture. Thus, in this approach, contrary to other works,
Contrastive Learning is used to enhance NAS rather than the other way around.

In this article, we propose to go further than the previous works listed above
by using differentiable NAS to directly enhance the Siamese (i.e., MLP)
architecture rather than improve the backbone CNN (which is similar to what
trending NAS frameworks such as DARTS [16] or FBNet [14] do). In Section IV,
we show that our NASiam approach is able to discover novel Siamese
architectures reaching higher performance than standard Contrastive SSL
frameworks such as SimCLR [2] or MoCo [4].

##  III Proposed Approach

This section highlights the key ideas behind our proposed approach: searching
for the multi-layer perceptron components of the encoder/predictor pair and
crafting an original search space specific to contrastive learning with
siamese neural networks.

###  III-A Searching for an Encoder/Predictor Pair

First, we focused on SimSiam [15] as a simple baseline upon which to build our
approach. SimSiam uses a siamese architecture consisting of an encoder $f$ and
a predictor $h$. The encoder $f$ is composed of a baseline CNN (e.g., ResNet50
[7]) and of a projector head (i.e., a three-layer MLP) that is duplicated on
twin branches that take variations of the same image as input. A two-layer MLP
$h$ is then added on top of one of the branches to act as a predictor head.
The discrepancy between the output feature vectors of the two branches is
computed using a contrastive loss as follows:

$\mathcal{L}=\frac{1}{2}(\mathcal{D}(p_{1},\texttt{stopgrad}(z_{2}))+\mathcal{D}(p_{2},\texttt{stopgrad}(z_{1}))$
(2)

where $z_{1}=f(x_{1})$, $z_{2}=f(x_{2})$, $p_{1}=h(z_{1})$, $p_{2}=h(z_{2})$
for input images $x_{1}$ and $x_{2}$, stopgrad is a mechanism that stops
gradient backpropagation (in other words, the argument inside stopgrad is
detached from the gradient computation), and $\mathcal{D}$ is the negative
cosine similarity defined as follows:

$\mathcal{D}(p,z)=-\frac{p}{||p||_{2}}.\frac{z}{||z||_{2}}$ (3)

where $||.||_{2}$ is the $l_{2}$ norm.

In our proposed approach, we kept most of the global structure of the
underlying siamese framework. However, we used a Differentiable NAS method to
search for an encoder projector head architecture up to $n$ layers and a
predictor architecture up to $m$ layers. More specifically, we consider a set
$O=\\{o_{1},...,o_{K}\\}$ of candidate operations. We search for two cells
(see Section II-A) $C_{e}$ and $C_{p}$ for the encoder and decoder
respectively. Contrary to DARTS [16], each cell is structured as a linear
sequence of layers where each layer is a mixed output of $|O|=K$ operations.
Each operation $o$ in each layer $i$ is weighted by a parameter
$\alpha_{i}^{o}$. The sets of architectural parameters for $C_{e}$ and $C_{p}$
are denoted $\alpha_{e}$ and $\alpha_{p}$ respectively. Similarly to Eq. 1,
operation values in each layer are discretized as follows:

$\overline{o}_{i}(x)=\sum^{K}_{k=1}\sigma_{SM}(\alpha^{k}_{i})o_{k}(x)$ (4)

where $\overline{o}_{i}$ is the mixed operation of layer $i$, $\alpha^{k}_{i}$
is the architectural weight assigned to $o_{k}\in O$ for layer $i$, and
$\sigma_{SM}$ denotes the softmax operation. The supernet encompassing $f$ and
$h$ is trained on a portion of a dataset while $C_{e}$ and $C_{p}$ are
simultaneously searched on another portion of the same dataset. Hence, we
solve a bi-level optimization problem formulated as

$\begin{split}\underset{\alpha_{e},\alpha_{p}}{\text{min
}}\mathcal{L}_{val}(w^{*}(\alpha_{e},\alpha_{p}),\alpha_{e},\alpha_{p}),\\\
\text{s.t.}w^{*}(\alpha_{e},\alpha_{p})=\underset{w}{\text{argmin
}}\mathcal{L}_{train}(w,\alpha_{e},\alpha_{p}),\end{split}$ (5)

where $w$ denotes the supernet weights,
$\mathcal{L}_{train}(w,\alpha_{e},\alpha_{p})=\mathcal{L}(w,\alpha_{e},\alpha_{p})$
is the training loss, and
$\mathcal{L}_{val}(w^{*},\alpha_{e},\alpha_{p})=\mathcal{L}(w^{*},\alpha_{e},\alpha_{p})$
is the validation loss.

Once the search phase is complete, for each layer $i$ of each cell, we select
the best-performing operation according to the discretized weights
$\alpha_{e}$ and $\alpha_{p}$ to form the encoder/predictor architecture
genotype $G$. The whole neural architecture search process is detailed in
Algorithm 1.

Algorithm 1 Algorithm describing the differentiable neural architecture search
process of NASiam

0: Object: $C_{e}$, encoder cell containing architectural weights

0: Object: $C_{p}$, predictor cell containing architectural weights

0: List: $O$, list of operations

0: List: $D_{t}$, train dataset

0: List: $D_{v}$, validation dataset

0: Object: $model$, backbone CNN model

0: Object: $opt$, model optimizer

0: Object: $search\\_opt$, search optimizer

0: Integer: $E$, number of epochs

for $e$ in [0, $E$[ do

for $(x_{1},x_{2})$ in $D_{v}$ do

$(x_{1},x_{2})\leftarrow model(x_{1},x_{2})$

$(z_{1},z_{2})\leftarrow C_{e}(x_{1},x_{2})$

$(p_{1},p_{2})\leftarrow C_{p}(x_{1},x_{2})$

$\text{stop\\_grad}(z_{1},z_{2})$

$loss\leftarrow-0.5(\text{cosine\\_similarity}(p_{1},z_{2})+\text{cosine\\_similarity}(p_{2},z_{1}))$

$search\\_opt.\text{optimization\\_step}(loss,C_{e}.weights,C_{p}.weights)$

end for

for $(x_{1},x_{2})$ in $D_{t}$ do

$(x_{1},x_{2})\leftarrow model(x_{1},x_{2})$

$(z_{1},z_{2})\leftarrow C_{e}(x_{1},x_{2})$

$(p_{1},p_{2})\leftarrow C_{p}(x_{1},x_{2})$

$\text{stop\\_grad}(z_{1},z_{2})$

$loss\leftarrow-0.5(\text{cosine\\_similarity}(p_{1},z_{2})+\text{cosine\\_similarity}(p_{2},z_{1}))$

$opt.\text{optimization\\_step}(loss,C_{e}.weights,C_{p}.weights)$

end for

end for

$A\leftarrow\text{empty\\_list()}$

for $C$ in $\\{C_{e},C_{p}\\}$ do

$C_{f}\leftarrow\text{empty\\_list()}$

$n\leftarrow|C|$

for $i$ in $[0,n[$ do

$op\leftarrow\text{argmax}(C.weights[i])$

$\text{append}(op,C_{f})$

end for

$\text{append}(C_{f},A)$

end for

return $A$

Our approach, dubbed NASiam (Neural Architecture Search for Siamese Networks),
is summarized in Fig. 1. In addition to SimSiam, we also experimented NASiam
on other Siamese frameworks such as SimCLR[2] and MoCo V2 [4]. However, those
frameworks (SimCLR and MoCo V2) do not rely on a predictor. Hence, in that
case, we only performed NAS for the MLP projector head of the encoder (i.e.,
only searching for cell $C_{e}$).

In Section IV, we show that NASiam can consistently improve the performance of
popular siamese frameworks (SimSiam, SimCLR, MoCo, and BYOL) in both small-
scale (CIFAR-10 and CIFAR-100 [17]) and large-scale (ImageNet [18]) image
classification datasets.

Figure 1: Layout of the NASiam architecture. Siamese network encoder/predictor
(projection MLPs) architectures are searched using differentiable NAS wrapped
around a Siamese framework such as SimSiam [2], which is the baseline used in
the present figure. First, an input image x is augmented to produce two
variations $x_{1}$ and $x_{2}$. Each of these two inputs is then fed into one
of the two branches of the siamese network. While both $x_{1}$ and $x_{2}$ go
through an encoder equipped with an MLP projection head, $x_{2}$ is further
processed by an MLP predictor. Finally, a negative cosine contrastive loss is
computed and backpropagated to minimize the similarity between the two
branches’ output feature maps. Both the encoder and the decoder contain cells
(i.e., $C_{e}$ and $C_{p}$, respectively) that are designed using a
differentiable NAS approach. Architectural parameters for $C_{e}$ and $C_{p}$
are denoted $\alpha_{e}$ and $\alpha_{p}$ respectively.

###  III-B Crafting a Contrastive Learning-Specific Search Space

To accompany our novel NASiam approach (see Section III-A), we crafted an
original search space $S$ specifically designed for MLPs. $S$ comprises the
following 7 operation blocks: linear \+ batch_norm \+ ReLU, linear \+
batch_norm \+ Hardswish, linear \+ batch_norm \+ SiLU, linear \+ batch_norm \+
ELU, max_pool_3x3 (1-dimensional) + batch_norm, avg_pool_3x3 (1-dimensional) +
batch_norm, and Identity (skip connection). Hence, $S$ includes several types
of fully connected layers, each featuring a different activation function. The
motivation behind adding activation functions to the search space is to
increase diversity among the candidate architectures and explore alternatives
to the classic ReLU function (e.g., Hardswish [32], or Mish [33]). To that
end, it also makes it possible to mix different activation functions according
to the type of network (i.e., projector or predictor) and the location inside
that network. In contrast, previous baselines [2, 15, 4] only relied on a
single activation for both networks regardless of their respective
architectures.

While unconventional, including pooling layers in the search space is helpful,
as we show in Section IV that they can help prevent collapsing. Moreover, the
authors of SimSiam [15] indicated that insufficient or too many Batch
Normalization (BN) layers could cause the model to underperform severely or
become unstable. They empirically demonstrate that the optimal setting for
SimSiam is to place BNs after every layer except for the predictor’s output
layer. Hence, we follow this assertion by adding BNs after every linear and
pooling operation except for the predictor’s final layer. Finally, we also
included the identity operation so that the search algorithm can modulate the
number of layers in the architecture. This way, we can indicate a maximum
number of layers $n$, and the search algorithm can craft an architecture of
size $m<n$ by skipping layers.

##  IV Experiments

This section presents the results of our image classification experiments on
small-scale (CIFAR-10, CIFAR-100) and large-scale (ImageNet) datasets.

###  IV-A Experimental Settings

We used RTX 3090 and Tesla V100 Nvidia GPUs to conduct our experiments. We
searched for predictor/encoder pairs for 100 epochs on CIFAR-10, and CIFAR-100
[17] using the SGD optimizer with $lr=0.06$, $wd=5e-4$, and a batch size of
512. We set a maximum of 6 layers for the encoder. If the baseline siamese
framework relies on a predictor, we search for a 4-layer predictor
architecture. The whole search process on these settings takes around 2.3 GPU
hours on a single GPU. We did not search on the full ImageNet [18] dataset as
it is prohibitively expensive (i.e., it takes around 12 GPU days on a single
GPU). Instead, we transferred our best CIFAR-searched architecture to
ImageNet. For the pre-training and linear classification phases, we kept the
same settings as [15]. Our code is based on PyTorch 1.12.

###  IV-B Ablation Study on the Importance of Pooling Layers

We conducted an ablation study on the importance of including pooling layers
in our novel space search $S$ (see Section III-B). To this end, we simply
removed max_pool_3x3 and avg_pool_3x3 from $S$ to form $S^{\prime}$. When
comparing the results in Table I, we can observe that, when searching on
$S^{\prime}$ rather than $S$, the validation top-1 accuracy of NASiam drops
significantly (by around 3 %). Moreover, when analyzing the genotypes searched
on CIFAR-10 and CIFAR-100 using $S$, it appears that the predictor
architectures always contain pooling layers (making up to 40 % of the total
architecture). In addition, Fig. 2 shows that the model (see Eq. 2) achieved
better similarity and faster convergence when searched on $S$ rather than
$S^{\prime}$. Thus, these findings highlight the critical role pooling layers
play in ensuring high performance and preventing collapse, especially
concerning the encoder architecture.

TABLE I: Results on CIFAR-10 linear classification of two NASiam models using search space $S$ and $S^{\prime}$ respectively. Both models were pre-trained for 100 epochs. The baseline framework is SimSiam with a ResNet18 backbone. Search Space |  Search Epochs |  Pre-Train Epochs |  Validation Top-1 (%) |  Validation Top-5
---|---|---|---|---
$S$ | 50 | 100 | 66.6 | 91.3
$S^{\prime}$ | 50 | 100 | 63.7 | 86.1
Figure 2: Plot of the negative cosine contrastive loss while pretraining two
NASiam models on CIFAR-10. The baseline framework is SimSiam with a ResNet18
backbone. The two models are searched on search spaces $S$ (blue line) and
$S^{\prime}$ (red line) respectively. The model searched on $S$ achieves
better similarity, thus making the relevance of pooling layers clear.

###  IV-C Incidence of Data Augmentations on the NAS process

In self-supervised learning, data augmentations are paramount to prevent the
model from overfitting and the contrastive loss (see Eq. 2) from saturating to
-1. In contrast, differentiable neural architecture search methods [16, 12,
22] scarcely employ data augmentation as they only train the supernet for a
small number of epochs (e.g., 50). Thus, a legitimate interrogation is how the
strong data augmentation policy used in SSL frameworks can interfere with the
differentiable search process. To answer this question, we searched for two
different SimSiam [15] models with and without the data augmentation policy
activated on CIFAR-10 and compared the resulting architectures.

Table II shows that deactivating the data augmentation policy leads to a
degenerated architecture with a dominance of skip connections (50 % of the
architecture) associated with performance collapse. Furthermore, Fig. shows
that, during the search phase, the similarity loss converges significantly
faster towards -1, thus presenting a collapsing behavior. This observation
correlates with the architectural collapse described in numerous
differentiable NAS studies [12, 34, 13]. This collapsing behavior is akin to
overfitting for NAS and is caused by the high prominence of skip connections
due to their unfair advantage (compared to parametric operations). Thus, data
augmentation clearly has a positive impact on the differentiable search
process and should not be deactivated, in contrast with supervised learning.

TABLE II: Results on CIFAR-10 linear classification of two NASiam models using either SimSiam data augmentation policy or no data augmentation. Both models were pre-trained for 800 epochs. The baseline framework is SimSiam with a ResNet18 backbone. Data Augmentation |  Search Epochs |  Pre-Train Epochs |  Validation Top-1 (%) |  % of Skip Connections
---|---|---|---|---
Yes | 100 | 800 | 91.2 | 20
No | 100 | 800 | 10.0 | 50
Figure 3: Plot of the negative cosine contrastive loss while pretraining two
NASiam models on CIFAR-10. The baseline framework is SimSiam with a ResNet18
backbone. The two models are searched with and without data augmentation
respectively.

###  IV-D Preliminary Results on CIFAR

To quickly assess the behavior of our novel approach NASiam, we first
conducted preliminary experiments on small-scale CIFAR datasets [17]. We
searched NASiam architectures for 100 epochs on CIFAR-10 and CIFAR-100 using
the CIFAR version of ResNet18 [7] as the encoder backbone. Then, we performed
unsupervised pretraining for 800 epochs with a cosine annealing schedule
before training a linear classifier using frozen features for 100 epochs. In
these settings, NASiam overperforms SimSiam by 1.4 % and 0.4 % on CIFAR-10 and
CIFAR-100 respectively (see Table III and Table IV). In addition, Fig. 4 shows
us that NASiam can achieve better similarity than SimSiam without saturating
the contrastive loss to $-1$ (i.e., a collapsing behavior). Furthermore,
results were also positive when using alternative siamese frameworks, with
NASiam overperforming both MoCo V2 [4] and SimCLR [2].

TABLE III: Results of pre-training for 800 epochs on CIFAR-10 linear classification with SGD. The backbone is the CIFAR version of ResNet18. $\dagger$: Result obtained by running the official implementation with the hyperparameters suggested by the authors for CIFAR-10. Model |  Batch Size |  Pre-Train Epochs |  Train Epochs |  Validation Top-1 (%)
---|---|---|---|---
MoCo V2 [4]$\dagger$ | 256 | 800 | 100 | 89.8
SimCLR [2]$\dagger$ | 256 | 800 | 100 | 91.1
SimSiam [15]$\dagger$ | 256 | 800 | 100 | 89.5
Ours NASiam (SimSiam) | 256 | 800 | 100 | 91.2
NASiam (MoCo V2) | 256 | 800 | 100 | 90.4
NASiam (SimCLR) | 256 | 800 | 100 | 92.1
TABLE IV: Results of training for 800 epochs on CIFAR-100 linear classification with SGD. The backbone is the CIFAR version of ResNet18. $\dagger$: Result obtained by running the official implementation with the hyperparameters suggested by the authors for CIFAR-10. Model |  Batch Size |  Pre-Train Epochs |  Train Epochs |  Validation Top-1 (%)
---|---|---|---|---
MoCo V2 [4]$\dagger$ | 256 | 800 | 100 | 62.9
SimCLR [2]$\dagger$ | 256 | 800 | 100 | 63.6
SimSiam [15]$\dagger$ | 256 | 800 | 100 | 63.7
Ours NASiam (SimSiam) | 256 | 800 | 100 | 64.1
NASiam (MoCo V2) | 256 | 800 | 100 | 65.0
NASiam (SimCLR) | 256 | 800 | 100 | 68.9
Figure 4: Plot of the negative cosine contrastive loss when pretraining
SimSiam and NASiam for 800 epochs on CIFAR-10. NASiam converges faster without
collapsing and achieves better similarity than SimSiam.

###  IV-E Results on ImageNet

We conducted image classification experiments on ImageNet [18] as a standard
practice to evaluate the performance of our novel approach on large-scale
datasets. As stated in Section IV-A, we transferred our best CIFAR
architecture instead of searching directly on ImageNet to save computational
resources. Then, we performed unsupervised pretraining on ImageNet for 100
epochs before training a linear classifier with frozen features for 100
epochs. The results are presented in detail in Table V. As for CIFAR (see
IV-D), NASiam consistently achieves better linear classification results than
the baseline frameworks, thus validating its usefulness.

TABLE V: Results of training for 100 epochs on ImageNet linear classification with SGD. The backbone is ResNet50. Models were pre-trained for 100 epochs on ImageNet. $\dagger$: Score obtained by the authors of SimSiam by using an improved version of the model. Model |  Batch Size |  Pre-Train Epochs |  Train Epochs |  Validation Top-1 (%)
---|---|---|---|---
MoCo V2 [4]$\dagger$ | 256 | 100 | 100 | 67.4
BYOL [26]$\dagger$ | 4096 | 100 | 100 | 66.5
SimCLR [2]$\dagger$ | 4096 | 100 | 100 | 66.5
SimSiam [15] | 256 | 100 | 100 | 67.1
Ours NASiam (SimSiam) | 256 | 100 | 100 | 67.4
NASiam (MoCo V2) | 256 | 100 | 100 | 67.4
NASiam (SimCLR) | 4096 | 100 | 100 | 67.2

###  IV-F Object Detection and Instance Segmentation Results on COCO

Table VI displays the results of transferring our NASiam models pretrained on
ImageNet [18] to Microsoft COCO [35] object detection and instance
segmentation tasks. We can see that NASiam consistently overperforms
handcrafted SSL architectures in both tasks. Hence, NASiam architectures can
successfully generalize to computer vision tasks other than image
classification.

TABLE VI: Comparison of backbone models for MaskRCNN [36] on COCO [35] using a 1x schedule and ResNet50 [7] as the baseline CNN. All models are pretrained for 200 epochs on ImageNet, finetuned for 12 epochs on COCO 2017 train set, and evaluated on COCO 2017 val set. Models |  $AP_{50}$ (%) |  $AP$ (%) |  $AP_{75}$ (%) |  $AP^{mask}_{50}$ (%) |  $AP^{mask}$ (%) |  $AP^{mask}_{75}$ (%)
---|---|---|---|---|---|---
ImageNet supervised | 58.2 | 38.2 | 41.2 | 54.7 | 33.3 | 35.2
SimCLR | 57.7 | 37.9 | 40.9 | 54.6 | 33.3 | 35.3
SimSiam | 57.5 | 37.9 | 40.9 | 54.2 | 33.2 | 35.2
MoCo V2 | 58.8 | 39.2 | 42.5 | 55.5 | 34.3 | 36.6
BYOL | 57.8 | 37.9 | 40.9 | 54.3 | 33.2 | 35.0
NASiam (SimSiam) | 58.6 | 39.0 | 42.1 | 55.2 | 34.1 | 36.3

###  IV-G Discussion on the Composition of the Architectures

Some facts are noteworthy when comparing encoder/predictor architectures
discovered on CIFAR-10 by our novel approach (see Section III-A) with those of
SimSiam [15].

First, in Fig. 5, we can see that both ResNet50 and ResNet18[17] NASiam
architectures are significantly deeper than the original SimSiam architecture.
Furthermore, a remarkable fact is that the ReLU activation function is in
minority in the discovered architectures (and even disappeared completely from
the ResNet50 one). Instead, a mix of different activation functions is
preferred, with SiLU and Hardswish having a high prominence. Thus, this may
indicate that ReLU, despite its popularity, is not the optimal activation
function for performing contrastive learning. In addition, the optimizer
always selected at least one AvgPool3x3+BN layer to be part of the predictor
architecture, hence validating the relevance of including pooling layers in
the search space (as already highlighted in Section IV-B).

Finally, when comparing both NAS-discovered architectures, we can observe that
the ResNet50 one possesses a deeper encoder than the ResNet18-based
architecture (i.e., 6 vs. 4 layers), with additional Linear+BN+Swish and
Linear+BN+SiLU blocks. However, the two predictor architectures retain the
same depth and a similar composition. This is coherent with the
recommendations of the authors of SimSiam [15], where they selected a
shallower architecture when training on CIFAR-10 with ResNet18 rather than
ResNet50. One hypothesis to explain this discrepancy in architectural sparsity
is that ResNet18, being a shallower model than ResNet50, has a less powerful
innate ability to extract representations and hence produces less complex
feature maps that would not require a deep projector to be analyzed. Using a
deeper architecture could even lead to adverse effects. To confirm this
hypothesis, we tried to fit a ResNet18 model on CIFAR-10 with the deeper
encoder/predictor pair discovered for ResNet50. Fig. 6 clearly shows that this
architectural setting quickly led to a collapsing behavior (with the
contrastive loss rapidly saturating to -1 as soon as epoch 350) with a higher
variance than the ResNet18-searched architecture. Hence, this validates the
ability of our NASiam approach to discover backbone-specific architectures.

Figure 5: Composition of encoder/predictor pair architectures. (Top) SimSiam
model. (Bottom left) NASiam model searched for 100 epochs on CIFAR-10 using
SimSiam as the baseline framework with ResNet18 as the backbone CNN. (Bottom
right) NASiam model searched for 100 epochs on CIFAR-10 using SimSiam as the
baseline framework with ResNet50 as the backbone CNN. ResNet18-searched and
ResNet50-searched architectures are clearly different, with ResNet50 needing a
deeper encoder. Figure 6: Plot of the negative cosine similarity loss while
pretraining NASiam with ResNet18 using architectures searched either with
ResNet18 or ResNet50 as backbone. The ResNet50-searched architecture quickly
collapses towards -1 and has high variance while the ResNet18-searched one
converges as expected.

##  V Conclusion

In this article, we presented NASiam, a novel approach for contrastive
learning with siamese networks that searches for efficient encoder/predictor
pairs using differentiable neural architecture search (see Section III-A).
This universal method can enhance many existing siamese frameworks while
preserving their underlying structure. In addition, NASiam is efficient as it
only costs a few GPU hours. Section IV showed that NASiam discovers
encoder/predictor pair architectures that efficiently learn robust
representations and overperform previous baselines in small-scale and large-
scale image classification datasets. These empirical results support our
intuition that the encoder and predictor architectural designs play a decisive
role in representation learning. We hope this work will pave the way to
further improvements for MLP-headed siamese networks.

## References

  * [1] J. Bromley, I. Guyon, Y. LeCun, E. Säckinger, and R. Shah, “Signature verification using a” siamese” time delay neural network,” Advances in neural information processing systems, vol. 6, 1993.
  * [2] T. Chen, S. Kornblith, M. Norouzi, and G. Hinton, “A simple framework for contrastive learning of visual representations,” in International conference on machine learning, pp. 1597–1607, PMLR, 2020.
  * [3] X. Chen, L. Xie, J. Wu, and Q. Tian, “Progressive darts: Bridging the optimization gap for nas in the wild,” International Journal of Computer Vision, vol. 129, no. 3, pp. 638–655, 2021.
  * [4] X. Chen, H. Fan, R. Girshick, and K. He, “Improved baselines with momentum contrastive learning,” arXiv preprint arXiv:2003.04297, 2020.
  * [5] F. Chollet, “Xception: Deep learning with depthwise separable convolutions,” in Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 1251–1258, 2017.
  * [6] C. Szegedy, S. Ioffe, V. Vanhoucke, and A. A. Alemi, “Inception-v4, inception-resnet and the impact of residual connections on learning,” in Thirty-first AAAI conference on artificial intelligence, 2017.
  * [7] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 770–778, 2016.
  * [8] B. Zoph and Q. V. Le, “Neural architecture search with reinforcement learning,” in International Conference on Learning Representations, 2017\.
  * [9] B. Zoph, V. Vasudevan, J. Shlens, and Q. V. Le, “Learning transferable architectures for scalable image recognition,” in Proceedings of the IEEE conference on Computer Vision and Pattern Recognition, pp. 8697–8710, 2018\.
  * [10] Y. Chen, G. Meng, Q. Zhang, S. Xiang, C. Huang, L. Mu, and X. Wang, “Renas: Reinforced evolutionary neural architecture search,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 4787–4796, 2019.
  * [11] E. Real, A. Aggarwal, Y. Huang, and Q. V. Le, “Regularized evolution for image classifier architecture search,” in Proceedings of the aaai conference on artificial intelligence, vol. 33, pp. 4780–4789, 2019.
  * [12] X. Chu, T. Zhou, B. Zhang, and J. Li, “Fair darts: Eliminating unfair advantages in differentiable architecture search,” in European conference on computer vision, pp. 465–480, Springer, 2020.
  * [13] P. Ye, B. Li, Y. Li, T. Chen, J. Fan, and W. Ouyang, “$\beta$-darts: Beta-decay regularization for differentiable architecture search,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2022.
  * [14] X. Dai, A. Wan, P. Zhang, B. Wu, Z. He, Z. Wei, K. Chen, Y. Tian, M. Yu, P. Vajda, et al., “Fbnetv3: Joint architecture-recipe search using predictor pretraining,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 16276–16285, 2021.
  * [15] X. Chen and K. He, “Exploring simple siamese representation learning,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 15750–15758, 2021.
  * [16] H. Liu, K. Simonyan, and Y. Yang, “DARTS: Differentiable architecture search,” in International Conference on Learning Representations, 2019\.
  * [17] A. Krizhevsky, G. Hinton, et al., “Learning multiple layers of features from tiny images,” 2009.
  * [18] O. Russakovsky, J. Deng, H. Su, J. Krause, S. Satheesh, S. Ma, Z. Huang, A. Karpathy, A. Khosla, M. Bernstein, A. C. Berg, and L. Fei-Fei, “ImageNet Large Scale Visual Recognition Challenge,” International Journal of Computer Vision (IJCV), vol. 115, no. 3, pp. 211–252, 2015.
  * [19] J. S. Bridle, “Probabilistic interpretation of feedforward classification network outputs, with relationships to statistical pattern recognition,” in Neurocomputing (F. F. Soulié and J. Hérault, eds.), (Berlin, Heidelberg), pp. 227–236, Springer Berlin Heidelberg, 1990.
  * [20] Y. Xu, L. Xie, X. Zhang, X. Chen, G.-J. Qi, Q. Tian, and H. Xiong, “Pc-darts: Partial channel connections for memory-efficient architecture search,” arXiv preprint arXiv:1907.05737, 2019.
  * [21] X. Chu, X. Wang, B. Zhang, S. Lu, X. Wei, and J. Yan, “Darts-: Robustly stepping out of performance collapse without indicators,” in International Conference on Learning Representations, 2021.
  * [22] A. Heuillet, H. Tabia, H. Arioui, and K. Youcef-Toumi, “D-darts: Distributed differentiable architecture search,” arXiv preprint arXiv:2108.09306, 2021\.
  * [23] Y. Taigman, M. Yang, M. Ranzato, and L. Wolf, “Deepface: Closing the gap to human-level performance in face verification,” in Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 1701–1708, 2014\.
  * [24] Y. LeCun, Y. Bengio, et al., “Convolutional networks for images, speech, and time series,” The handbook of brain theory and neural networks, vol. 3361, no. 10, p. 1995, 1995.
  * [25] K. He, H. Fan, Y. Wu, S. Xie, and R. Girshick, “Momentum contrast for unsupervised visual representation learning,” in Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp. 9729–9738, 2020.
  * [26] J.-B. Grill, F. Strub, F. Altché, C. Tallec, P. Richemond, E. Buchatskaya, C. Doersch, B. Avila Pires, Z. Guo, M. Gheshlaghi Azar, et al., “Bootstrap your own latent-a new approach to self-supervised learning,” Advances in Neural Information Processing Systems, vol. 33, pp. 21271–21284, 2020.
  * [27] S. Kaplan and R. Giryes, “Self-supervised neural architecture search,” arXiv preprint arXiv:2007.01500, 2020.
  * [28] N. Nguyen and J. M. Chang, “Csnas: Contrastive self-supervised learning neural architecture search via sequential model-based optimization,” IEEE Transactions on Artificial Intelligence, vol. 3, no. 4, pp. 609–624, 2021.
  * [29] I. Misra and L. van der Maaten, “Self-supervised learning of pretext-invariant representations,” in IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2020.
  * [30] Z. Li, Y. Gao, Z. Zha, Z. Hu, Q. Xia, S. Zhang, and D. N. Metaxas, “Towards self-supervised and weight-preserving neural architecture search,” arXiv preprint arXiv:2206.04125, 2022.
  * [31] Y. Chen, Y. Guo, Q. Chen, M. Li, W. Zeng, Y. Wang, and M. Tan, “Contrastive neural architecture search with neural architecture comparators,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 9502–9511, 2021.
  * [32] A. Howard, M. Sandler, G. Chu, L.-C. Chen, B. Chen, M. Tan, W. Wang, Y. Zhu, R. Pang, V. Vasudevan, et al., “Searching for mobilenetv3,” in Proceedings of the IEEE/CVF international conference on computer vision, pp. 1314–1324, 2019.
  * [33] D. Misra, “Mish: A self regularized non-monotonic activation function,” British Machine Vision Conference, 2020.
  * [34] A. Zela, T. Elsken, T. Saikia, Y. Marrakchi, T. Brox, and F. Hutter, “Understanding and robustifying differentiable architecture search,” in International Conference on Learning Representations, 2020.
  * [35] T.-Y. Lin, M. Maire, S. Belongie, J. Hays, P. Perona, D. Ramanan, P. Dollár, and C. L. Zitnick, “Microsoft coco: Common objects in context,” in European conference on computer vision, pp. 740–755, Springer, 2014.
  * [36] K. He, G. Gkioxari, P. Dollár, and R. Girshick, “Mask r-cnn,” in Proceedings of the IEEE international conference on computer vision, pp. 2961–2969, 2017.


| 28 | $0.22\cdot 10^{33}$ | $0.13\cdot 10^{8}$ | $0.62\cdot 10^{7}$ | 24 | $0.34\cdot 10^{-27}$ | 24 | $0.28\cdot 10^{-32}$
| 100 | $0.94\cdot 10^{33}$ | $0.29\cdot 10^{3}$ | $0.14\cdot 10^{1}$ | 24 | $0.53\cdot 10^{-31}$ | 24 | $0.37\cdot 10^{-32}$
34 | 100 | $0.70\cdot 10^{33}$ | $0.26\cdot 10^{3}$ | $0.15\cdot 10^{1}$ | 34 | $0.21\cdot 10^{-31}$ | 34 | $0.68\cdot 10^{-32}$
44 | 100 | $0.53\cdot 10^{33}$ | $0.23\cdot 10^{3}$ | $0.16\cdot 10^{1}$ | 44 | $0.54\cdot 10^{-32}$ | 44 | $0.40\cdot 10^{-32}$
54 | 100 | $0.42\cdot 10^{33}$ | $0.16\cdot 10^{3}$ | $0.18\cdot 10^{1}$ | 54 | $0.28\cdot 10^{-31}$ | 54 | $0.73\cdot 10^{-32}$
Table 2:  The results of computing the roots of the Wilkinson polynomial
$p_{\text{wilk}}(x)$ in extended precision, using our algorithm and eig, with
$\delta=10^{-3}$. †The error was so large here that some roots were outside of
the region (310).

######  Remark 5.1.

The remarkable stability of eig for many of the examples in Tables 1 and 2 is
explained by the following observation. The colleague matrix is the sum of a
tridiagonal matrix and a matrix that is all zeros except for the last row,
which is essentially equal to the coefficient vector $c$ (see formula (31)).
When there are large elements of $c$ near the tail of the vector, the
corresponding large entries in the colleague matrix cannot be balanced away,
since they are very close to the diagonal of the matrix. On the other hand,
when all of the elements of $c$ near the tail of the vector are relatively
small, and the large elements of $c$ appear near the head of the vector, these
large elements can be easily balanced away, since they are far from the
diagonal, and the corresponding elements on the other side of the diagonal are
all zero. The coefficient vector $c$ is usually large only because the last
coefficient of the corresponding non-monic Chebyshev expansion is small. If
the function being approximated by this non-monic Chebyshev expansion has been
adequately represented, then taking additional terms in the expansion will
result in corresponding expansion coefficients which are all machine epsilon
in size. Thus, adding terms to the Chebyshev expansion has the effect of
adding elements of size approximately one to the tail of the coefficient
vector $c$; if enough such elements are added, then all the large elements of
$c$ will be closer to the head of the vector, and can be balanced away. We
also observe that, not unexpectedly, the size of the largest eigenvalue of the
colleague matrix is approximately the same size as the norm of the colleague
matrix after balancing. Thus, if enough terms are taken in a Chebyshev
expansion, all of the eigenvalues of the colleague matrix will eventually be
small. All of this indicates that, provided enough terms are taken, a dense
eigensolver combined with balancing can result in a backward stable
rootfinding algorithm. Of course, we note that balancing the colleague matrix
destroys the Hermitian plus rank-1 structure, which bars the use of structured
QR algorithms depending on this property.

###  5.3 $f_{\text{sin}}(x)$: A Smooth Function

Here we construct an order-$n$ Chebyshev expansion of the smooth function

$\displaystyle f_{\text{sin}}(x)=\sin(2+20(x+0.222)^{2}).$ (312)

Since $f_{\text{sin}}(x)$ is analytic, its expansion coefficients decay
exponentially. When $i\geq 80$, the coefficients $a_{i}$ are around $10^{-14}$
in size (see Figure 5). If the coefficients are computed in extended
precision, then, when $i\geq 125$, the coefficients $a_{i}$ are around
$10^{-34}$ in size. Since the function is approximated accurately by a
Chebyshev expansion, its roots can be computed from the corresponding
colleague matrix (see, for example, [12] for a nice discussion). The results
of our numerical experiment are shown in Tables 3 and 4. Plots of the real and
complex roots of the order-$100$ Chebyshev expansion are shown in Figure 6.

Figure 5:  The function $f_{\text{sin}}(x)$, shown on the left, and the magnitude of its Chebyshev expansion coefficients $a_{i}$, shown on the right.  |  |  |  | eig |  Algorithm 4
---|---|---|---|---|---
$n$ | ${\lVert c\rVert}$ | ${\lVert\texttt{bal}(C)\rVert}$ | $\max_{i}{\lvert z_{i}\rvert}$ | $n_{\text{roots}}$ | $\max_{i}{\lvert\eta(p;\widehat{x}_{i})\rvert}$ | $n_{\text{roots}}$ | $\max_{i}{\lvert\eta(p;\widehat{x}_{i})\rvert}$
80 | $0.89\cdot 10^{15}$ | $0.16\cdot 10^{2}$ | $0.29\cdot 10^{1}$ | 14 | $0.74\cdot 10^{-14}$ | 14 | $0.10\cdot 10^{-13}$
100 | $0.14\cdot 10^{15}$ | $0.14\cdot 10^{2}$ | $0.11\cdot 10^{1}$ | 14 | $0.84\cdot 10^{-14}$ | 14 | $0.26\cdot 10^{-13}$
Table 3:  The results of computing the roots of the Chebyshev expansion of the function $f_{\text{sin}}(x)$, using our algorithm and eig, with $\delta=10^{-3}$.  |  |  |  | eig |  Algorithm 4
---|---|---|---|---|---
$n$ | ${\lVert c\rVert}$ | ${\lVert\texttt{bal}(C)\rVert}$ | $\max_{i}{\lvert z_{i}\rvert}$ | $n_{\text{roots}}$ | $\max_{i}{\lvert\eta(p;\widehat{x}_{i})\rvert}$ | $n_{\text{roots}}$ | $\max_{i}{\lvert\eta(p;\widehat{x}_{i})\rvert}$
125 | $0.92\cdot 10^{33}$ | $0.30\cdot 10^{2}$ | $0.25\cdot 10^{1}$ | 14 | $0.53\cdot 10^{-32}$ | 14 | $0.50\cdot 10^{-31}$
200 | $0.49\cdot 10^{32}$ | $0.29\cdot 10^{2}$ | $0.11\cdot 10^{1}$ | 14 | $0.12\cdot 10^{-31}$ | 14 | $0.60\cdot 10^{-31}$
Table 4:  The results of computing the roots of the Chebyshev expansion of the
function $f_{\text{sin}}(x)$ in extended precision, using our algorithm and
eig, with $\delta=10^{-3}$.  Figure 6:  The roots of the Chebyshev expansion
of order $100$ of the function $f_{\text{sin}}(x)$. The complex roots
$\widehat{z}_{i}$ are plotted with purple crosses ($+$) and the real roots
$\widehat{x}_{i}$ are plotted with blue stars ($\star$).

###  5.4 $p_{\text{mult}}(x)$: A Polynomial with Multiple Roots

Here we construct an order $n$ Chebyshev expansion of the degree $m$
polynomial

$\displaystyle
f_{\text{mult}}(x)=(x+\tfrac{1}{2})(x+\tfrac{1}{3})(x+0.61)(x-0.121)\prod_{i=1}^{m-4}(x-(1-10^{-3})).$
(313)

This polynomial has four simple roots on the interval $[-1,1]$, and a root of
multiplicity $(m-4)$ at the point $1-10^{-3}$ (see Figure 7). The results of
our numerical experiments are shown in Tables 5 and 6. We observe that, in
double precision, when the multiplicity of the root is greater than or equal
to 5, not all real roots are found. This is because the error in these roots
is approximately equal to $\epsilon^{\frac{1}{5}}$, and when $\epsilon\approx
10^{-14}$, we have that $\epsilon^{\frac{1}{5}}\approx 1.6\cdot 10^{-3}$; when
$\delta=10^{-3}$, this means that some of these roots can be outside the
region (310). Likewise, since when $\epsilon\approx 10^{-34}$,
$\epsilon^{\frac{1}{12}}\approx 1.4\cdot 10^{-3}$, it follows that in extended
precision, real roots are missed when their multiplicity is greater than or
equal to approximately 12. See the excellent discussion in [13] for more
details.

Figure 7:  The polynomial $p_{\text{mult}}(x)$ of order $7$ on the left, and order $9$ on the right. The roots are indicated with blue stars ($\star$).  |  |  |  |  | eig |  Algorithm 4
---|---|---|---|---|---|---
Degree | $n$ | ${\lVert c\rVert}$ | ${\lVert\texttt{bal}(C)\rVert}$ | $\max_{i}{\lvert z_{i}\rvert}$ | $n_{\text{roots}}$ | $\max_{i}{\lvert\eta(p;\widehat{x}_{i})\rvert}$ | $n_{\text{roots}}$ | $\max_{i}{\lvert\eta(p;\widehat{x}_{i})\rvert}$
7 | 100 | $0.10\cdot 10^{15}$ | $0.38\cdot 10^{2}$ | $0.11\cdot 10^{1}$ | 7 | $0.82\cdot 10^{-14}$ | 7 | $0.14\cdot 10^{-14}$
8 | 8 | $0.12\cdot 10^{3}$ | $0.50\cdot 10^{1}$ | $0.10\cdot 10^{1}$ | 8 | $0.59\cdot 10^{-15}$ | 8 | $0.93\cdot 10^{-15}$
| 9 | $0.54\cdot 10^{15}$ | $0.22\cdot 10^{13}$ | $0.22\cdot 10^{13}$ |  5† | $0.33\cdot 10^{-4}$ | 8 | $0.11\cdot 10^{-14}$
| 10 | $0.61\cdot 10^{15}$ | $0.16\cdot 10^{7}$ | $0.11\cdot 10^{7}$ |  6† | $0.14\cdot 10^{-9}$ | 8 | $0.88\cdot 10^{-15}$
| 11 | $0.79\cdot 10^{15}$ | $0.17\cdot 10^{5}$ | $0.93\cdot 10^{4}$ | 8 | $0.20\cdot 10^{-11}$ | 8 | $0.83\cdot 10^{-15}$
| 100 | $0.99\cdot 10^{14}$ | $0.32\cdot 10^{2}$ | $0.11\cdot 10^{1}$ | 8 | $0.64\cdot 10^{-14}$ | 8 | $0.26\cdot 10^{-15}$
9 | 100 | $0.97\cdot 10^{14}$ | $0.32\cdot 10^{2}$ | $0.11\cdot 10^{1}$ |  8⋄ | $0.99\cdot 10^{-14}$ |  8⋄ | $0.88\cdot 10^{-14}$
10 | 100 | $0.96\cdot 10^{14}$ | $0.26\cdot 10^{2}$ | $0.11\cdot 10^{1}$ |  8⋄ | $0.73\cdot 10^{-15}$ |  8⋄ | $0.38\cdot 10^{-15}$
13 | 100 | $0.92\cdot 10^{14}$ | $0.22\cdot 10^{2}$ | $0.11\cdot 10^{1}$ |  12⋄ | $0.12\cdot 10^{-14}$ |  12⋄ | $0.88\cdot 10^{-15}$
Table 5:  The results of computing the roots of the polynomial $p_{\text{mult}}(x)$, using our algorithm and eig, with $\delta=10^{-3}$. †The error was so large here that some roots were outside of the region (310). ⋄The multiplicity of the rightmost root was so large here that some roots were outside of the region (310).  |  |  |  |  | eig |  Algorithm 4
---|---|---|---|---|---|---
Degree | $n$ | ${\lVert c\rVert}$ | ${\lVert\texttt{bal}(C)\rVert}$ | $\max_{i}{\lvert z_{i}\rvert}$ | $n_{\text{roots}}$ | $\max_{i}{\lvert\eta(p;\widehat{x}_{i})\rvert}$ | $n_{\text{roots}}$ | $\max_{i}{\lvert\eta(p;\widehat{x}_{i})\rvert}$
10 | 100 | $0.69\cdot 10^{33}$ | $0.29\cdot 10^{3}$ | $0.13\cdot 10^{1}$ | 10 | $0.94\cdot 10^{-32}$ | 10 | $0.51\cdot 10^{-33}$
11 | 11 | $0.75\cdot 10^{4}$ | $0.98\cdot 10^{1}$ | $0.10\cdot 10^{1}$ | 11 | $0.76\cdot 10^{-33}$ | 11 | $0.12\cdot 10^{-32}$
| 12 | $0.16\cdot 10^{34}$ | $0.11\cdot 10^{30}$ | $0.11\cdot 10^{30}$ |  6† | $0.28\cdot 10^{-5}$ | 11 | $0.95\cdot 10^{-33}$
| 13 | $0.37\cdot 10^{34}$ | $0.53\cdot 10^{15}$ | $0.35\cdot 10^{15}$ |  7† | $0.63\cdot 10^{-19}$ | 11 | $0.66\cdot 10^{-33}$
| 14 | $0.25\cdot 10^{34}$ | $0.68\cdot 10^{10}$ | $0.35\cdot 10^{10}$ | 11 | $0.51\cdot 10^{-25}$ | 11 | $0.37\cdot 10^{-33}$
| 100 | $0.70\cdot 10^{33}$ | $0.37\cdot 10^{3}$ | $0.13\cdot 10^{1}$ | 11 | $0.11\cdot 10^{-31}$ | 11 | $0.15\cdot 10^{-32}$
12 | 100 | $0.71\cdot 10^{33}$ | $0.29\cdot 10^{3}$ | $0.13\cdot 10^{1}$ | 12 | $0.20\cdot 10^{-31}$ | 12 | $0.81\cdot 10^{-33}$
13 | 100 | $0.72\cdot 10^{33}$ | $0.30\cdot 10^{3}$ | $0.13\cdot 10^{1}$ | 13 | $0.17\cdot 10^{-31}$ | 13 | $0.11\cdot 10^{-32}$
14 | 100 | $0.72\cdot 10^{33}$ | $0.28\cdot 10^{3}$ | $0.13\cdot 10^{1}$ | 14 | $0.24\cdot 10^{-31}$ | 14 | $0.12\cdot 10^{-32}$
15 | 100 | $0.73\cdot 10^{33}$ | $0.28\cdot 10^{3}$ | $0.13\cdot 10^{1}$ |  9⋄ | $0.21\cdot 10^{-31}$ |  11⋄ | $0.21\cdot 10^{-31}$
Table 6:  The results of computing the roots of the polynomial
$p_{\text{mult}}(x)$ in extended precision, using our algorithm and eig, with
$\delta=10^{-3}$. †The error was so large here that some roots were outside of
the region (310). ⋄The multiplicity of the rightmost root was so large here
that some roots were outside of the region (310).

###  5.5 $p_{\text{yuji}}(x)$: A Pathological Example from [26]

Here we consider the order-$8$ polynomial

$\displaystyle p_{\text{yuji}}(x)=\sum_{i=0}^{8}a_{i}T_{i}(x),$ (314)

where the coefficient vector $a$ is given by

$\displaystyle
a=\left(\begin{array}[]{ccccccccc}-\frac{1}{10}&-\frac{1}{10}&-\frac{1}{10}&-\frac{1}{10}&-\frac{1}{10}&-\frac{1}{10}&10^{-10}&1&10^{-15}\end{array}\right),$
(316)

described in §6.1 of [26] (in [26], the authors set the last element of the
coefficient vector to $10^{-20}$; we set it close to machine epsilon instead).
Observe that the entry in the bottom right corner of the corresponding
colleague matrix is around $10^{15}$ in size. This polynomial has seven real
roots on the interval $[-1,1]$, and a single large imaginary root. We report
the results of our numerical experiment in Table 7. Clearly, eig struggles to
produce any accuracy at all, while our algorithm returns all the roots to
machine precision.

|  |  |  |  | eig |  Algorithm 4
---|---|---|---|---|---|---
Degree | $n$ | ${\lVert c\rVert}$ | ${\lVert\texttt{bal}(C)\rVert}$ | $\max_{i}{\lvert z_{i}\rvert}$ | $n_{\text{roots}}$ | $\max_{i}{\lvert\eta(p;\widehat{x}_{i})\rvert}$ | $n_{\text{roots}}$ | $\max_{i}{\lvert\eta(p;\widehat{x}_{i})\rvert}$
8 | 8 | $0.10\cdot 10^{16}$ | $0.50\cdot 10^{15}$ | $0.50\cdot 10^{15}$ | 7 | $0.21\cdot 10^{-1}$ | 7 | $0.77\cdot 10^{-14}$
Table 7:  The results of computing the roots of the polynomial
$p_{\text{yuji}}(x)$, using our algorithm and eig, with $\delta=10^{-3}$.

###  5.6 $f_{\text{cas}}(x)$: A Pathological Example from [14]

Here we consider the order-$n$ Chebyshev expansion of the smooth function

$\displaystyle f_{\text{cas}}(x)=\sin\Bigl{(}\frac{1}{x^{2}+10^{-2}}\Bigr{)},$
(317)

described in [14]. The first 1430 Chebyshev expansion coefficients of
$f_{\text{cas}}(x)$ are shown in Figure 8. This function is highly
oscillatory, and requires a Chebyshev expansion of order at least 1430 to
resolve it. Our numerical experiments are shown in Table 8. Plots of the real
and complex roots of the order-$1600$ Chebyshev expansion are shown in Figure
9.

Figure 8:  The magnitudes of the first 1430 Chebyshev expansion coefficients of $f_{\text{cas}}(x)$.  |  |  |  | eig |  Algorithm 4
---|---|---|---|---|---
$n$ | ${\lVert c\rVert}$ | ${\lVert\texttt{bal}(C)\rVert}$ | $\max_{i}{\lvert z_{i}\rvert}$ | $n_{\text{roots}}$ | $\max_{i}{\lvert\eta(p;\widehat{x}_{i})\rvert}$ | $n_{\text{roots}}$ | $\max_{i}{\lvert\eta(p;\widehat{x}_{i})\rvert}$
1430 | $0.16\cdot 10^{14}$ | $0.39\cdot 10^{2}$ | $0.10\cdot 10^{1}$ | 62 | $0.25\cdot 10^{-13}$ | 62 | $0.98\cdot 10^{-12}$
Table 8:  The results of computing the roots of the Chebyshev expansion of the
function $f_{\text{cas}}(x)$, using our algorithm and eig, with
$\delta=10^{-4}$.  Figure 9:  The roots of the Chebyshev expansion of order
$1600$ of the function $f_{\text{cas}}(x)$. The complex roots
$\widehat{z}_{i}$ are plotted with purple crosses ($+$) and the real roots
$\widehat{x}_{i}$ are plotted with blue stars ($\star$).

###  5.7 CPU Times

The CPU times of our algorithm are compared to the times of MATLAB’s eig in
Figure 10. These timing experiments were performed on polynomials with random,
independent, normally distributed Chebyshev expansion coefficients, with the
last coefficient chosen so that the vector $c$ has the desired norm (see
Section 5.1). We found that the CPU times do not depend on ${\lVert c\rVert}$,
so we report the results only for ${\lVert c\rVert}=2$. We observe that our
algorithm is strictly faster than eig, even for small inputs, except perhaps
for $n=7$, for which our algorithm and eig cost about the same. The growth in
CPU times taken by our algorithm agrees nicely with the expected asymptotic
cost of $O(n^{2})$, while eig shows a growth of $O(n^{3})$.

Figure 10:  The CPU times of our algorithm, plotted with purple crosses ($+$),
and the CPU times of eig, plotted with green x’s ($\times$), for various
values of $n$, where $n$ is the dimensionality of the colleague matrix.

##  6 Conclusions and Generalizations

In this manuscript, we describe an explicit, $O(n^{2})$ structured QR
algorithm for colleague matrices (more generally, for Hessenberg matrices that
are the sum of a Hermitian matrix and a rank-1 matrix), and prove that it is
componentwise backward stable. These results can be generalized in several
directions, of which we describe four. First, the algorithm can be modified in
a fairly straightforward way to work on Hessenberg matrices that are the sum
of a Hermitian matrix and a rank-$k$ perturbation (as opposed to a rank-1
perturbation). Like in the rank-1 case, most of the entries in the Hermitian
part are inferred from the low rank part, except that they are inferred from a
rank-$k$ matrix instead of a rank-1 matrix. The QR iteration proceeds
similarly, the main difference being that the correction in Line 12 of
Algorithm 1 becomes a correction to a row of an $n\times k$ matrix.

Second, the extension of this algorithm to an implicit, $O(n^{2})$ structured
QR algorithm that is also componentwise backward stable is fairly
straightforward. The key observation of this manuscript (that, to maintain
componentwise error bounds, a correction must be applied to the rank-1 part
whenever an entry of the matrix is eliminated) can be applied to a bulge-
chasing algorithm where the matrix is similarly represented by generators.

Third, we observe that this algorithm can be used to accelerate the
calculation of eigenvalues of general matrices, not necessarily in Hessenberg
form, that are representable as the sum of a Hermitian matrix and a rank-1 (or
rank-$k$) matrix. Such matrices can be quickly reduced to Hessenberg form in
$O(n^{3})$ operations, and once they are in Hessenberg form, our $O(n^{2})$
algorithm can be used to compute the eigenvalues. Thus, the cost of the
algorithm is dominated by the reduction to Hessenberg form, which will have a
much smaller constant than the standard algorithm for the evaluation of the
eigenvalues of the original dense matrix. Furthermore, if the reduction to
Hessenberg form can be done in a componentwise backward stable fashion, then
this scheme results in a componentwise backward stable eigensolver for general
matrices of the form Hermitian plus rank-1 (or rank-$k$).

Fourth, we observe that our algorithm can be used to find the roots of
polynomials expressed in other bases besides Chebyshev polynomials. It was
observed in [5] that, given any orthogonal polynomial basis that satisfies a
three-term recurrence relation, and given a polynomial expressed in that
basis, it is possible to construct an analogue of the colleague matrix from
the expansion coefficients. This matrix is a Hessenberg matrix that is the sum
of a (not necessarily symmetric) tridiagonal matrix and a rank-1 matrix;
matrices of this form are called _comrade matrices_. For all classical
orthogonal polynomials, the tridiagonal part can made symmetric by balancing,
without making any entries of the matrix much larger or much smaller. Our
algorithm can then be applied to this new matrix, which is a Hessenberg matrix
that is the sum of a symmetric tridiagonal matrix and a rank-1 matrix.

## References

  * [1] Advanpix Multiprecision Computing Toolbox for MATLAB. Yokohama, Japan: Advanpix LLC. See `http://www.advanpix.com/`.
  * [2] Aurentz, J., T. Mach, L. Robol, R. Vandebril, D.S. Watkins. Core-Chasing Algorithms for the Eigenvalue Problem. SIAM, 2018.
  * [3] Aurentz, J.L., T. Mach, L. Robol, R. Vandelbril, and D.S. Watkins. “Fast and backward stable computation of roots of polynomials, part II: backward error analysis; companion matrix and companion pencil.” SIAM J. Matrix Anal. Appl., 39.3 (2018): 1245–1269.
  * [4] Aurentz, J.L., T. Mach, R. Vandelbril, and D.S. Watkins. “Fast and backward stable computation of roots of polynomials.” SIAM J. Matrix Anal. Appl., 36.3 (2015): 942–973.
  * [5] Barnett, S. “A companion matrix analogue for orthogonal polynomials.” Linear Algebra Appl. 12 (1975): 197–208.
  * [6] Battles, Z. and L.N. Trefethen. “An extension of Matlab to continuous functions and operators.” SIAM J. Sci. Comput. 25 (2004): 1743–1770.
  * [7] Bauer, F.L., C.T. Fike. “Norm and exclusion theorems.” Numer. Math. 2 (1960): 137–141.
  * [8] Van Barel, M., R. Vandebril, P. Van Dooren, K. Frederix. “Implicit double shift QR-algorithm for companion matrices.” Numer. Math. 116 (2010): 177–212.
  * [9] Bini, D.A., P. Bonito, Y. Eidelman, L. Gemignani, and I. Gohberg. “A fast implicit QR eigenvalue algorithm for companion matrices.” Linear Algebra Appl. 432 (2010): 2006–2031.
  * [10] Bini, D.A., Y. Eidelman, L. Gemignani, and I. Gohberg. “Fast QR eigenvalue algorithms for Hessenberg matrices which are rank-one perturbations of unitary matrices.” SIAM J. Matrix Anal. Appl. 29.2 (2007): 566–585.
  * [11] Bini, D.A., L. Gemignani, V.Y. Pan. “Fast and stable QR eigenvalue algorithms for generalized companion matrices and secular equations.” Numer. Math. 100 (2005): 373–408.
  * [12] Boyd, J.P. “Finding the Zeros of a Univariate Equation: Proxy Rootfinders, Chebyshev Interpolation, and the Companion Matrix.” SIAM Rev. 55.2 (2013): 375–396.
  * [13] Boyd, J.P and D.H. Gally. “Numerical experiments on the accuracy of the Chebyshev-Frobenius companion matrix method for finding the zeros of a truncated series of Chebyshev polynomials.” J. Comput. Appl. Math. 205 (2007): 281–295.
  * [14] Casulli, A. and L. Robol. “Rank-structured QR for Chebyshev rootfinding.” arXiv:2010.11416v1 [math.NA], Oct. 2020.
  * [15] Chandrasekaran, S., M. Gu, J. Xia, J. Zhu. “A fast QR algorithm for companion matrices.” Recent Advances in Matrix and Operator Theory, Oper. Theory Adv. Appl. 179 (2008): 111–143.
  * [16] Chebfun for MATLAB. See `http://www.chebfun.org/`.
  * [17] Corless, Robert M. “Generalized companion matrices in the Lagrange basis.” Proceedings EACA (pp. 317–322). Santander, Spain: Universidad de Cantabria, 2004.
  * [18] Corless, R.M. Personal Communication. 13 July 2020.
  * [19] Edelman, A. and H. Murakami. “Polynomial roots from companion matrix eigenvalues.” Math. Comput. 64.210 (1995): 763–776.
  * [20] Eidelman, Y., L. Gemignani, and I. Gohberg. “Efficient eigenvalue computation for quasiseparable Hermitian matrices under low rank perturbations.” Numer. Algor. 47 (2008): 253–273.
  * [21] Good, I.J. “The colleague matrix, a Chebyshev analogue of the companion matrix.” Q. J. Math., 2.12 (1961): 61–68.
  * [22] Higham, N.J. Accuracy and Stability of Numerical Algorithms. 2nd ed. Philadelpha: SIAM, 2002.
  * [23] Lawrence, P.W., M. Van Barel, P. Van Dooren. “Backward error analysis of polynomial eigenvalue problems solved by linearization.” SIAM J. Matrix Anal. Appl. 37.1 (2016): 123–144.
  * [24] McNamee, J.M. Numerical Methods for Roots of Polynomials, Part I. Elsevier, 2007.
  * [25] McNamee, J.M. and V.Y. Pan. Numerical Methods for Roots of Polynomials, Part II. Elsevier, 2013.
  * [26] Nakatsukasa, Y. and V. Noferini. “On the stability of computing polynomial roots via confederate linearizations.” Math. of Comput., 85.301 (2016): 2391–2425.
  * [27] Noferini, V., L. Robol, and R. Vandebril. “Structured backward errors in linearizations.” arXiv:1912.04157v1 [math.NA], Dec. 2019.
  * [28] Pan, V.Y. “Solving a polynomial equation: some history and recent progress.” SIAM Rev. 39.2 (1997): 187–220.
  * [29] Parlett, B.N. and C. Reinsch. “Balancing a matrix for calculation of eigenvalues and eigenvectors.” Numer. Math. 13 (1963), 293–304.
  * [30] Pérez, J. and V. Noferini. “Chebyshev rootfinding via computing eigenvalues of colleague matrices: what is it stable?” Math. Comput. 86.306 (2017): 1741–1767.
  * [31] Peters, G. and J.H. Wilkinson. “Practical problems arising in the solution of polynomial equations.” J. Inst. Maths. Appl. 8 (1971), 16–35.
  * [32] Serkh, K. “A fast, simple, and remarkably stable QR method for colleague matrices, and connections to event detection for the numerical solution of ODEs.” 2020 SIAM/CAIMS: 2nd Joint Annual Meeting, July, 2020. See `https://meetings.siam.org/sess/dsp_talk.cfm?p=106559`.
  * [33] Specht, W., “Die Lage der Nullstellen eines Polynoms III.” Math. Nach. 16 (1957): 369–389.
  * [34] Specht, W. “Die Lage der Nullstellen eines Polynoms IV.” Math. Nach. 21 (1960): 201–222.
  * [35] Tisseur, F. “Backward stability of the QR algorithm.” Equipe d’Analyse Numerique, Université Jean Monnet de Saint-Etienne Technical Report 239 (1996).
  * [36] Trefethen, L.N. Approximation Theory and Practice. SIAM, 2013.
  * [37] Wilkinson, J.H. The Algebraic Eigenvalue Problem. Oxford: Oxford University Press, 1965.