Dissertation submitted in fulfillment of the requirements for

the degree of Doctor of Philosophy

under the joint supervision of

Benjamin Negrevergne, Yann Chevaleyre and Jamal Atif

Pour commencer, je souhaiterais remercier Teddy Furon et Alain Rakotomamonjy d'avoir accepté d'être examinateurs de cette thèse ainsi qu'Élisa Fromont, Rémi Gribonval et Krzysztof Choromanski de s'être intéressés à mes travaux de recherche et d'avoir accepté d'être membres du jury.
Nos échanges pendant la relecture ainsi que la soutenance ont été très enrichissants.

J'éprouve une profonde gratitude envers mes trois encadrants de thèse, Jamal, Yann et Benjamin durant ces quatre dernières années où j'ai pu découvrir le monde de la recherche et me former à devenir un bon chercheur.
En effet, le métier de chercheur est bien le plus beau métier du monde !
Merci d'avoir cru en moi et de m'avoir offert cette chance, merci pour tous ces échanges, désaccords, réunions, sessions de travail et enfin, merci de m'avoir toléré quand j'étais pénible.
Vous m'avez vraiment offert un encadrement exceptionnel: peu de doctorants peuvent se targuer de pouvoir échanger avec leurs encadrants tous les jours !

Cette expérience a été d'autant plus riche grâce aux autres doctorants du laboratoire.
Ainsi, je souhaite remercier Rafael et Laurent pour les collaborations réalisées ensemble, et également Geovani, Raphaël, Alexandre V., Éric, Virginie et Céline pour avoir participé à l'émulation et la bonne ambiance du labo.
Merci également à Florian, Clément et Alexandre A. pour les échanges que l'on a pu avoir.

Cette thèse n'aurait pas été possible sans les financements de Wavestone.
Ainsi, je souhaite remercier Cyril pour m'avoir donné cette opportunité ainsi que Nicolas pour nos longs échanges sur mes travaux de recherche.
Également, je souhaite remercier David, Hugo et Cédric pour m'avoir aidé au cours de cette expérience.

Je souhaite également remercier ma famille, belle-famille et mes amis pour m'avoir toléré, accompagné, conseillé, encouragé et pour avoir essayé de comprendre ce que je faisais ces quatre dernières années.
Désolé d'avoir autant travaillé pendant les week-ends, j'espère ne pas avoir été trop désagréable ou trop “dans mon monde” pendant cette période.
Un grand merci à Othmane qui, à travers nos échanges, m'a donné l'idée de faire une thèse et de me lancer dans la recherche.

Pour finir, je souhaite exprimer ma profonde gratitude envers Chuthima, ma partenaire du quotidien.
Merci pour tous les compromis faits pour me laisser réaliser cette thèse, pour le soutien quotidien, les encouragements et les relectures.
La qualité de cette thèse est énormément due à ta présence à mes côtés.

Ecrire Acknowledgments en français

Deep neural networks are state-of-the-art in a wide variety of tasks, however, they exhibit important limitations which hinder their use and deployment in real-world applications.
When developing and training neural networks, the accuracy should not be the only concern, neural networks must also be cost-effective and reliable.
Although accurate, large neural networks often lack these properties.

In this thesis, we leverage the properties of structured matrices from the Toeplitz family to build compact and secure neural networks.
Our contributions are twofold.

First, we propose a new neural network architecture that is not only accurate but also compact and easy to train.
The purpose of this contribution is to study deep diagonal-circulant neural networks, which are deep neural networks in which weight matrices are the product of diagonal and circulant ones.

We perform a theoretical analysis of their expressivity and propose an initialization procedure and an intelligent use of nonlinearity functions to facilitate training.
Furthermore, we show that these networks outperform recently introduced deep networks with other types of structured layers.
We conduct a thorough experimental study to compare the performance of deep diagonal-circulant networks with state-of-the-art models based on structured matrices and with dense models.
We show that our models achieve better accuracy than other structured approaches while requiring 2x fewer weights than the next best approach.
Finally, we train compact and accurate deep diagonal-circulant networks on a real-world video classification dataset with over 3.8 million training examples.

Secondly, we propose an approach to build robust neural networks to adversarial examples.
In this contribution, we introduce a new Lipschitz regularization for Convolutional Neural Networks that improves the robustness of neural networks.

Lipschitz regularity is now established as a key property of modern deep learning with implications in training stability, generalization, robustness against adversarial examples, etc.
However, computing the exact value of the Lipschitz constant of a neural network is known to be NP-hard.
Recent attempts from the literature introduce upper bounds to approximate this constant that are either efficient but loose or accurate but computationally expensive.
In this work, by leveraging the properties of doubly-block Toeplitz matrices, we introduce a new upper bound of the singular values of convolution layers that is both tight and easy to compute.
Based on this result we devise an algorithm to train Lipschitz-regularized Convolutional Neural Networks.

Les réseaux de neurones profonds sont considérés comme étant état de l'art dans une grande variété de tâches, mais ils présentent des limites importantes qui entravent leur utilisation et leur déploiement.
Lors du développement et l'entraînement de réseaux de neurones, la précision ne devrait pas être la seule préoccupation, ils se doivent aussi d'être efficaces et sécurisés.
Bien que précis, les réseaux de neurones dotés de nombreux paramètres n'ont souvent pas ces propriétés.
Dans cette thèse, nous exploitons les propriétés des matrices structurées de la famille de Toeplitz pour construire des réseaux de neurones compacts et sécurisés.
Nous réalisons deux contributions sur ces thématiques.

Premièrement, nous proposons une nouvelle architecture de réseau de neurones précise, mais également compacte et facile à entraîner.
L'objectif de cette contribution est d'étudier les réseaux de neurones diagonaux-circulants, qui sont des réseaux de neurones profonds pour lesquels les matrices de poids sont le produit des matrices diagonales et circulantes.
Nous effectuons une analyse théorique de leur expressivité et proposons une procédure d'initialisation et une utilisation intelligente des fonctions de non-linéarité qui facilitent leur entraînement.

Nous montrons que nos modèles atteignent une meilleure précision que les autres approches structurées tout en nécessitant deux fois moins de paramètres.
Enfin, nous entraînons des réseaux de neurones diagonaux-circulants sur un ensemble de données de classification vidéo qui contient plus de 3,8 millions d'exemples.

Deuxièmement, en plus d'être compacts et précis, les réseaux de neurones se doivent d'être sécurisés.
Pour améliorer leur robustesse, nous proposons une nouvelle régularisation pour les réseaux convolutifs basée sur la constante de Lipschitz.
La régularisation Lipschitz est maintenant établie comme une propriété clé de l'apprentissage profond avec des implications en stabilité, généralisation et robustesse contre les attaques adversariales, etc.
Cependant, le calcul de la constante de Lipschitz d'un réseau de neurones est connu pour être un problème NP-complet.
De récentes tentatives introduisent des bornes supérieures pour approximer cette constante qui sont soit efficaces, mais peu précises, soit précises mais coûteuses.
Dans cette thèse, en exploitant les propriétés des matrices de Toeplitz à bloc de Toeplitz, nous introduisons une nouvelle borne supérieure de cette constante pour les couches convolutionnelles qui est à la fois précise et facile à calculer.
Sur la base de ce résultat, nous concevons un algorithme pour entraîner des réseaux de neurones convolutifs avec une régularisation Lipschitz.

tocchapterTable of Contents
*Table of Contents
[title=List of Symbols]

CHAPTER: INTRODUCTION

§ CONTEXT AND MOTIVATION

The neural network architecture (AlexNet) proposed by <cit.> which won the ImageNet Large-Scale Visual Recognition Challenge in 2012.

Since the dawn of computer science, researchers have been trying to emulate intelligence through computers.
Alan Turing was the first, in a paper called Computing Machinery and Intelligence <cit.>, to lay the foundation for what we now call Artificial Intelligence.
In the last 20 years, with the surge in data collection and computing resources, the interest and use cases for Machine Learning have grown exponentially.
More specifically, Deep Learning, a subfield of Machine Learning, consisting of training Deep Neural Networks on high-level data (images, sounds, texts) have shown great achievements, even outperforming humans on certain tasks.

One of the most remarkable breakthroughs of Deep Learning happened in 2012 during the ImageNet Large-Scale Visual Recognition Challenge <cit.>.
The challenge aims at evaluating different algorithms for object detection and image classification.
In 2012,  obtained 1 place and beat every other participant by a 10.8% margin with a neural network architecture called AlexNet.
The main reasons for this success are twofold.
First, they used a convolutional neural network (CNN) with more than 60 million parameters which was one of the largest models of the time.
Secondly, they designed a specific architecture to exploit dual programmable graphics processing units (GPUs) to speed up the arithmetic operations, which enabled them to significantly reduce training time.
<Ref> shows the AlexNet architecture which consists of five convolution layers with two fully connected layers at the end.

table-number-alignment = center,
table-space-text-pre =    ,

Authors     Models     2c#Params     TOP-5 Acc.

<cit.>     AlexNet                  61     M     84.7%

<cit.>           VGG                     144     M     92.0%

<cit.>                 ResNet-152               60     M     93.8%

<cit.>       Inception-ResNet-v2      56     M     95.1%

<cit.>          ResNeXt-101              84     M     95.6%

<cit.>              SENet                   146     M     96.2%

<cit.>        AmoebaNet-A             469     M     96.7%

<cit.>             AmoebaNet-B             556     M     97.0%

Computer Vision Models

Authors     Models     2c#Params

<cit.>             ELMo                 94     M

<cit.>       GPT                 110     M

<cit.>             BERT                340     M

<cit.>              XLNet (Large)       340     M

<cit.>             RoBERTa (Large)     355     M

<cit.>        GPT-2                 1     B

<cit.>        MegatronLM            8     B

<cit.>        T5-11B               11     B

<cit.>        T-NLG                17     B

<cit.>          GPT-3               175     B

<cit.>            Switch Transformers     1     T

Natural Language Processing Models

Evolution of the number of parameters for Computer Vision and Natural Language Processing models developed in the years after AlexNet.

Following this result, many architectures with an increasing number of parameters have been developed.
This growth in the number of parameters has led to an substantial gains in accuracy, exceeding even human performance, on the ImageNet dataset <cit.>.
<Ref> shows a list of the different state-of-the-art architectures along with their size and accuracy.
As we can see, the accuracy of the models generally improves at the cost of the model size.
For computer vision models, <cit.> have empirically shown that the relationship between model size and accuracy seems to obey a power law.

This relationship has also been observed for neural networks designed for Natural Language Processing (NLP) <cit.> aided by the availability of large-scale datasets such as the Common Crawl dataset <cit.> which constitutes nearly a trillion words.

As a result of their size and improved accuracy, deep neural networks now achieve state-of-the-art performances in a variety of domains such as image recognition <cit.>, object detection <cit.>, natural language processing <cit.>, speech recognition <cit.>, games <cit.>, etc.
Specifically, computer vision and natural language processing models have achieved sufficient performance for being used in real-world applications such as autonomous vehicles <cit.>, translation <cit.>, vocal assistants <cit.>, etc.

However, accuracy is not the only concern, when implemented in a critical decision process, neural networks need to be compact, cost-effective and secure.
Although accurate, large neural networks often lack these properties.
Indeed, training state-of-the-art models on computer vision or natural language processing tasks requires gigabytes of memory and can take several months on a single GPU <cit.>.
For example, the GPT-3 model proposed by <cit.>, culminates at 175 billion parameters and requires 355 years of training on a single GPU and $4600000 to train on a cloud-computing platform <cit.>.
It has also been estimated by <cit.> that the training and development costs of the large Transformer model proposed by <cit.> with neural architecture search emits an estimated 284019 kg of $\mathrm{CO}_2$ whereas a human life will consume an average of 5000 kg of $\mathrm{CO}_2$ for one year.
Furthermore, with the rise of smartphones and “Internet of things” devices (IoT) with limited computational and memory resources, neural networks also need to be efficient during the inference phase.
In addition, with the growing concern over data privacy, methods such as federated learning are gaining ground.
Federated learning involves training a model across multiple decentralized devices (, smartphones) with local data samples.
This avoids the step of centralizing all users' data into one server, thus addressing, even modestly, the issue of data privacy.
Thus, building compact and cost-effective neural networks have been an important goal in order to reduce training time, reduce cost and allow for faster research and development.

In addition to being compact and cost-effective, neural networks also need to be secure.
Due to their high complexity and expressivity, large neural networks exhibit instability to small perturbations of their inputs.
Unstable neural networks tend to be vulnerable to adversarial examples, , imperceptible variations of natural examples, crafted to deliberately mislead the models <cit.>.
<Ref> gives an example of an adversarial attack on an image.
The small perturbation (center) is added to the original image (left) leading to an adversarial image (right).
This behavior can cause serious security problems when neural networks are used for critical decision-making (, self-driving cars, predictive justice, etc.).

This thesis focuses on the problem of training neural networks which are not only accurate but also compact, easy to train, reliable and robust to adversarial examples.

§ PROBLEM STATEMENT AND CONTRIBUTIONS

Neural networks, which find their roots in the work of <cit.>, can be analytically described as a composition of multi-dimensional linear functions interlaced with nonlinear functions (also called activation functions).
More formally, a neural network is a function $N_{\Omega} : \Rbb^n \rightarrow \Rbb^m$ parameterized by a set of weights $\Omega$ of the form:

\begin{equation} \label{equation:ch1-neural_network}
N_{\Omega}(\xvec) = \psi^{(\depth)} \circ \rho \circ \psi^{(\depth-1)} \cdots \circ \psi^{(2)} \circ \rho \circ \psi^{(1)} (\xvec) \enspace.
\end{equation}

Here, $\depth$ corresponds to the depth of the network (, the number of layers) and $\rho$ is a nonlinear function.
Each $\psi^{(i)}$ is a multi-dimensional linear function $\psi^{(i)}: \xvec \mapsto \Wmat^{(i)} \xvec + \bvec^{(i)}$ parameterized by a weight matrix $\Wmat^{(i)}$ and a bias vector $\bvec^{(i)}$ and $\Omega$ is the union of the parameters of all the layers.

Classical neural networks typically have a large number of parameters to train.

If they have no restrictions on the weight matrices $\Wmat^{(i)}$, the layers are said to be fully connected.
Typically, fully connected neural networks have a large number of parameters.
For example, a fully connected neural network with $\depth$ layers and $n$ neurons on each layer ($\Wmat^{(i)} \in \Rbb^{n \times n}$) will have $pn (n + 1)$ parameters.
Since the input and output dimensions are generally large (, ImageNet has an input dimension of $224^2 \times 3$ and an output of 1000), simple fully connected neural networks with few layers accumulate over hundreds of millions of parameters.
Generally, this type of neural network has been shown to perform poorly due to a large search space or due to the important expressivity of the model which leads to overfitting[For Machine Learning models, overfitting is a well understood phenomenon. However, it has been discovered that large deep neural networks exhibit a “double descent” phenomenon (see <cit.>), where the performance first gets worse (overfits) then gets  better with longer training.].
Moreover, they are computationally expensive, which makes them impractical for a number of use cases (smartphones, IoT devices, etc.).
To reduce the number of parameters on each layer, researchers have devised specific linear operations that reduce the number of parameters and have better properties for the problem at hand.

Example of Adversarial Attack on an image.

An example of widely used neural networks with specialized and more compact linear operations are Convolutional Neural Networks (CNN) <cit.> which achieve state-of-the-art results for computer vision tasks.
Convolutional neural networks, which find their roots in the work of <cit.>, use specific weight matrices which encode the translation invariant property often desirable to process images.

Whereas a classical linear layer with a dense matrix will have $n \times n$ parameters, a convolution layer only has $k \times k$ parameters where $k \ll n$ is the kernel size and is usually small (, 3 or 5 for classical convolution layers).
A convolutional neural network is the most common type of structured neural networks.
Indeed, the convolution operation can be represented by a structured matrix , a matrix that can be represented with less than $n^2$ parameters.

In addition to offering a more compact representation, the structure of certain matrices can be exploited to obtain better algorithms for the matrix-vector product, thus optimizing memory and computing operations.
Based on the success of convolutional neural networks, researchers have studied and proposed other types of neural networks based on weight matrices with different structures (,  <cit.>).
<Ref> shows different types of structured matrices that have been used for deep learning.
Although convolutional neural networks have been state-of-the-art for computer vision tasks, it remains unclear whether other types of structured networks can be beneficial to other types of applications and which type of structure can provide both accuracy and efficient computation.

The contributions of this thesis lie at the intersection of linear algebra, Fourier analysis and deep learning.
As a result, we build compact and secure neural networks by leveraging the properties of structured matrices from the Toeplitz family.
Hereafter, we summarize our contributions.

§.§ Training Compact Neural Networks

As a first contribution, we use circulant matrices, which are a particular type of matrix from the Toeplitz family, to devise a new compact architecture replacing fully connected neural networks.
More precisely, we study deep diagonal-circulant neural networks, which are deep neural networks in which weight matrices are replaced by the product of diagonal and circulant ones.
Besides making a theoretical analysis of their expressivity, we introduce principled techniques for training these models: we devise an initialization scheme and propose a smart use of nonlinearity functions in order to train deep diagonal-circulant networks.
Furthermore, we show that these networks outperform recently introduced deep networks with other types of structured layers.
We conduct a thorough experimental study to compare the performance of these networks with state-of-the-art models.

We show that our models achieve better accuracy than other structured approaches while requiring 2x fewer weights than the next best approach.
Finally, we train accurate deep diagonal-circulant networks on a real-world video classification dataset with over 3.8 million training examples.

This contribution has been the subject of the following publications:

* Training Compact Deep Learning Models for Video Classification using Circulant Matrices
in the European Conference on Computer Vision Workshops on Video Classification

* Understanding and Training Deep Diagonal Circulant Neural Networks in the
24th European Conference on Artificial Intelligence.

§.§ Training Robust Neural Networks

As a second contribution, we build robust neural networks by studying the properties of the structure of convolutions.
We devise a new upper bound on the largest singular value of convolution layers that is both tight and easy to compute.
Our work is based on the result of <cit.> which states that an upper bound on the singular value of Toeplitz matrices can be computed from the inverse Fourier transform of the characteristic sequence of these matrices.
From our analysis immediately follows an algorithm for bounding the Lipschitz constant of a convolution layer, and by extension the Lipschitz constant of the whole network.
Finally, we illustrate our approach to adversarial robustness.
Recent work has shown that empirical methods such as adversarial training offer poor generalization <cit.> and can be improved by applying Lipschitz regularization <cit.>.
To illustrate the benefit of our new method, we train neural networks with Lipschitz regularization and show that it offers a significant improvement over adversarial training alone.

Additional joint contributions have also been made on the topic of robust neural networks.
A first work studied the effectiveness of noise injection at training and inference time in neural networks to protect against adversarial attacks.
In this work, we have shown that noise drawn from the Exponential family offers a provable protection against adversarial attacks.
A follow-up work conducts a geometrical analysis of defense mechanisms designed to protect neural networks against several types of attacks.
This work shows that neural networks designed to be robust against one type of adversarial example offer poor protection against other types of attacks.

The contribution on adversarial robustness has been the subject of the following publications:

* On Lipschitz Regularization of Convolutional Layers using Toeplitz Matrix Theory
in the 35 AAAI Conference on Artificial Intelligence

* Theoretical evidence for adversarial robustness through randomization in the
Advances in Neural Information Processing Systems.

* Advocating for Multiple Defense Strategies against Adversarial Examples in the
European Conference on Machine Learning Workshop for CyberSecurity

§ OUTLINE OF THE THESIS

This thesis is organized in six chapters.
First, <Ref> gives an introduction to the theory of Toeplitz matrices and on supervised learning and neural networks.
This chapter presents the necessary technical tools we will need for presenting the related work and for our contributions.
<Ref> is dedicated to discussing the state-of-the-art approaches related to our contributions.
The chapter is divided into two parts.
First, we review some techniques to build compact neural networks with an important focus on techniques that use structured matrices.
The second part focuses on presenting regularization methods for improving the robustness of neural networks.
<Ref> and <Ref> constitute our main contributions.
<Ref> presents results on compact neural networks built from diagonal and circulant matrices.
<Ref> presents our new regularization scheme to improve the robustness of neural networks based on the properties of doubly-block Toeplitz matrices.
<Ref> proposes a discussion and some perspectives on our contributions.
Appendix <ref> constitutes some complements to <Ref>.
It provides additional experiments on video classification with compact neural networks.
Finally, Appendices <ref> and  <ref> provide further work on the robustness of neural networks done during this Ph.D. thesis.

CHAPTER: BACKGROUND

This chapter gives an overview on the theory of Toeplitz matrices and on supervised learning with neural networks.
The first section describes the mathematical properties of Toeplitz matrices and some known theorems that we use in this thesis.
A Toeplitz matrix, named after Otto Toeplitz, is a matrix in which each descending diagonal, from left to right, is constant.
This simple property has led to interesting theoretical results and numerous applications.
We will use a number of these results in the context of neural networks.

The second section of this chapter is divided into four parts.
First, we review notions of supervised learning which refer to the problem of optimizing the parameters of a function in order to map an input to an output based on a series of input-output pairs.
Then, we formally define neural networks and recall some of their properties.
We pursue by introducing the concept of adversarial examples which we will use in <Ref>.
Finally, we present some recent theoretical results on neural networks that allow a better understanding of the contributions of this thesis.

§ A PRIMER ON CIRCULANT AND TOEPLITZ MATRICES

§ SUPERVISED LEARNING AND NEURAL NETWORKS

§ SUMMARY OF THE CHAPTER

As explained in the Introduction (<Ref>), our contributions lie at the intersection between neural networks and structured matrices.
In this chapter, we have reviewed the necessary concepts to present our contributions and some related works.

First, <Ref> introduced circulant and Toeplitz matrices which are the main mathematical objects used in this thesis.
Circulant and Toeplitz matrices are structured matrices in which each descending diagonal, from left to right, is constant.
These structured matrices are the building blocks of our contribution on compact neural networks (<Ref>) and enable fast approximation of the Lipschitz constant of convolution layers leading to a new regularization scheme (<Ref>).

Finally, in <Ref>, we gave a quick overview of the concept of supervised learning, which presents the mathematical tools for optimizing a parameterized function in order to map an input to an output based on a series of input-output pairs.
Although the statistical learning framework considers generic hypothesis space, in this work we use a class of functions called neural networks presented in <Ref>.
We also presented, in <Ref>, the concept of adversarial attacks and robustness of neural networks.
We showed how a neural network can be sensitive to small perturbations to its input and thus vulnerable to adversarial examples.
Reducing the sensitivity and therefore increasing the robustness of neural networks is the central theme of our second contribution presented in <Ref>.
Finally, in <Ref>, we have presented some recent results on the theory of neural networks.
These results give us insights on how neural networks generalize and a theoretical justification of the regularization scheme that we propose in <Ref>.

CHAPTER: RELATED WORK

This chapter, divided into two parts, is intended to provide an overview of the state of the art related to our contributions.
First, we present current methods for building compact neural networks.
Since the scope of application of these techniques is broad, we have chosen to focus mainly on work that uses linear algebra tools and more particularly structured matrices.
We present in the first subsection an overview of general techniques for building compact neural networks.
In the next subsection, we present in more detail the current methods for building compact neural networks with structured matrices.
Finally, we discuss these techniques with respect to our contribution to compact neural networks.

The second part of this chapter presents current methods for regularizing the Lipschitz constant of neural networks with the aim of improving their robustness.
This section is divided into four parts.
First, we present techniques that focus on the computation of the Lipschitz constant of neural networks.
Although theoretically and empirically interesting, we will see how these techniques do not scale and therefore cannot be applied to current neural network architectures.
The following subsection presents the approach of Lipschitz regularization via the Lipschitz constant of individual layers of the networks.
We describe the advantages and disadvantages of this approach.
Moreover, in the third subsection, we focus our presentation on current techniques that compute the singular values of convolution layers.
Finally, we discuss these techniques with respect to our contribution on Lipschitz regularization of convolutional neural networks.

§ RELATED WORK ON COMPACT NEURAL NETWORKS

§.§ Discussion

In this section, we have shown current methods and techniques for designing compact neural networks with structured matrices.
Our contributions on Deep Diagonal Circulant Neural Networks are a direct follow-up to the work of <cit.> focusing on compact neural networks with structured matrices.
More precisely, we extend the work of <cit.> by training fully structured networks (, networks with structured layers only) hence demonstrating that diagonal-circulant layers are able to model complex relations between inputs and outputs.
Although, this diagonal-circulant layers fit in the low displacement rank framework, we demonstrate much better performances in practice.
Indeed, thanks to a solid theoretical analysis and thorough experiments, we were able to train deep (up to 40 layers) circulant neural networks, and apply, for the first time, this structured architecture in the context of large-scale video classification.
This contrasts with previous experiments in which only one or a few dense layers were replaced inside a large redundant network such as VGG <cit.>.

§ RELATED WORK ON LIPSCHITZ REGULARIZATION

§.§ Discussion

We have presented state-of-the-art methods for regularizing the Lipschitz constant of neural networks with the aim to improve their robustness against adversarial attacks.
The power method <cit.> is a popular technique for approximating the maximal singular value of a matrix.

Recent works in deep learning use this method in a wide variety of settings, for example, robustness <cit.>, generalization <cit.> or to stabilize the training of Generative Adversarial Networks (GANs) <cit.>.
Despite a number of interesting results, using the power method is expensive and results in prohibitive training times.
Other approaches to regularize the Lipschitz constant of neural networks have been proposed by <cit.> and  <cit.>.
These methods exploit the properties of circulant matrices to approximate the maximal singular value of a convolution layer.
Although interesting, theses method results in a loose approximation of the maximal singular value.
Our work is positioned at the intersection between these works, we will introduce a new approach for regularizing the Lipschitz constant of neural networks, that is more efficient than the power method and more accurate than methods relying on the structure of convolutions.

CHAPTER: DIAGONAL AND CIRCULANT MATRICES FOR COMPACT NEURAL NETWORKS

§ INTRODUCTION

As seen in the previous chapters, structured matrices are at the very core of most of the work on compact networks.

Despite substantial efforts  (, <cit.>), the performance of compact models is still far from achieving an acceptable accuracy motivating their use in real-world scenarios.
This raises several questions about the effectiveness of such models and about our ability to train them.
In particular two main questions call for investigation:

* What is the expressive power of structured layers compared to dense layers?

* How to efficiently train deep neural networks with a large number of structured layers?

We aim at answering these questions by studying deep diagonal-circulant neural networks (DCNNs), which are deep neural networks in which weight matrices are the product of diagonal and circulant ones.

To answer the first question, we propose an analysis of the expressivity of DCNNs by extending the results obtained by <cit.> which states that any matrix can be decomposed into the product of $2n-1$ alternating diagonal and circulant matrices.
We introduce a new bound on the number of diagonal-circulant products required to approximate a matrix that depends on its rank.
Building on this result, we demonstrate that a DCNN with bounded width and small depth can approximate any dense neural networks with ReLU activations.

To answer the second question, we first describe a theoretically sound initialization procedure for DCNN which allows the signal to propagate through the network without vanishing or exploding.
Furthermore, we provide a number of empirical insights to explain the behavior of DCNNs and show the impact on the number of nonlinearities in the network on the convergence rate and the accuracy of the network.
By combining all these insights, we are able (for the first time) to train large and deep DCNNs and demonstrate the good performance of these networks on a large-scale application (the video classification problem) and obtain very competitive accuracy.

The chapter is organized as follows:
<Ref> introduces our new result extending the one from <cit.>.
<Ref> proposes a theoretical analysis on the expressivity of DCNNs.
<Ref> describes two efficient techniques for training deep diagonal-circulant neural networks.
<Ref> presents extensive experiments to compare the performance of deep diagonal-circulant neural networks in different settings with respect to other state-of-the-art approaches.
Finally, <Ref> provides concluding remarks.

§ DIAGONAL AND CIRCULANT MATRICES FOR MATRIX DECOMPOSITION

As seen in the Background (<Ref>), circulant matrices exhibit several interesting properties from the perspective of numerical computations.
Most importantly, any $n \times n$ circulant matrix $\Cmat$ can be represented using only $n$ coefficients instead of the $n^2$ coefficients required to represent classical unstructured matrices.
In addition, the matrix-vector product is simplified from $\bigO(n^2)$ to $\bigO(n \log n)$ using the  convolution theorem.
As we will show in this chapter, circulant matrices combined with diagonal matrices can have a strong expressive power.

They can also be used as building blocks to represent any linear transform <cit.> with an arbitrary precision.

We are interested in the relation between the product of diagonal and circulant matrices and the expressivity of low-rank matrices.
<cit.> were able to bound the number of factors that is required to approximate any matrix $\Amat$ with arbitrary precision.
We recall this result in <Ref> as it is the starting point of our theoretical analysis.

For every matrix $\Mmat \in \Cbb^{n \times n}$, there exists a sequence of matrices $\Amat^{(1)} \ldots \Amat^{(2n-1)}$ where $\Amat^{(i)}$ is a circulant matrix if $i$ is odd, and a diagonal matrix otherwise, such that for any $\epsilon > 0$, we have

\begin{equation}
\norm{\Amat^{(1)} \ldots \Amat^{(2n-1)} - \Mmat}_\fro < \epsilon \enspace.
\end{equation}

Unfortunately, this theorem is of little use to understand the expressive power of diagonal-circulant matrices when they are used in deep neural networks for two reasons:

* the bound only depends on the dimension of the matrix $\Mmat$, not on the matrix itself;

* the theorem does not provide any insights regarding the expressive power of $m$ diagonal-circulant factors when $m$ is much lower than $2n - 1$ as it is the case in most practical scenarios we consider in this chapter.

In the following theorem, we improve the result of <cit.> by expressing the number of factors required to approximate $\Mmat$, as a function of the rank of $\Mmat$.
This is useful when one deals with low-rank matrices, which is common in machine learning problems.
Note that in this chapter, our results hold for complex matrices.
This is due to the fact that they are based on <Ref> which holds for complex matrices.

[Rank-based diagonal-circulant decomposition]
Let $\Mmat \in \Cbb^{n \times n}$ be a matrix of rank at most $k$.
Assume that $n$ can be divided by $k$.
There exists a sequence of $4k+1$ matrices $\Amat^{(1)} \ldots \Amat^{(4k+1)}$, where $\Amat^{(i)} \in \Cbb^{n \times n}$ is a circulant matrix if $i$ is odd, and a diagonal matrix otherwise, such that for any $\epsilon > 0$, we have

\begin{equation}
\norm{\Amat^{(1)} \ldots \Amat^{(4k+1)} - \Mmat}_{\fro} < \epsilon \enspace.
\end{equation}

Let $\Umat \mathbf{\Sigma} \Vmat^*$ be the SVD decomposition of $\Mmat$ where $\Umat,\Vmat$ and $\mathbf{\Sigma}$ are $n \times n$ matrices.
Because $\Mmat$ is of rank $k$, the last $n-k$ columns of $\Umat$ and $\Vmat$ are null.
In the following, we will first decompose $\Umat$ into a product of matrices $\Wmat\Rmat\Omat$, where $\Rmat$ and $\Omat$ are respectively circulant and diagonal matrices, and $\Wmat$ is a matrix which will be further decomposed into a product of diagonal and circulant matrices.
Then, we will apply the same decomposition technique to $\Vmat$.
Ultimately, we will get a product of $4k+1$ matrices alternatively diagonal and circulant.

Let $\Rmat = \circulant(r_1\ldots r_n)$. Let $\Omat$ be an $n \times n$ diagonal matrix where $\leftmat \Omat \rightmat_{i,i} = 1$ if $i \le k$ and $0$ otherwise.
The $k$ first columns of the product $\Rmat\Omat$ will be equal to that of $\Rmat$, and the $n-k$ last columns of $\Rmat\Omat$ will be zeros. For example, if $k=2$, we have:

\begin{equation}
\Rmat\Omat = \leftmatrix
r_{1} & r_{n} & 0 & \cdots & 0\\
r_{2} & r_{1}\\
r_{3} & r_{2} & \vdots &  & \vdots\\
\vdots & \vdots\\
r_{n} & r_{n-1} & 0 & \cdots & 0
\rightmatrix
\end{equation}

Let us define $k$ diagonal matrices $\Dmat^{(i)} = \diagonal(d_1^{(i)} \ldots d_n^{(i)})$ for $i \in [k]$.
For now, the values of $d_{j}^{(i)}$ are unknown, but we will show how to compute them.
Let $\Zmat_{1, k}$ and $\Zmat_{1, n}$ be $1$-unit-circulant matrix respectively of size $k \times k$ and $n \times n$ as defined in <Ref> and let $\Wmat = \sum_{i=1}^{k} \Dmat^{(i)} \Zmat_{1,n}^{i-1}$.

Note that the $n-k$ last columns of the product $\Wmat\Rmat\Omat$ will be zeros.
For example, with $k=2$, we have:

\begin{equation}
\Wmat = \leftmatrix
d_{1}^{(1)} &  &  &  & d_{1}^{(2)} \\
d_{2}^{(2)} & d_{2}^{(1)} \\
& d_{3}^{(2)} & \ddots \\
&  & \ddots & \ddots \\
&  &  & d_{n}^{(2)} & d_{n}^{(1)}
\rightmatrix
\end{equation}

\begin{equation}
\Wmat\Rmat\Omat = \leftmatrix
r_{1} d_{1}^{(1)} + r_{n}   d_{1}^{(2)} & r_{n}   d_{1}^{(1)} + r_{n-1} d_{1}^{(2)} & 0 & \cdots & 0 \\
r_{2} d_{2}^{(1)} + r_{1}   d_{2}^{(2)} & r_{1}   d_{2}^{(1)} + r_{n}   d_{2}^{(2)} & 0 & \cdots & 0 \\
\vdots & \vdots & \vdots &  & \vdots \\
r_{n} d_{n}^{(1)} + r_{n-1} d_{n}^{(2)} & r_{n-1} d_{n}^{(1)} + r_{n-2} d_{n}^{(2)} & 0 & \cdots & 0
\rightmatrix
\end{equation}

We want to find the values of $d_{j}^{(i)}$ such that $\Wmat \Rmat \Omat = \Umat$.
We can formulate this as a linear equation system.
In case $k=2$, we get:

\begin{equation}
\leftmatrix
r_{n} & r_{1}\\
r_{n-1} & r_{n}\\
&  & r_{1} & r_{2}\\
&  & r_{n} & r_{1}\\
&  &  &  & r_{2} & r_{3}\\
&  &  &  & r_{1} & r_{2}\\
&  &  &  &  &  & \ddots\\
&  &  &  &  &  &  & \ddots
\rightmatrix \times \leftmatrix
d_{1}^{(2)} \\
d_{1}^{(1)} \\
d_{2}^{(2)} \\
d_{2}^{(1)} \\
d_{3}^{(2)} \\
d_{3}^{(1)} \\
\vdots\\
\vdots
\rightmatrix = \leftmatrix
(\Umat)_{0,0} \\
(\Umat)_{0,1} \\
(\Umat)_{1,0} \\
(\Umat)_{1,1} \\
\\
\\
\vdots\\
\\
\rightmatrix
\end{equation}

The $i^{th}$ block of this block-diagonal matrix is a Toeplitz matrix induced by a contiguous subsequence of length $k+1$ of $(r_1,\ldots r_n,r_1 \ldots r_n)$.
Set $r_{j}=1$ for all $j\in\{k,2k,3k,\ldots n\}$ and set $r_{j}=0$ for all other values of $j$.
Then it is easy to see that each block is equal to $\Zmat^\alpha_{1,k}$ for some $\alpha$.

Note that the matrices $\Zmat^\alpha_{1,k}$ are invertible.

This entails that the block diagonal matrix above is also invertible.
So by solving this set of linear equations, we could find $d_1^{(1)} \ldots d_n^{(k)}$ such that $\Wmat\Rmat\Omat=\Umat$.
We can apply the same idea to factorize $\Vmat = \Wmat' \Rmat \Omat$ for some matrix $\Wmat'$.
Finally, we get

\begin{equation}
\Amat = \Umat \mathbf{\Sigma} \Vmat^* = \Wmat\Rmat\Omat \mathbf{\Sigma} \Omat^* \Rmat^* \Wmat^{'*}
\end{equation}

Note that the matrix $\Wmat\Rmat$ can be decomposed as follows:

\begin{equation}
\Wmat\Rmat = \left( \sum_{i=1}^{k} \Dmat^{(i)} \Zmat_{1,n}^{(i-1)} \right) \Rmat
\end{equation}

where the last matrix $\Zmat_{1,n}^{(k-1)} \Rmat$ is a circulant matrix because both matrices are circulant (see <Ref>).
The same reasoning can be applied with the matrix $\Rmat^*\Wmat'^*$.
Therefore, by construction, the matrices $\Wmat\Rmat$ and $\Rmat^*\Wmat'^*$ can both be factorized by $2k$ circulant and diagonal matrices.
Also note that $\Omat \mathbf{\Sigma} \Omat^*$ is a diagonal matrix, because $\Omat$ and $\mathbf{\Sigma}$ are diagonal matrices.
Overall, $\Amat$ can be represented with a product of $4k+1$ matrices, alternatively diagonal and circulant.

A direct consequence of <Ref>, is that if the number of diagonal-circulant factors is set to a value $K$, we can represent all linear transforms $\Mmat$ whose rank is $\frac{K - 1}{4}$.
Compared to <cit.>, this result shows that structured matrices with fewer than $2n$ diagonal-circulant matrices (as it is the case in practice) can still represent a large class of matrices.

In the following section, we will analyze the expressivity of neural networks based on diagonal and circulant matrices.
In order to characterize the expressivity, we will decompose the matrices of a dense neural network with diagonal and circulant matrices based on <Ref>.

§ ANALYSIS OF DIAGONAL CIRCULANT NEURAL NETWORKS

<cit.> have shown that circulant networks with 2 layers and unbounded width are universal approximators.
However, results on unbounded networks offer weak guarantees and two important questions have remained open until now:

* Can we approximate any function with a bounded-width diagonal-circulant network?

* What function can we approximate with a diagonal-circulant neural network that has a bounded width and a small depth?

We answer these two questions in this section.

First, we present two lemmas that establish a link between the matrix decomposition presented in <Ref> and DCNNs and allow us to present our answer to the first question (<Ref>).
Then, we analyze the expressive power of small depth diagonal-circulant neural networks by comparing them to dense neural networks.
As in the previous section, we still work in the complex domain.

Therefore, we need to extend the definition of neural networks to the complex domain.
First, let us introduce an extension of the function.

Let us define the complex function $\rho: \Cbb^n \rightarrow \Cbb^n$ by: $\rho(\zvec) = \max\left(0, \mathfrak{R}(\zvec)\right) + \ci \max\left(0, \mathfrak{I}(\zvec) \right)$

Given a depth $\depth \in \Nbb$,
let us define $\Omega = \left\{ \left( \Wmat^{(i)}, \bvec^{(i)} \right) \right\}_{i \in [\depth]}$ a set of weights matrices and bias vectors
such that $\Wmat^{(i)} \in \Cbb^{n \times n}$ and $\bvec^{(i)} \in \Cbb^n$.
Let $\Xset \subset \Cbb^n$ and $\Yset \subset \Cbb^n$ be the input space and output space respectively.
A dense neural network is a function $\nn_\Omega: \Xset \rightarrow \Yset$ such that

\begin{equation}
\nn^\act_\Omega (\xvec) \triangleq \layer^\act_{\Wmat^{(\depth)}, \bvec^{(\depth)}} \circ \cdots \circ \layer^\act_{\Wmat^{(1)}, \bvec^{(1)}}(\xvec)
\end{equation}

where $\rho$ is the complex function, $\layer^{\act}_{\Wmat^{(i)},\bvec^{(i)}}: \Cbb^n \rightarrow \Cbb^n$ is a layer parameterized by the weight matrix $\Wmat^{(i)}$ and the bias vector $\bvec^{(i)}$, which can be expressed as follows:

\begin{equation}
\layer^{\act}_{\Wmat^{(i)},\bvec^{(i)}} (\xvec) \triangleq \act \left(\Wmat^{(i)}\xvec + \bvec^{(i)}\right) \enspace,
\end{equation}

Given a depth $\depth \in \Nbb$,
let us define $\Pi = \left\{ \left( \Dmat^{(i)}, \Cmat^{(i)}, \bvec^{(i)} \right) \right\}_{i \in [\depth]}$ a set of weight matrices and bias vectors
such that $\Dmat^{(i)} \in \Cbb^{n \times n}$ is diagonal, $\Cmat^{(i)} \in \Cbb^{n \times n}$ is circulant and $\bvec^{(i)} \in \Cbb^n$.
Let $\Xset \subset \Cbb^n$ and $\Yset \subset \Cbb^n$ be the input space and output space respectively.
Let us denote the product of $\Dmat^{(i)}$ and $\Cmat^{(i)}$ by $\Dmat\Cmat^{(i)}$.
A diagonal-circulant neural network is a function $\nn_\Pi: \Xset \rightarrow \Yset$ such that

\begin{equation}
\nn^\act_\Pi (\xvec) \triangleq \layer^\act_{\Dmat\Cmat^{(\depth)}, \bvec^{(\depth)}} \circ \cdots \circ \layer^\act_{\Dmat\Cmat^{(1)}, \bvec^{(1)}}(\xvec)
\end{equation}

where $\layer^{\act}_{\Dmat\Cmat^{(i)},\bvec^{(i)}}: \Cbb^n \rightarrow \Cbb^n$ is a layer parameterized by the weight matrix $\Dmat\Cmat^{(i)}$, the bias vector $\bvec^{(i)}$ and can be expressed as follows:

\begin{equation}
\layer^{\act}_{\Dmat\Cmat^{(i)},\bvec^{(i)}} (\xvec) \triangleq \act \left(\Dmat\Cmat^{(i)}\xvec + \bvec^{(i)}\right) \enspace,
\end{equation}

where $\rho$ is the complex function.

Diagonal-circulant neural networks are compact due to the layer being parameterized by diagonal and circulant matrices.
Indeed, diagonal and circulant matrices of size  $n \times n$ can be represented with only $n$ values.
Therefore, the layer $\layer^{\act}_{\Dmat\Cmat^{(i)},\bvec^{(i)}}$ is parameterized by $3n$ complex values.

Diagonal-circulant neural networks can have more parameters than a dense neural networks but their depth need to be scaled accordingly.
Let $\depth_1$ and $\depth_2$ be the depth of a dense neural network and a diagonal-circulant neural network respectively, then $\depth_2$ needs to be higher than $\depth_1 \frac{n+1}{3}$ to have more parameters than the dense network.

§.§ From Matrix Decomposition to Neural Networks

The purpose of this section is to extend the matrix decomposition presented in <Ref> to neural networks (<Ref>) and show that bounded-width diagonal-circulant neural networks can approximate any dense neural network (<Ref>).

Let $\Wmat^{(1)} \ldots \Wmat^{(\depth)} \in \Cbb^{n\times n}$, $\bvec \in \Cbb^n$ and let $\Xset \subset \Cbb^n$ be a bounded set.
There exists $\cvec^{(1)} \ldots \cvec^{(\depth)} \in \Cbb^n$ such that for all $\xvec \in \Xset$ we have

\begin{equation}
\rho\left(\Wmat^{(\depth)} \ldots \Wmat^{(1)} \xvec + \bvec \right) = \layer^\rho_{\Wmat^{(\depth)},\cvec^{(\depth)}} \circ \ldots \circ \layer^\rho_{\Wmat^{(1)},\cvec^{(1)}}(\xvec)
\end{equation}

where $\layer^\rho_{\Wmat^{(i)},\cvec^{(i)}} = \rho( \Wmat^{(i)} \xvec + \cvec^{(i)} )$ and $\rho$ is the complex function.

Let $\pmb{\mathcal{W}}(j) = \prod_{k=1}^{j} \Wmat^{(k)}$ and let us define the following set:

\begin{equation}
\mathcal{S} = \left\{ \left(\pmb{\mathcal{W}}(j) \xvec \right)_{t} \mid \ \xvec \in \Xset, t \in [n], j \in [\depth] \right\}
\end{equation}

and let $\Xi = \max\left\{ |\mathfrak{R}(v)|: v \in \Sset \right\} + \ci \max\left\{ |\mathfrak{I}(v)|:v \in \Sset \right\}$.
Intuitively, the real and imaginary parts of $\Xi$ are the largest any activation in the network can have.
Define $\psi_{\Wmat^{(i)}, \cvec^{(i)}}(\xvec) = \Wmat^{(i)} \xvec + \cvec^{(i)}$. Let $\cvec^{(1)} = \Xi \onevec{n}$.
Clearly, for all $\xvec \in \Xset$ we have $\psi_{\Wmat^{(1)}, \cvec^{(1)}}(\xvec) \ge 0$, so $\act(\psi_{\Wmat^{(1)}, \cvec^{(1)}}(\xvec)) = \psi_{\Wmat^{(1)}, \cvec^{(1)}}(\xvec)$ where $\act$ is the complex function.
More generally, for all $j < \depth-1$ define $\cvec^{(j+1)} = \onevec{n} \Xi - \Wmat^{(j+1)} \cvec^{(j)}$.
It is easy to see that for all $j < \depth$ we have

\begin{equation}
\psi_{\Wmat^{(j)}, \cvec^{(j)}}\circ \ldots \circ \psi_{\Wmat^{(1)}, \cvec^{(1)}}(\xvec) = \pmb{\mathcal{W}}(j) \xvec + \onevec{n} \Xi \enspace.
\end{equation}

This guarantees that for all $j < \depth$,

\begin{equation}
\psi_{\Wmat^{(j)}, \cvec^{(j)}} \circ \ldots \circ \psi_{\Wmat^{(1)}, \cvec^{(1)}}(\xvec) = \act \circ \psi_{\Wmat^{(j)}, \cvec^{(j)}} \circ \ldots \circ \act \circ \psi_{\Wmat^{(1)}, \cvec^{(1)}}(\xvec) \enspace.
\end{equation}

Finally, define $\cvec^{(\depth)} = \bvec - \Wmat^{(\depth)} \cvec^{(\depth-1)}$.
We have,

\begin{equation}
\act \circ \psi_{\Wmat^{(\depth)}, \cvec^{(\depth)}} \circ \ldots \circ \act \circ \psi_{\Wmat^{(1)}, \cvec^{(1)}}(\xvec) = \act \big( \pmb{\mathcal{W}}(\depth) \xvec + \bvec \big) \enspace,
\end{equation}

which concludes the proof.

The following lemma is our first result on the expressivity of diagonal-circulant neural networks.
It states that a diagonal-circulant neural network with bounded width and depth can approximate any dense neural network.
To prove this result, we use the matrix decomposition from <Ref> and <Ref> to decompose the dense matrices of the layers of a dense network and unfold it.

Let $\nn_\Omega$ be a dense neural network of width $n$ and depth $\depth$,
and let $\Xset \subset \Cbb^n$ be a bounded set.
There exists a diagonal-circulant neural network $\nn_\Pi$ of width $n$ and of depth $(2n-1)\depth$ such that for any $\epsilon > 0$, we have

\begin{equation}
\norm{\nn_\Omega(\xvec) - \nn_\Pi(\xvec)}_2 < \epsilon, \quad \forall \xvec \in \Xset \enspace.
\end{equation}

Let us assume $\nn_\Omega = \layer_{\Wmat^{(\depth)}, \bvec^{(\depth)}} \circ \ldots \circ \layer_{\Wmat^{(1)}, \bvec^{(1)}}$.
By <Ref>, for any $\epsilon' > 0$, any matrix $\Wmat^{(i)}$, there exists a sequence of $2n-1$ diagonal, $\{\Dmat^{(i, j)}\}_{i \in [\depth], j \in [2n-1]}$, and circulant matrices, $\{\Cmat^{(i, j)}\}_{i, \in [\depth], j \in [2n-1]}$, such that for all $i \in [\depth]$,

\begin{equation}
\norm{\prod_{j=1}^{2n-1} \Dmat^{(i,2n-j)} \Cmat^{(i,2n-j)} - \Wmat^{(i)}}_\fro < \epsilon' \enspace.
\end{equation}

For simplicity, let us denote the product of the two matrices $\Dmat^{(i,j)} \Cmat^{(i,j)}$ by $\Dmat\Cmat^{(i,j)}$.
By <Ref>, we know that there exists a sequence of bias vectors $\left\{ \cvec^{(i,j)} \right\}_{i \in [\depth], j \in [2n-1]}$ such that for all $i\in[\depth]$,

\begin{equation}
\layer^\act_{\Dmat\Cmat^{(i, 2n-1)}, \cvec^{(i, 2n-1)}} \circ \ldots \circ \layer^\act_{\Dmat\Cmat^{(i, 1)}, \cvec^{(i, 1)}}(\xvec) = \act \left(\Dmat\Cmat^{(i,2n-1)} \ldots \Dmat\Cmat^{(i,1)} \xvec + \bvec^{(i)} \right).
\end{equation}

Now if $\epsilon'$ tends to zero,

\begin{equation}
\norm{ \layer^\act_{\Dmat\Cmat^{(i,2n-1)},\cvec^{(i,2n-1)}} \circ \ldots \circ \layer^\act_{\Dmat\Cmat^{(i,1)}, \cvec^{(i,1)}} - \act \left(\Wmat^{(i)} \xvec + \bvec^{(i)} \right)}_2
\end{equation}

will also tend to zero for any $\xvec \in \Xset$, because the function is continuous and $\Xset$ is bounded.
Let $\nn_\Pi = \layer^\act_{\Dmat\Cmat^{(\depth,2n-1)},\cvec^{(\depth,2n-1)}} \circ \ldots \circ \layer^\act_{\Dmat\Cmat^{(1,1)}, \cvec^{(1,1)}}$, because all functions are continuous, for all $\xvec \in \Xset$, $\norm{\nn_\Omega(\xvec) - \nn_\Pi(\xvec)}_2$ tends to zero as $\epsilon'$ tends to zero which concludes the proof.

Now that we know that diagonal-circulant neural networks can approximate any dense neural networks with arbitrary precision, we can extend is result to any function, thus demonstrating that they are universal approximators.
First, let us present universal approximation results for neural networks.
<cit.> have shown that neural networks with a single hidden layer and sigmoid activation can approximate any function if the hidden layer is allowed to be arbitrary large.
However, arbitrary large neural networks lack practical applications.

More recently, the universal approximation results have been extended to bounded width neural networks with arbitrary depth <cit.>.
More formally, we have the following result for neural networks with $\relu$ activations:

For any continuous function $f:[0,1]^{n} \rightarrow \Rbb_+$ of bounded supremum norm, for any $\epsilon>0$, there exists a neural network $\nn_\weights$ parameterized by $\weights$ with an input layer of width $n$, an output layer of width $1$, hidden layers of width $n+3$ and ReLU activations such that

\begin{equation}
\forall x \in [0,1]^n, \quad \left| f(\xvec) - \nn_\weights(\xvec) \right| < \epsilon \enspace.
\end{equation}

From <Ref> and <Ref> by <cit.> which states that dense neural networks are universal approximators, we can prove that bounded-width diagonal-circulant neural networks are also universal approximators.

Diagonal circulant neural networks with bounded width are universal approximators in the following sense:
for any continuous function $f:[0,1]^n \rightarrow \Rbb_+$ of bounded supremum norm, for any $\epsilon > 0$, there exists a complex-valued diagonal-circulant neural network $\nn_\Pi$ of width $n+3$ such that $\forall \xvec \in [0,1]^{n+3}$, $\left| f(\xvec_{1} \ldots \xvec_{n}) - \left( \nn_\Pi \left( \xvec \right) \right)_0 \right| < \epsilon$ where $|\cdot|$ refer to the complex modulus.

From <Ref>, we know that there exists a dense neural network $\nn_\Omega$ with an input layer of width $n$, an output layer of width $1$, hidden layers of width $n+3$ and ReLU activations such that $\forall \xvec \in [0,1]^n, \left| f(\xvec) - \nn_\Omega (\xvec) \right| < \epsilon$.
From $\nn_\Omega$, we can easily build a dense neural networks $\widetilde{\nn}_\Omega$ of width exactly $n+3$, such that $\forall \xvec \in [0,1]^{n+3}$, $\left| f(\xvec_{1} \ldots \xvec_{n}) - \left(\widetilde{\nn}_\Omega (\xvec) \right)_{0}\right| < \epsilon$.
Thanks to <Ref>, this last network can be approximated arbitrarily well by a diagonal-circulant neural network of width $n+3$.
Note that the matrices in the diagonal-circulant neural network are complex, even though we are approximating a real-valued function.

The previous result shows that diagonal-circulant neural networks are universal approximators of real-valued functions.
However the depth needed is in $\bigO(n)$ where $n$ is the width of the network (size of the input).
The depth needed to reach universal approximation is not small, in our experiments, $n$ can be over 3000.
Nonetheless, <cit.> have provided empirical evidence that diagonal-circulant neural networks with small depth can offer good performance.
In the following subsection, we study the theoretical expressivity of diagonal-circulant neural networks with bounded-width and small depth.
This study allows us to better understand why DCNNs show good empirical performances with limited depth.

§.§ The Expressive Power of Diagonal-Circulant Neural Networks

In this subsection, we study the expressive power of diagonal-circulant neural networks with small depth.
To assess the expressivity of DCNNs, we compare the depth needed to approximate dense neural networks with low total rank which we define as the sum of ranks of each weights matrix.
With the concept of total rank, we present in the following,  our result on the expressive power of DCNNs with respect to the total rank of dense neural networks.

The total rank $\kappa(\nn_{\Omega})$ of the neural network $\nn_{\Omega}$ corresponds to the sum of the ranks of the matrices $\Wmat^{(1)} \ldots \Wmat^{(\depth)}$ as follows

\begin{align}
\kappa(\nn_{\Omega}) \triangleq \sum_{i \in [\depth]} \rank(\Wmat^{(i)}) \enspace.
\end{align}

[Rank-based expressive power of DCNNs]
Let $\nn_\Omega$ be a dense neural network of width $n$, depth $\depth$ and a total rank $K$, and assume $n$ is a power of $2$.
Let $\Xset \subset \Cbb^n$ be a bounded set.
Then, for any $\epsilon > 0$, there exists a diagonal-circulant neural network $\nn_\Pi$ of width $n$ such that $\norm{\nn_\Omega(\xvec) - \nn_\Pi(\xvec)}_2 < \epsilon$ for all $\xvec \in \Xset$ and the depth of $\nn_\Pi$ is bounded by $5K$.

Let $\nn_\Omega$ be a dense neural networks parameterized by $\Omega = \{( \Wmat^{(i)}, \bvec^{(i)} \}_{i \in [\depth]}$ of width $n$, depth $\depth$.
Let $k^{(1)} \ldots k^{(\depth)}$ be the ranks of matrices $\Wmat^{(1)} \ldots \Wmat^{(\depth)}$, which are $n \times n$ matrices. Assume that $\forall i$, $n$ can be divided by $k^{i}$.

By <Ref>, for any $\epsilon > 0$, any matrix $\Wmat^{(i)}$ of rank $k^{(i)}$, there exists a sequence of diagonal matrices $\{\Dmat^{(i, j)}\}_{i \in [\depth], j \in [4k^{i}+1]}$ and circulant matrices, $\{\Cmat^{(i, j)}\}_{i, \in [\depth], j \in [4k^{i}+1]}$, such that for all $i \in [\depth]$,

\begin{equation}
\norm{\prod_{j=1}^{4k^{i}+1} \Dmat^{(i,4k^{i}+2-j)} \Cmat^{(i,4k^{i}+2-j)} - \Wmat^{(i)}}_\fro < \epsilon' \enspace.
\end{equation}

Using the exact same technique as in <Ref>, we can build a diagonal-circulant neural network $\nn_\Pi$, such that

\begin{equation}
\norm{ \nn_\Omega (\xvec) - \nn_\Pi (\xvec)}_2 < \epsilon, \quad \forall \xvec \in \Xset,
\end{equation}

for which the total number of layers is bounded as follows:

\begin{equation}
\sum_{i \in [\depth]} \left( 4 k^{(i)} + 1 \right) \le \depth + 4 \sum_{i \in [\depth]} k^{(i)} \le \depth + 4 \kappa(\nn_\Omega) \le 5 \kappa(\nn_\Omega) \enspace.
\end{equation}

where $\kappa(\nn_\Omega)$ is the total rank of the dense neural network $\nn_\Omega$.

Remark that in the theorem, we require that $n$ is a power of $2$.
We conjecture that the result still holds even without this condition.
This result refines <Ref> and answers our second question: does DCNN of bounded width and small depth can approximate a dense neural network of low total rank?
Note that the converse is not true because an $n \times n$ circulant matrix can be of full rank, therefore approximating a DCNN of depth $1$ can require a dense network of total rank equal to $n$.

Finally, what if we choose to use diagonal-circulant networks with a small depth to approximate a dense neural network whose matrices are not of lower rank?
To answer this question, we present three results. First, we characterize the negative impact of replacing matrices by their low rank approximation.
Then, we extend this result to neural networks and bound the error between a dense neural network with full total rank and one with low total rank.
Finally, <Ref> presents our result which bounds the error between a dense neural network with full total rank and a diagonal-circulant neural network.

Let $\Wmat \in \Cbb^{n \times n}$ with singular values $\sigma_1 \ldots \sigma_n$, and let $\bvec, \xvec, \yvec \in \Cbb^n$.
Let $\widetilde{\Wmat}$ be the matrix obtained by an SVD approximation of rank $k$ of matrix $\Wmat$.
Then we have:

\begin{equation}
\norm{ \act \big( \Wmat\xvec + \bvec \big) - \act \big( \widetilde{\Wmat}\yvec + \bvec \big)}_2 \le \sigma_{1} \norm{\xvec - \yvec}_2 + \sigma_{k+1} \norm{\xvec}_2
\end{equation}

Let us denote $\sigma_j$ be the $j^{th}$ singular value of $\Wmat$ and recall that $\sigma_1(\Wmat) = \norm{\Wmat}_2$ by the definition of the spectral norm.
Furthermore, we have $\sigma_1(\Wmat) = \sigma_1(\widetilde{\Wmat})$ because the greatest singular values are equal for both $\Wmat$ and $\widetilde{\Wmat}$.

Also, note that $\norm{\Wmat - \widetilde{\Wmat}}_{2} = \sigma_{k+1}$.
First, let us bound the formula without ReLUs:

\begin{align}
\norm{\big(\Wmat\xvec+\bvec\big) - \big(\widetilde{\Wmat} \yvec + \bvec \big)}_2 &= \norm{\Wmat\xvec - \widetilde{\Wmat} \xvec - \widetilde{\Wmat} \big(\yvec - \xvec \big)}_2 \\
&\le \norm{\big(\Wmat - \widetilde{\Wmat}\big)\xvec}_2 + \norm{\widetilde{\Wmat}}_{2} \norm{\xvec - \yvec}_2 \\
&\le \norm{\xvec}_2 \sigma_{k+1} + \sigma_1 \norm{\xvec - \yvec}_2
\end{align}

Finally, it is easy to see that for any pair of vectors $\xvec,\yvec \in \Cbb^n$, we have

\begin{equation}
\norm{ \act(\xvec) - \act(\yvec)}_2 \le \norm{\xvec - \yvec}_2 \enspace,
\end{equation}

because the complex function is $1$-Lipschitz.
This concludes the proof.

The lemma above bound the error between a linear transform and its equivalent low rank approximation.
In the following, we extend this result to neural networks.

Let $\nn_\Omega: \Cbb^n \rightarrow \Cbb^n$ be a dense neural network, with activation, parameterized by $\Omega = \left\{ \big(\Wmat^{(i)},\bvec^{(i)} \big) \right\}_{i \in [\depth]}$ with $\Wmat^{(i)} \in \Cbb^{n \times n}, \bvec^{(i)} \in \Cbb^n$ for all $i \in [\depth]$ and $\nn_\Omega = \layer_{\Wmat^{(\depth)},\bvec^{(\depth)}} \circ \ldots \circ \layer_{\Wmat^{(1)},\bvec^{(1)}}$ of depth $\depth$ and width $n$.
Let $\widetilde{\Omega} = \left\{ \big( \widetilde{\Wmat}^{(i)},\bvec^{(i)} \big) \right\}_{i \in [\depth]}$ where $\widetilde{\Wmat}^{(i)}$ is the matrix obtained by the SVD approximation of rank $k$ of matrix $\Wmat^{(i)}$.
Define the network $\nn_{\widetilde{\Omega}}$ and let $\sigma_{j}^{(i)}$ be the $j^{th}$ singular value of $\Wmat^{(i)}$ and denote $\sigma_j^{(\max)} = \max_i \sigma_j^{(i)}$, the largest $j^{th}$ singular value across layers.
Then, for any $\xvec \in \Cbb^n$, we have:

* if $\sigma_{1}^{(\max)} = 1$:

\begin{equation}
\norm{ \nn_\Omega(\xvec) - \nn_{\widetilde{\Omega}}(\xvec)}_2 \le \depth \left(  R \sigma_{k+1}^{(\max)} \right) \enspace.
\end{equation}

* if $\sigma_{1}^{(\max)} \neq 1$:

\begin{equation}
\norm{ \nn_\Omega(\xvec) - \nn_{\widetilde{\Omega}}(\xvec)}_2 \le \frac{\left( \big(\sigma_{1}^{(\max)}\big)^\depth - 1 \right) R \sigma_{k+1}^{(\max)}}{\sigma_{1}^{(\max)} - 1}
\end{equation}

where $R$ is an upper bound on the norm of the output of any layer in $\nn_\Omega$.

Let $\xvec^{(0)} \in \Cbb^n$ and $\yvec^{(0)} = \xvec^{(0)}$.
For all $i \in [\depth]$, define $\xvec^{(i)} = \act \left(\Wmat^{(i)} \xvec^{(i-1)} + \bvec^{(i)} \right)$ and $\yvec^{(i)} = \act \left( \widetilde{\Wmat}^{(i)} \yvec^{(i-1)} + \bvec^{(i)} \right)$.
We aim to upper bound the difference in norm of $\xvec^{(i)}$ and $\yvec^{(i)}$.
First, let us consider the linear transform within $\xvec^{(i)}$ and $\yvec^{(i)}$:
The difference in norm between $\xvec^{(i)}$ and $\yvec^{(i)}$ can be upper bounded as follows:

\begin{align}
\norm{ \xvec^{(i)} - \yvec^{(i)}}_2 &\le \sigma_{1}^{(i)} \norm{ \xvec^{(i-1)} - \yvec^{(i-1)}}_2 + \sigma_{k+1}^{(i)} \norm{ \xvec^{(i-1)}}_2 \\
&\leq \sigma_1^{(\max)} \norm{ \xvec^{(i-1)} - \yvec^{(i-1)}}_2 + \sigma_{k+1}^{(\max)} R
\end{align}

where the first inequality stems from <Ref> and the second by setting $\sigma_j^{(\max)} = \max_i \sigma_j^{(i)}$ where $\sigma_j^{(\max)}$ is the largest $j^{th}$ singular value across layers and $R = \max_i \norm{\xvec^{(i)}}_2$.
From there, we need to consider two cases:

* If $\sigma_1^{(\max)} = 1$: we have a recurrence relation of the form $a_n = a_{n-1} + s$ with $a_0 = 0$ which can unfold as follows: $a_n = ns$.
We can apply this formula to bound our error as follows:

\begin{equation}
\norm{ \xvec^{(\depth)} - \yvec^{(\depth)}}_2 \le \depth \left(  R \sigma_{k+1}^{(\max)} \right) \enspace.
\end{equation}

* If $\sigma_1^{(\max)} \neq 1$: we have a recurrence relation of the form $a_n = r a_{n-1} + s$ with $a_0= 0$ which can unfold as follows: $a_n = \frac{s \left(r^n - 1\right)}{r-1}$.
We can apply this formula to bound our error as follows:

\begin{equation}
\norm{ \xvec^{(\depth)} - \yvec^{(\depth)} }_2 \le \frac{ \left( \big(\sigma_1^{(\max)}\big)^\depth - 1 \right) R \sigma_{k+1}^{(\max)}}{\sigma_1^{(\max)}-1} \enspace,
\end{equation}

which concludes the proof.

<Ref> bounds the error between a dense neural network and a neural network whose matrices are the low rank approximations of the first ones.
Two cases are presented. If the largest singular value across the network is equal to 1, then the error is polynomial with the depth of the network.
In the case where the largest singular value across the network is different from 1, the error is exponential with respect to the depth of the network.

Now, we can easily extend <Ref> to diagonal-circulant neural networks.
By <Ref>, we can replace the layers with low rank approximation by the product of diagonal and circulant matrices leading to a diagonal-circulant neural network with a higher depth.

Let $\nn_\Omega$ be a dense neural network of depth $\depth$ and width $n$ and parameterized by $\Omega = \left\{(\Wmat^{(i)},\bvec^{(i)} ) \right\}_{i \in [\depth]}$.
Let $\sigma_1^{(i)}$ be the largest singular value of $\Wmat^{(i)}$.
Let $\Xset \subset \Cbb^n$ be a bounded set.
Let $k$ be an integer dividing $n$.
There  exists a diagonal-circulant neural network $\nn_\Pi$ of width $n$ and of depth $m = (4k+1)\depth$, parameterized by $\Pi = \left\{(\Dmat^{(i)}\Cmat^{(i)},\cvec^{(i)})\right\}_{i \in [m]}$ such that, for any $\xvec \in \Cbb^n$, we have:

* if $\sigma_{1}^{(\max)} = 1$:

\begin{equation}
\norm{ \nn_\Omega(\xvec) - \nn_{\Pi}(\xvec)}_2 \le \depth \left(  R \sigma_{k+1}^{(\max)} \right) \enspace.
\end{equation}

* if $\sigma_{1}^{(\max)} \neq 1$:

\begin{equation}
\norm{ \nn_\Omega(\xvec) - \nn_{\Pi}(\xvec)}_2 \le \frac{\left( \big(\sigma_{1}^{(\max)}\big)^\depth - 1 \right) R \sigma_{k+1}^{(\max)}}{\sigma_{1}^{(\max)} - 1} \enspace,
\end{equation}

where $R$ is an upper bound on the norm of the output of any layer in $\nn_\Omega$.

Let $\nn_{\widetilde{\Omega}}$ be a dense neural network of depth $\depth$ and width $n$ and let $\widetilde{\Omega} = \left\{ \big(\widetilde{\Wmat}^{(i)}, \bvec^{(i)} \big) \right\}_{i \in [\depth]}$ such that $\widetilde{\Wmat}^{(i)}$ is the matrix obtained by the SVD approximation of rank $k$ of matrix $\Wmat^{(i)}$.
With <Ref>, we have an error bound on $\norm{ \nn_\Omega (\xvec) - \nn_{\widetilde{\Omega}} (\xvec)}_2$.
Now each matrix $\widetilde{\Wmat}^{(i)}$ can be replaced by a product of $4k+1$ diagonal-circulant matrices.
By <Ref>, this product yields a diagonal-circulant neural network of depth $m = (4k+1)\depth$, strictly equivalent to $\nn_{\widetilde{\Omega}}$ on $\Xset$.
This concludes the proof.

Illustration of the expressivity of diagonal-circulant neural networks.

We highlight the significance of these results with the two following properties.

Let $\Rset_{K}$ be the set of all functions $\nn_\Omega:\Cbb^n \rightarrow \Cbb^n$ for all $n$, representable by a dense neural network with complex activation of total rank at most $K$ and let $\Cset_\depth$ be the set of all functions $\nn_\Pi:\Cbb^n \rightarrow \Cbb^n$ for all $n$, representable by deep diagonal-circulant networks of depth at most $\depth$, then:

\begin{align}
\label{property:eq1} \forall K,\exists \depth \, &\quad \Rset_{K} \subsetneq \Cset_{\depth} \\
\label{property:eq2} \forall \depth,\nexists K\, &\quad \Cset_{\depth} \subseteq \Rset_{K}
\end{align}

We illustrate the meaning of these properties using <Ref>.
As we can see, the set $\mathcal{R}_K$ of all the functions representable by a dense neural network of total rank $K$ is strictly included in the set $\mathcal C_{9K}$ of all diagonal-circulant neural networks of depth $9K$ (as by <Ref>).
These properties are interesting for many reasons.
First, <Ref> shows that diagonal-circulant networks are strictly more expressive than networks with low total rank.
Second and most importantly, in standard deep neural networks, it is known that the most of the singular values are close to zero. This phenomenon is called rank collapse and it has been first observed by <cit.> and confirmed later by <cit.>.
<Ref> shows that these networks can efficiently be approximated by diagonal-circulant networks.
Finally, several publications have shown that neural networks can be trained explicitly to have low-rank weight matrices <cit.>.
This opens the possibility of learning compact and accurate diagonal-circulant networks.

§ HOW TO TRAIN DEEP DIAGONAL CIRCULANT NEURAL NETWORKS?

Training diagonal-circulant neural networks has revealed to be a challenging problem.
Indeed, as discussed earlier, the expressivity of diagonal-circulant neural networks scale with depth.

In the following, we devise two techniques to facilitate the training of deep diagonal-circulant neural networks.
First, we propose an initialization procedure which guarantees the signal is propagated across the network without vanishing nor exploding.
Secondly, we study the behavior of DCNNs with different nonlinearity functions and determine the best parameters for different settings.
Note that we choose to train diagonal-circulant neural networks with real matrices instead of complex ones.
Indeed, the complex version of diagonal-circulant neural networks would have twice the number of parameters.

§.§ Initialization Scheme of Diagonal-Circulant Neural Networks

In order to facilitate the training of deep diagonal-circulant neural networks, we extend the Xavier initialization <cit.> which is an initialization scheme proposed for dense and convolutional neural networks.
First, for each circulant matrix $\Cmat = \circulant (\cvec)$ with $\cvec \in \Rbb^n$, each $\cvec_i$ is randomly drawn from $\mathcal{N} \left(0,\alpha^2\right)$, with $\alpha = \sqrt{\frac{2}{n}}$.
Next, for each diagonal matrix $\Dmat = \diagonal (\dvec)$ with $\dvec \in \Rbb^n$, each $\dvec_{i}$ is drawn randomly and uniformly from $\{-1,1\}$ for all $i$.
Finally, all biases in the network are randomly drawn from $\mathcal{N}\left(0,\alpha'^{2}\right)$, for some small value of $\alpha'$.

Let $\cvec, \dvec, \bvec$ be random variables in $\Rbb^n$ such that $\cvec \sim \mathcal{N}(0, \Imat_n\alpha^2)$, $\bvec \sim \mathcal{N}(0, \Imat_n\alpha'^2)$ and $\dvec_i \sim \{-1, 1\}, \forall i$ uniformly.
Define $\Cmat = \circulant(\cvec)$ and $\Dmat = \diagonal(\dvec)$
Define $\yvec = \Dmat \Cmat \uvec + \bvec$ for some vector $\uvec$ in $\Rbb^n$.
Then, for all $i$, the probability density function () of $\yvec_{i}$ is symmetric.

* Assume $\uvec_0, \ldots, \uvec_{n-1}$ are fixed. Then, we have for $i \neq i'$:

\begin{equation}
\Cov(\yvec_i, \yvec_{i'}) = 0 \quad \text{and} \quad
\Var(\yvec_i) = \alpha'^2 + \sum_{j} \uvec_{j}^{2} \alpha^{2}
\end{equation}

* Let $\xvec$ be random variables in $\Rbb^n$ such that the of $\xvec_i$ is symmetric for all $i$, and let $\uvec_i = \rho(\xvec_i)$.
We have for $i \neq i':$

\begin{equation}
\Cov(\yvec_i, \yvec_{i'}) = 0 \quad \text{and} \quad
\Var(\yvec_i) = \alpha'^2 + \frac{1}{2} \sum_{j} \Var(\xvec_i) \alpha^2
\end{equation}

By an abuse of notation, we write $\cvec_0 = \cvec_n, \cvec_{-1} = \cvec_{n-1}$ and so on.
First, note that: $\yvec_i = \sum_{j=0}^{n-1} \cvec_{j-i} \uvec_j \dvec_j + \bvec_i$.
Observe that the term $\cvec_{j-i}\uvec_j\dvec_j$ has symmetric because of $\dvec_j$.
Thus, $\yvec_i$ has a symmetric .
Now let us compute the covariance.

\begin{align}
\Cov(\yvec_i,\yvec_{i'}) &= \left( \sum_{j,j'=0}^{n-1} \Cov\left(\cvec_{j-i}\uvec_j\dvec_j,\cvec_{j'-i'}\uvec_{j'}\dvec_{j'}\right) \right) + \Cov(\bvec_i, \bvec_{i'}) \\
&= \sum_{j,j'=0}^{n-1} \Ebb\left[\cvec_{j-i}\uvec_j\dvec_j\cvec_{j'-i'}\uvec_{j'}\dvec_{j'}\right] - \Ebb\left[\cvec_{j-i}\uvec_j\dvec_j\right] \Ebb\left[ \cvec_{j'-i'}\uvec_{j'}\dvec_{j'} \right]
\end{align}

Observe that $\Ebb\left[\cvec_{j-i}\uvec_j\dvec_j\right] = \Ebb\left[\cvec_{j-i}\uvec_j\right] \Ebb\left[\dvec_j\right] = 0$ because $\dvec_j$ is independent from $\cvec_{j-i}\uvec_j$.
Also, observe that if $j\neq j'$ then $\Ebb\left[\dvec_j\dvec_{j'}\right]=0$ and

\begin{equation}
\Ebb\left[\cvec_{j-i} \uvec_j \dvec_j \cvec_{j'-i'} \uvec_{j'} \dvec_{j'} \right] = \Ebb\left[ \dvec_j \dvec_{j'} \right] \Ebb\left[ \cvec_{j-i} \uvec_j \cvec_{j'-i'} \uvec_{j'} \right] = 0 \enspace.
\end{equation}

Therefore, the only non null terms are those for which $j = j'$.
We get:

\begin{align}
\Cov(\yvec_i,\yvec_{i'}) & =\sum_{j=0}^{n-1} \Ebb\left[\cvec_{j-i}\uvec_j\dvec_j\cvec_{j-i'}\uvec_j\dvec_j\right] \\
& =\sum_{j=0}^{n-1} \Ebb\left[\cvec_{j-i}\cvec_{j-i'}\uvec_j^{2}\right]
\end{align}

Assume $\uvec$ is a fixed vector.
Then, $\Var(\yvec_i) = \sum_{j=0}^{n-1} \uvec_j^{2} \alpha^2$ and $\Cov(\yvec_i,\yvec_{i'}) = 0$ for $i\neq i'$ because $\cvec_{j-i}$ is independent from $\cvec_{j-i'}$.
Now assume that $\uvec_j = \rho(\xvec_j)$ where $\xvec_j$ is a random variables in $\Rbb^n$.
Clearly, $\uvec_j^2$ is independent from $\cvec_{j-i}$ and $\cvec_{j-i'}$, thus, we have:

\begin{align}
\Cov(\yvec_i,\yvec_{i'}) &= \sum_{j=0}^{n-1} \Ebb \left[ \cvec_{j-i}\cvec_{j-i'}\right] \Ebb\left[ \uvec_j^2 \right] \enspace.
\end{align}

For $i \neq i'$, then $\cvec_{j-i}$ and $\cvec_{j-i'}$ are independent, we have $\Ebb\left[\cvec_{j-i}c_{j-i'}\right] = \Ebb\left[\cvec_{j-i}\right] \Ebb\left[\cvec_{j-i'}\right] = 0$ and $\Cov(\yvec_i,\yvec_{i'}) = 0$ if $i \neq i'$.
Let us compute the variance.
We get $\Var(\yvec_i) = \sum_{j=0}^{n-1} \Var(\cvec_{j-i}) \Ebb\left[\uvec_j^2\right]$.
Because the of $\xvec_j$ is symmetric, $\Ebb\left[\xvec_j^2\right] = 2\Ebb\left[\uvec_j^2\right]$ and $\Ebb\left[\xvec_j\right] = 0$.

\begin{equation}
\Var(\yvec_i) = \frac{1}{2} \sum_{j=0}^{n-1} \Var(\cvec_{j-i}) \Ebb\left[\xvec_j^2\right] = \frac{1}{2} \sum_{j=0}^{n-1} \Var(\cvec_{j-i}) \Var(\xvec_j)
\end{equation}

which concludes the proof.

Now we can state our result on the initialization of diagonal-circulant neural networks.
The following proposition states that the covariance matrix at the output of any layer in a diagonal-circulant neural network is constant.
Moreover, note that the result of this proposition is independent of the depth of the network.

Let $\nn_\Pi$ be a diagonal-circulant neural network of depth $\depth$ initialized according to our procedure, with $\alpha'=0$.
Assume that all layers $1$ to $\depth-1$ have activation functions, and that the last layer has the identity activation function.
Then, for any $\xvec \in \Rbb^n$, the covariance matrix of $\nn_\Pi(\xvec)$ is $\frac{2}{n} \Imat_n \norm{\xvec}_2^2$.

Let $\nn_\Pi \triangleq \layer_{\Dmat^{(\depth)}, \Cmat^{(\depth)}} \circ \ldots \circ \layer_{\Dmat^{(1)}, \Cmat^{(1)}}$ be a $\depth$ layer diagonal-circulant neural network.
All matrices are initialized as described in the statement of the proposition.
Let $\yvec = \Dmat^{(1)} \Cmat^{(1)} \xvec$.
<Ref> shows that $\Cov(\yvec_{i}, \yvec_{i'}) = 0$ for $i \neq i'$ and $\Var(\yvec_{i}) = \frac{2}{n}\norm{x}_{2}^{2}$.
For any $j \le \depth$, define

\begin{equation}
\zvec^{(j)} = \layer_{\Dmat^{(j)}, \Cmat^{(j)}} \circ \ldots \circ \layer_{\Dmat^{(1)}, \Cmat^{(1)}}(\xvec) \enspace.
\end{equation}

By a recursive application of <Ref>, we get

\begin{equation}
\Cov(\zvec_i^{(j)}, \zvec_{i'}^{(j)}) = 0 \quad \text{and} \quad  \Var(\zvec_i^{(j)}) = \frac{2}{n} \norm{\xvec}_2^2
\end{equation}

which concludes the proof.

The effect of <Ref> is that the signal and the gradient will not vanish during the training, facilitating the convergence.
The fact that the result of <Ref> is independent of the depth of the network allows us to train very deep diagonal-circulant neural networks.

§.§ Analysis of the Use of Nonlinearities

Impact of increasing the slope of a Leaky-ReLU in DCNNs.

We empirically found that the ReLU activations had an impact on the training of deep diagonal-circulant neural networks on CIFAR10 dataset.
Indeed, the deeper the network, the more nonlinear it is, which makes convergence difficult.
In an experiment, we replace the ReLU activations with Leaky-ReLU activations (<Ref>) and vary the slope of the Leaky-ReLU (a higher slope means an activation function that is closer to a linear function).
The results of this experiment are presented in <Ref>.
In this experiment, we try different slopes for the Leaky-ReLU activation and train diagonal-circulant neural networks with different depth.
We can observe that a higher slope (making the network more linear) facilitates convergence, allowing us to train deeper networks.
This is an interesting result, since  we can use this technique to adjust the number of parameters in the network (increasing depth), without facing training difficulties.
We hence rely on this setting in the experimental section.

§ EXPERIMENTS

This experimental section aims at answering the following questions:

Q1 How do DCNNs compare to other approaches such as ACDC, LDR or other structured approaches?

Q2 How do DCNNs compare to other compression based techniques?

Q3 How do DCNNs perform in the context of large-scale real-world machine learning applications?

§.§ Comparison with Other Structured Approaches (Q1)

Comparison of the training loss of DCNNs and ACDC networks on a regression task with synthetic data.

Comparison of the training loss of DCNNs and ACDC networks on a CIFAR-10 dataset.

§.§.§ Comparison with ACDC

In <Ref>, we have discussed the differences between the ACDC framework and our approach from a theoretical perspective.
In this section, we conduct experiments to compare the performance of DCNNs with neural networks based on ACDC layers.
We first reproduce the experimental setting from <cit.>, and compare both approaches using only linear networks (, networks without any nonlinear activation).
The synthetic dataset has been created in order to reproduce the experiment on the regression linear problem proposed by <cit.>.
We draw $\Xmat$ and $\Wmat$ from a uniform distribution between [-1, +1] and $\epsilon$ from a normal distribution with mean 0 and variance $0.01$.
The relationship between $\Xmat$ and $\Ymat$ is defined by $\Ymat = \Xmat\Wmat + \epsilon$.
The results are presented in <Ref>.
In this simple setting, while both architectures demonstrate good performance, we can observe that DCNNs offer a better convergence rate.
In <Ref>, we compare neural networks with ReLU activations on CIFAR-10.

We found that networks which are based only on ACDC layers are difficult to train and offer poor accuracy on CIFAR-10 (we have tried different initialization schemes including the one from the original paper, and the one we introduce in this chapter).
<cit.> managed to train a large VGG network  however these networks are generally highly redundant and the contribution of the structured layer is difficult to quantify.
We also observe that adding a single dense layer improves the convergence rate of ACDC in the linear case, which explains the good results of <cit.>.
However, it is difficult to characterize the true contribution of the ACDC layers when the network has a large number of expressive layers.
In contrast, deep DCNNs can be trained and offer good performance without additional dense layers (these results are in line with our experiments on the dataset).

Network size Accuracy compared on Dense networks, DCNNs, DTNNs, neural networks based on Toeplitz matrices and neural networks based on Low Rank-based matrices.

Accuracy of different structured architecture given the number of trainable parameters.

§.§.§ Comparison with Dense Networks, Toeplitz Networks and Low Rank Networks

We now compare DCNNs with other state-of-the-art structured networks by measuring the accuracy on a flattened version of the CIFAR-10 dataset.
Our baseline is a dense feed-forward network with a fixed number of weights (9 million weights).
We compare with DCNNs and with DTNNs (see below), Toeplitz networks, and Low-Rank networks <cit.>.
We first consider Toeplitz networks which are stacked Toeplitz matrices interleaved with ReLU activations since Toeplitz matrices are closely related to circulant matrices.
However, Toeplitz networks have a different structure than DCNNs (they do not include diagonal matrices), therefore, we also experiment using DTNNs, a variant of DCNNs where all the circulant matrices have been replaced by Toeplitz matrices.
Finally we conduct experiments using networks based on low-rank matrices as they are also closely related to our work.
For each approach, we report the accuracy of several networks with a varying depth ranging from 1 to 40 (DCNNs, Toeplitz networks) and from 1 to 30 (from DTNNs).
For low-rank networks, we used a fixed depth network and increased the rank of each matrix from 7 to 40.
We also tried to increase the depth of low rank matrices, but we found that deep low-rank networks are difficult to train so we do not report the results here.
We compare all the networks based on the number of weights from 21K (0.2% of the dense network) to 370K weights (4% of the dense network) and we report the results in <Ref>.
First we can see that the size of the networks correlates positively with their accuracy which demonstrated successful training in all cases.
We can also see that the DCNNs achieves the maximum accuracy of 56% with 20 layers ($\sim$ 200K weights) which is as good as the dense networks with only 2% of the number of weights.
Other approaches also offer good trade-offs but they are not able to reach the accuracy of a dense network.

Architectures     #Parameters     Accuracy

Dense     9.4M	     0.562

DCNN $(5\ layers)$     49K	    0.543

DCNN $(2\ layers)$     21K     0.536

LDR–TD	$(r = 2)$	             64K	    0.511

LDR–TD	$(r = 3)$	             70K	    0.473

Toeplitz-like $(r=2)$	             46K	    0.483

Toeplitz-like $(r =3)$	             52K      0.496

Comparison of LDR and Diagonal-Circulant neural networks on a flattened version of CIFAR-10.

§.§.§ Comparison with LDR networks

We now compare DCNNs with the LDR framework using the network configuration proposed by <cit.>: a single LDR structured layer followed by a dense layer.
In the LDR framework, we can change the size of a network by adjusting the rank of the residual matrix, effectively capturing matrices with a structure that is close to a known structure but not exactly (in the LDR framework, Toeplitz matrices can be encoded with a residual matrix with rank=2, so a matrix that can be encoded with a residual of rank=3 can be seen as Toeplitz-like.).
The results are presented in <Ref> and demonstrate that DCNNs outperform all LDR networks both in terms of size and accuracy.

§.§ Comparison with Other Compression Based Approaches (Q2)

We provide a comparison with other compression based approaches such as HashNet <cit.>, Dark Knowledge <cit.>.
<Ref> shows the test error of DCNN against other known compression techniques on the MNIST datasets.
We can observe that DCNN outperforms HashNet <cit.> and Dark Knowledge <cit.> with fewer number of parameters.

§.§ Large-scale Video Classification on the Dataset (Q3)

1cArchitecture     1c#Params     Error (%)

LeNet <cit.>                 4.2M             0.61

HashNet <cit.>             46K              2.79

Dark Knowledge <cit.>     46K              6.32

DCNN                                  25K     1.74

Comparison with compression based approaches.

To understand the performance of deep DCNNs on large-scale applications, we conducted experiments on the video classification with 3.8 training examples introduced by <cit.>.
This section provides a summary of the results obtained on the dataset, the full experimental analysis is reported in Appendix <ref>.
Notice that we favor this experiment over ImageNet applications because modern image classification architectures involve a large number of convolution layers, and compressing convolution layers is out of our scope.
Also, as mentioned earlier, testing the performance of DCNN architectures mixed with a large number of expressive layers makes little sense.
The includes two datasets describing 8 million labeled videos.
Both datasets contain audio and video features for each video.
In the first dataset (aggregated) all audio and video features have been aggregated every 300 frames.
The second dataset (full) contains the descriptors for all the frames.
To compare the models we use the GAP score (Global Average Precision) proposed by <cit.>.
On the simpler aggregated dataset we compared off-the-shelf DCNNs with a dense baseline with 5.7 million weights.
On the full dataset, we designed three new compact architectures based on the state-of-the-art architecture introduced by <cit.>.

Experiments on the aggregated dataset with DCNNs
We compared DCNNs with a dense baseline with 5.7 million weights.
The goal of this experiment is to discover a good trade-off between depth and model accuracy.
To compare the models we use the GAP score (Global Average Precision) following the experimental protocol proposed by <cit.>, to compare our experiments.
<Ref> shows the results of our experiments on the aggregated dataset in terms of the number of weights and GAP score.
These results suggest that we can compress the baseline at the cost of a little decrease of GAP score.

Architecture     #Weights     GAP@20

Baseline        5.7M      0.773

4 DC            25K       0.599

32 DC           122K      0.685

4 DC + 1 FC     4.46M     0.747

GAP score on the dataset with DCNNs.

Architecture     #Weights     GAP@20

original     45M     0.846

DBoF with DC       36M (80)     0.838

FC with DC        41M (91)     0.845

MoE with DC       12M (26)     0.805

GAP score on the dataset with different layers represented with diagonal-circulant decomposition.

Experiments with DCNNs Deep Bag-of-Frames Architecture:
The Deep Bag-of-Frames architecture can be decomposed into three blocks of layers, as illustrated in <Ref>.
The first block of layers, composed of the Deep Bag-of-Frames embedding (DBoF), is meant to model an embedding of these frames in order to make a simple representation of each video.
A second block of fully connected layers (FC) reduces the dimensionality of the output of the embedding and merges the resulting output with a concatenation operation.
Finally, the classification block uses a combination of Mixtures-of-Experts (MoE) <cit.> and Context Gating <cit.> to calculate the final class probabilities.
<Ref> shows the results in terms of the number of weights, size of the model (MB) and GAP on the full dataset, replacing the DBoF block reduces the size of the network without impacting the accuracy.
We obtain the best compression ratio by replacing the MoE block with DCNNs (26%) of the size of the original dataset with a GAP score of 0.805 (95% of the score obtained with the original architecture).
We conclude that DCNN are both theoretically sound and of practical interest in real, large-scale applications.

Diagram of the state-of-the-art neural network architecture, initially proposed by <cit.> and later improved by <cit.>.

Architectures & Hyper-Parameters:
For the first set of our experiments (experiments on CIFAR-10), we train all networks for 200 epochs, a batch size of 200, Leaky ReLU activation with a different slope.
We minimize the Cross Entropy Loss with Adam optimizer and use a piecewise constant learning rate of $5 \times 10^{-5}$, $2.5\times10^{-5}$, $5\times10^{-6}$ and $1\times10^{-6}$ after respectively 40K, 60K and 80K steps.
For the dataset experiments, we built a neural network based on state-of-the-art architecture initially proposed by <cit.> and later improved by <cit.>.
Remark that no convolution layer is involved in this application since the input vectors are embeddings of video frames processed using state-of-the-art convolutional neural networks trained on ImageNet.
We trained our models with the CrossEntropy loss and used Adam optimizer with a 0.0002 learning rate and a 0.8 exponential decay every 4 million examples.
All fully connected layers are composed of 512 units.
DBoF, NetVLAD and NetFV are respectively 8192, 64 and 64 of cluster size for video frames and 4096, 32, 32 for audio frames.
We used 4 mixtures for the MoE Layer.
We used all the available 300 frames for the DBoF embedding.
In order to stabilize and accelerate the training, we used batch normalization before each nonlinear activation and gradient clipping.

§.§ Exploiting Image Features

Architectures     #Parameters     Accuracy

DC $(1\ layers)$     124K     0.757

DC $(3\ layers)$     217K     0.785

LDR-SD $(r=1)$     140K     0.701

LDR-SD $(r=10)$     420K     0.728

Toeplitz-like $(r=1)$     110K     0.711

Toeplitz-like $(r=10)$     388K     0.720

Accuracy of scattering models followed by LDR or DC layer on CIFAR-10 dataset.

Dense layers and DCNNs are not designed to capture task-specific features such as the translation invariance inherently useful in image classification.
We can further improve the accuracy of such general-purpose architectures on image classification without dramatically increasing the number of trained parameters by stacking them on top of fixed (ie non-trained) transforms such as the scattering transform <cit.>.
In this section we compare the accuracy of various structured networks, enhanced with the scattering transform, on an image classification task, and run comparative experiments on CIFAR-10.

Our test architecture consists of 2 depth scattering on the RGB images followed by a batch norm and LDR or DC layer.
To vary the number of parameters of Scattering+LDR architecture, we increase the rank of the matrix (stacking several LDR matrices quickly exhausted the memory).
The <Ref> and <Ref> show the accuracy of these architectures given the number of trainable parameters.

First, we can see that the DCNN architecture very much benefits from the scattering transform and is able to reach a competitive accuracy over 78%.
We can also see that scattering followed by a DC layer systematically outperforms scattering + LDR or scattering + Toeplitz-like with fewer parameters.

§ CONCLUDING REMARKS

This chapter dealt with the training and understanding of diagonal-circulant neural networks.
To the best of our knowledge, training such networks with a large number of layers had not been done before.
We also endowed this kind of architecture with theoretical guarantees, hence enriching and refining previous theoretical work from the literature.
More importantly, we showed that DCNNs outperform their competing structured alternatives, including the very recent general approach based on LDR networks.
Our results suggest that stacking diagonal-circulant layers with nonlinearities improves the convergence rate and the final accuracy of the network.
Formally proving these statements constitutes the future directions of this work.
We would like to generalize the good results of DCNNs to convolutional neural networks.
We also believe that circulant matrices deserve particular attention in deep learning because of their strong ties with convolutions: a circulant matrix operator is equivalent to the convolution operator with circular padding.
This fact makes any contribution to the area of circulant matrices particularly relevant to the field of deep learning with impacts beyond the problem of designing compact models.

CHAPTER: BOUND ON THE LIPSCHITZ CONSTANT OF CONVOLUTION LAYERS

§ INTRODUCTION

In this chapter we introduce a new upper bound on the largest singular value of convolution layers that is both tight and easy to compute.
Instead of using the power method to iteratively approximate this value, we study the properties of doubly-block Toeplitz matrices and its links with Fourier analysis.

Our work is based on the result of <cit.> which states that an upper bound on the singular value of Toeplitz matrices can be computed from the inverse Fourier transform of the characteristic sequence of these matrices.
We first extend this result to doubly-block Toeplitz matrices (, block Toeplitz matrices where each block is Toeplitz) and then to convolutional operators, which can be represented as stacked sequences of doubly-block Toeplitz matrices.
From our analysis immediately follows an algorithm for bounding the Lipschitz constant of a convolution layer, and by extension the Lipschitz constant of the whole network.
We theoretically study the approximation of this algorithm and show experimentally that it is more efficient and accurate than competing approaches.

Finally, we illustrate our approach on adversarial robustness.
Recent work has shown that empirical methods such as adversarial training (AT) offer poor generalization <cit.>, and can be improved by applying Lipschitz regularization <cit.>.
To illustrate the benefit of our new method, we train a large Wide ResNet with Lipschitz regularization and show that it offers a significant improvement over adversarial training alone, and over other methods for Lipschitz regularization.
In summary, we make the three following contributions:

* We devise an upper bound on the singular values of the operator matrix of convolution layers by leveraging Toeplitz matrix theory and its links with Fourier analysis.

* We propose an efficient algorithm to compute this upper bound which enables its use in the context of Convolutional Neural Networks.

* We use our method to regularize the Lipschitz constant of neural networks for adversarial robustness and show that it offers a significant improvement over AT alone.

§ RESULTS ON THE SPECTRUM OF MATRICES FROM THE TOEPLITZ FAMILY

§.§ Upper-Bounds on the Largest Singular Value of Toeplitz and Block Toeplitz Matrices

Doubly-block Toeplitz matrices inherit the properties of Toeplitz and block Toeplitz matrices.
Recall that for Toeplitz and block Toeplitz matrices, there exist no closed-form expression to compute their eigenvalues.
However, we can represent Toeplitz and block Toeplitz matrices with a 2$\pi$-periodic function which can describe very precisely the spectrum of the matrices.

Let $\{a_h\}_{h \in \Iset_n}$ be the characteristic sequence of a Toeplitz matrix $\Amat \in \Rbb^{n \times n}$ and let $\{\Bmatsf^{(h)}\}_{h \in \Iset_n}$ be the characteristic sequence of $m \times m$ blocks of a block Toeplitz matrix $\Bmat \in \Rbb^{nm \times nm}$ such that $\Amat = (a_{j-i})_{i,j\in \Iset_n}$ and $\Bmat = (\Bmatsf^{(j-i)})_{i,j \in \Iset_n}$ with $\Iset_n = \left\{ -n+1, \cdots, n-1 \right\}$.

Building on the results presented in the Background (<Ref>), we can define two trigonometric polynomials $f: \Rbb \rightarrow \Cbb$ and $F: \Rbb \rightarrow \Cbb^{m \times m}$ as follows:

\begin{equation}
f(\omega) \triangleq \sum_{h \in \Iset_n} a_h e^{\ci h \omega} \quad \quad F(\omega) \triangleq \sum_{h \in \Iset_n} \Bmatsf^{(h)} e^{\ci h \omega} \enspace.
\end{equation}

$f(\omega)$  and $F(\omega)$ are the inverse Fourier transforms of the sequences $\{a_h\}_{h \in \Iset_n}$ and $\{\Bmatsf^{(h)}\}_{h \in \Iset_n}$ respectively.
From there, inspired by the work done by <cit.>, we recall from <Ref> the operator $\Tmat$ mapping integrable functions to Toeplitz matrices:

\begin{equation} \label{equation:expression_toeplitz_matrix}
\Tmat(f) \triangleq\leftmat\frac{1}{2\pi} \int_{0}^{2\pi} e^{-\ci(i-j)\omega} f(\omega) \,\diff \omega \rightmat_{i,j \in \Iset^+_n} \enspace,
\end{equation}

with this operator, we have $\Tmat(f) = \Amat$ and $\Tmat(F) = \Bmat$.

Now, we can state two known theorems which upper bound the largest singular value of Toeplitz and block Toeplitz matrices with respect to their generating functions.

Let $f: \Rbb \rightarrow \Cbb$, be a continuous and $2\pi$-periodic function, then $\Tmat(f) \in \Rbb^{n \times n}$ is a Toeplitz matrix generated by the function $f$
We can bound the largest singular value of the Toeplitz matrix $\Tmat(f)$ as follows:

\begin{align}
\sigma_1 \left( \Tmat(f) \right) \leq \sup_{\omega \in [0, 2\pi]} |f(\omega)|.
\end{align}

<Ref> is a direct application of Lemma 4.1 in <cit.> for real Toeplitz matrices.

Let $F: \Rbb \rightarrow \Cbb^{m \times m}$ be a continuous and $2 \pi$-periodic matrix-valued function, then, $\Tmat(F) \in \Rbb^{mn \times mn}$ is a block Toeplitz matrix generated by the function $F$.
We can bound the largest singular value of the block Toeplitz matrix $\Tmat(F)$ as follows:

\begin{align}
\sigma_1 \left( \Tmat(F) \right) \leq \sup_{\omega \in [0, 2\pi]} \sigma_1 \left( F\left( \omega \right) \right) .
\end{align}

§.§ Upper-Bound on the Largest Singular Value of Doubly-Block Toeplitz Matrices

We extend the reasoning from Toeplitz and block Toeplitz matrices to doubly-block Toeplitz matrices (block Toeplitz matrices where each block is also a Toeplitz matrix).
A doubly-block Toeplitz matrix can be generated by a function $f: \Rbb^2 \rightarrow \Cbb$ using the 2-dimensional inverse Fourier transform.
For this purpose, we define an operator $\Dmat$ which maps a function $f: \Rbb^2 \rightarrow \Cbb$ to a doubly-block Toeplitz matrix of size $nm \times nm$.
For the sake of clarity, the dependence of $\Dmat(f)$ on $m$ and $n$ is omitted.
Let $\Dmat(f) \triangleq \leftmat \Dmat_{i,j}(f)\rightmat_{i,j \in \Iset^+_n}$ where $\Dmat_{i,j}(f)$ is a $m \times m$ matrix defined as:

\begin{equation} \label{equation:doubly_block_toeplitz_operator}
\Dmat_{i,j}(f) \triangleq \leftmat \frac{1}{4\pi^{2}} \int_{0}^{2\pi} \int_{0}^{2\pi} e^{- \ci \left((i-j)\omega_{1}+(k-l)\omega_{2}\right)}f(\omega_{1},\omega_{2}) \,\diff \omega_{1} \,\diff \omega_{2} \rightmat_{k,l \in \Iset^+_m} \enspace.
\end{equation}

We are now able to combine <Ref> to bound the largest singular value of doubly-block Toeplitz matrices with respect to their generating functions.
Note that in the following, we only consider generating functions as trigonometric polynomials with real coefficients therefore the matrices generated by $\Dmat(f)$ are real.

[Bound on the largest singular value of a Doubly-Block Toeplitz Matrix]
Let $f: \Rbb^2 \rightarrow \Cbb$ be a multivariate trigonometric polynomial of the form:

\begin{equation}\label{equation:muli_variate_poly_on_M}
f(\omega_1, \omega_2) \triangleq \sum_{h_1 \in \Iset_n} \sum_{h_2 \in \Iset_m} d_{h_1, h_2} e^{\ci (h_1 \omega_1 + h_2 \omega_2)}.
\end{equation}

Then, $\Dmat(f) \in \Rbb^{nm \times nm}$ is a doubly-block Toeplitz matrix where $d_{h_{1},h_{2}}$ is the ${h_2}^\textrm{th}$ scalar of the ${h_1}^\textrm{th}$ block of the matrix.
We can bound the largest singular value of the matrix $\Dmat(f)$ as follows:

\begin{align}
\sigma_{1} \left( \Dmat(f) \right) &\leq \sup_{\omega_1, \omega_2 \in [0, 2\pi]^2}|f(\omega_1,\omega_2)|
\end{align}

By definition, a doubly-block Toeplitz matrix is a block matrix where each block is a Toeplitz matrix.
Let $\Amat$ be a $mn \times mn$ doubly-block Toeplitz matrices determined by the sequence of blocks $\{\Amatsf^{(-n+1)}, \dots, \Amatsf^{(n-1)} \}$ of size $m \times m$ where the blocks $\Amatsf$ are Toeplitz matrices such that $\Amatsf^{(i)}$ is determined by the sequence of scalars $\{d_{i, -m+1}, \dots, d_{i, m-1} \}$.
Therefore the matrix $\Amat$ can be expressed with the operator $\Tmat$ with the matrix-valued generating function $F: \Rbb \rightarrow \Cbb^{n \times n}$ such that:

\begin{equation}
F(\omega_1) = \sum_{h_1 \in \Iset_n} \Amatsf^{(h_1)} e^{\ci h_1\omega_1}
\end{equation}

From <Ref> we have:

\begin{equation} \label{equation:th_bound_block_toeplitz}
\sigma_1\left(\Tmat(F) \right) \leq \sup_{\omega_1 \in [0,2\pi] } \sigma_{1}\left( F(\omega_1) \right)
\end{equation}

Note that because the function $F$ is a linear combination of the Toeplitz matrices $\Amatsf$ and that Toeplitz matrices are closed under addition and scalar product, $F(\omega_1)$ is also a Toeplitz matrix of size $m \times m$.
Therefore, we can define a function $f: \Rbb^2 \rightarrow \Cbb$ such that:

\begin{align}
f(\omega_1,\omega_2) &= \sum_{h_2 \in \Iset_m} \left[ F(\omega_1) \right]_{h_2} e^{\ci h_2 \omega_2} \\
f(\omega_1,\omega_2) &= \sum_{h_2 \in \Iset_m} \left[ \sum_{h_1 \in \Iset_n} \Amatsf^{(h_1)} e^{\ci h_1\omega_1} \right]_{h_2} e^{\ci h_2 \omega_2} \\
f(\omega_1,\omega_2) &= \sum_{h_2 \in \Iset_m} \left[ \sum_{h_1 \in \Iset_n} \Amatsf^{(h_1)} \right]_{h_2} e^{\ci \left( h_1\omega_1 + h_2 \omega_2 \right)} \\
f(\omega_1,\omega_2) &= \sum_{h_1 \in \Iset_n} \sum_{h_2 \in \Iset_m} d_{h_1, h_2} e^{\ci \left( h_1\omega_1 + h_2 \omega_2\right)} \enspace,
\end{align}

From <Ref>, we can write:

\begin{align}
\sigma_{1}\left( F(\omega_1) \right) &\leq \sup_{\omega_2 \in [0,2\pi]} \left| f(\omega_1, \omega_2) \right| \\
\Rightarrow \sup_{\omega_1 \in [0,2\pi]} \sigma_{1}\left( F(\omega_1) \right) &\leq  \sup_{\omega_1, \omega_2 \in [0,2\pi]^2} \left| f(\omega_1, \omega_2) \right| \\
\Rightarrow \sigma_1\left(\Tmat(F) \right) &\leq \sup_{\omega_1, \omega_2 \in [0,2\pi]^2} \left| f(\omega_1, \omega_2) \right|
\end{align}

Because the function $f(\omega_1,\ \cdot\ )$ is the generating function of $F(\omega_1)$ it is easy to show that the function $f$ is also the generating function of the matrix $\Tmat(F)$. Therefore, $\Tmat(F) = \Dmat(f)$ which concludes the proof.

§ EXTENDING THE BOUND TO CONVOLUTIONAL LAYERS

From now on, without loss of generality, we will assume that $n=m$ to simplify notations.
A discrete convolution operation with a 2-dimensional kernel applied on a 2-dimensional signal (, an image) is equivalent to a matrix multiplication with a doubly-block Toeplitz matrix <cit.>.
In practice, the input signal often has 3 or more dimensions called channels (for example, RGB images have 3 channels, one for each color).

If we denote $\cin$, the number of channels of the input signal, then, the input signal is a tensor of size $\cin \times n \times n$.
Moreover, when we perform multiple convolutions on the same signal the output signal will have multiple channels denoted $\cout$.
Therefore, the kernel is defined as a 4-dimensional tensor of size: $\cout \times \cin \times s \times s$.
The operation performed by a 4-dimensional kernel on a 3-dimensional signal can be formulated as the concatenation (horizontally and vertically) of doubly-block Toeplitz matrices.
Hereafter, we bound the singular value of multiple vertically stacked doubly-block Toeplitz matrices which corresponds to the operation performed by a 3-dimensional kernel with $\cout = 1$ on a 3-dimensional signal.

Consider doubly-block Toeplitz matrices $\Dmat(f_1), \dots, \Dmat(f_{\cin})$ where each $f_i: \Rbb^2 \rightarrow \Cbb$ is a multivariate polynomial of the same form as <Ref>.
Construct a matrix $\Mmat$ with $\cin\times n^2$ rows and $n^2$ columns, as follows:

\begin{equation}
\Mmat \triangleq \leftmat \Dmat^\top(f_1), \dots, \Dmat^\top(f_{\cin}) \rightmat^\top .
\end{equation}

Then, we can bound the largest singular value of the matrix $\Mmat$ as follows:

\begin{align} \label{equation:bound_asymptotic_equiv}
\sigma_1 \left( \Mmat \right) &\leq \sup_{\omega_1, \omega_2 \in \left[0, 2\pi\right]^2} \sqrt{ \sum_{i=1}^{\cin} \left|f_i\right (\omega_1, \omega_2)|^2} \enspace.
\end{align}

In order to prove <Ref>, we will need the following lemmas:

Let $f:\Rbb^2 \rightarrow \Cbb$ and $g:\Rbb^2 \rightarrow \Cbb$ be two continuous and $2\pi$-periodic functions. Let $\Dmat(f)$ and $\Dmat(g)$ be doubly-block Toeplitz matrices generated by the functions $f$ and $g$ respectively.

* $\Dmat^\top(f) = \Dmat(f^*)$

* $\Dmat(f) + \Dmat(g) = \Dmat(f + g)$

If the doubly-block Toeplitz matrix $\Dmat(f)$ is generated by a function $f: \Rbb^2 \rightarrow \Rbb$, then the matrix $\Dmat(f)$ is Hermitian.

If the doubly-block Toeplitz matrix $\Dmat(f)$ is generated by a non-negative function $f$ not identically zero, then the matrix $\Dmat(f)$ is positive definite.

Let $\Amat$ and $\Bmat$ be Hermitian positive semi-definite matrices. If $\Amat - \Bmat$ is positive semi-definite, then:

\begin{equation}
\lambda_1 \left( \Bmat \right) \leq \lambda_1 \left( \Amat \right)
\end{equation}

We need now to extend the well known Widom identity <cit.> which expresses the relation between Toeplitz and Hankel matrices to doubly-block Toeplitz and Hankel matrices.
Let us first generalize the doubly-block Toeplitz operator presented in <Ref>.

Given a function $f:\Rbb^2 \rightarrow \Cbb$, let $\Gmat^{\alpha_p} (f)$ be a matrix such that $\Gmat^{\alpha_p} (f) = \leftmat \Gmat^{\alpha_p}_{i,j}(f) \rightmat_{i,j \in \Iset^+_n}$ where $\Gmat^{\alpha_p}_{i,j}$ is defined as:

\begin{equation}
\Gmat^{\alpha_p}_{i,j}(f) =\leftmat \frac{1}{4\pi^{2}} \int_{0}^{2\pi} \int_{0}^{2\pi} e^{-\ci \alpha_p(i, j, k, l, \omega_1, \omega_2)} f(\omega_{1},\omega_{2}) \,\diff \omega_{1} \,\diff \omega_{2}
\rightmat_{k,l \in \Iset^+_n} \enspace.
\end{equation}

Note that as with the operator $\Dmat(f)$ we only consider generating functions as trigonometric polynomials with real coefficients therefore the matrices generated by $\Gmat(f)$ are real.
And as with the operator $\Dmat(f)$, the matrices generated by the operator $\Gmat^{\alpha_p}$ are of size $n^2 \times n^2$.

We will use the following $\alpha$ functions:

$\alpha_0(i, j, k, l, \omega_1, \omega_2) = (-j-i-1)\omega_1 + (k-l)\omega_2$

$\alpha_1(i, j, k, l, \omega_1, \omega_2) = (i-j)\omega_1 + (-l-k-1)\omega_2$

$\alpha_2(i, j, k, l, \omega_1, \omega_2) = (-j-i-1)\omega_1 + (-l-k-1)\omega_2$

$\alpha_3(i, j, k, l, \omega_1, \omega_2) = (-j-i+n)\omega_1 + (-l-k-1)\omega_2$

We now present the generalization of the Widom identity for Doubly-Block Toeplitz matrices below:

Let $f:\Rbb^2 \rightarrow \Cbb$ and $g:\Rbb^2 \rightarrow \Cbb$ be two continuous and $2\pi$-periodic functions.
Let $fg$ be the product of the functions $f$ and $g$ such that $(fg)(\omega_1, \omega_2) = f(\omega_1, \omega_2) g(\omega_1, \omega_2)$.
We can decompose the Doubly-Block Toeplitz matrix $\Dmat(fg)$ as follows:

\begin{equation}
\Dmat(fg) = \Dmat(f)\Dmat(g) + \sum_{p=0}^3 \Gmat^{\alpha_p \top}(f^*) \Gmat^{\alpha_p}(g) + \Jmat_{n^2} \left( \sum_{p=0}^3 \Gmat^{\alpha_p \top}(f) \Gmat^{\alpha_p }(g^*) \right) \Jmat_{n^2}.
\end{equation}

where $\Jmat_{n^2}$ is the reflection of the identity matrix of size $n^2 \times n^2$.

The proof of this Lemma is delayed to Appendix <ref>.
Now we have all the elements to prove <Ref> which bounds the largest singular value of vertically stacked doubly-block Toeplitz matrices with their generating functions.

Consider doubly-block Toeplitz matrices $\Dmat(f_1), \dots, \Dmat(f_{\cin})$ where each $f_i: \Rbb^2 \rightarrow \Cbb$ is a multivariate polynomial of the same form as <Ref>.
Construct a matrix $\Mmat$ with $\cin\times n^2$ rows and $n^2$ columns, such that:

\begin{equation}
\Mmat \triangleq \leftmat \Dmat^\top(f_1), \dots, \Dmat^\top(f_{\cin}) \rightmat^\top .
\end{equation}

First, let us observe the following equality which relates the largest singular value of the matrix $\Mmat$ and the largest eigenvalue of the sum of the doubly-block Toeplitz matrices composing $\Mmat$:

\begin{equation}
\sigma_1^2 \left( \Mmat \right) = \lambda_1 \left( \Mmat^{\top} \Mmat \right) = \lambda_1 \left( \sum_{i=1}^{\cin} \Dmat^{\top} \left(f_i \right) \Dmat (f_i) \right).
\end{equation}

Secondly, let us bound the largest eigenvalue of the sum of doubly-block Toeplitz generated by $|f_i|^2$:

\begin{align}
\lambda_1 \left( \sum_{i=1}^{\cin} \Dmat \left(|f_i|^2 \right) \right) \quad &= \quad \lambda_1 \left( \Dmat \left( \sum_{i=1}^{\cin} |f_i|^2 \right) \right) \\
\quad &= \quad \sigma_1 \left( \Dmat \left( \sum_{i=1}^{\cin} |f_i|^2 \right) \right) \\
\quad &\leq \quad \sup_{\omega_1, \omega_2 \in [0, 2\pi]^2} \sum_{i=1}^{\cin} |f_i(\omega_1, \omega_2)|^2.
\end{align}

where the first equality is due to <Ref>, the second equality is due to <Ref> and the last inequality is due to <Ref>.
To finalize the proof, we need to demonstrate that the following inequality holds:

\begin{equation} \label{equation:ch5-eq2}
\lambda_1 \left( \sum_{i=1}^{\cin} \Dmat^{\top} \left(f_i \right) \Dmat (f_i) \right) \leq \lambda_1 \left( \Dmat \left( \sum_{i=1}^{\cin} |f_i|^2 \right) \right).
\end{equation}

In order to prove the inequality above, let us study the positive definiteness of the following matrix:

\begin{equation}
\Dmat \left( \sum_{i=1}^{\cin} |f_i|^2 \right) - \sum_{i=1}^{\cin} \Dmat^{\top} \left(f_i \right) \Dmat (f_i),
\label{equation:ch5-eq3}
\end{equation}

One can observe that the term $\Dmat \left( \sum_{i=1}^{\cin} |f_i|^2 \right)$ of <Ref> is a real symmetric positive definite matrix by <Ref>.
Furthermore, the term $\sum_{i=1}^{\cin} \Dmat^{\top} \left(f_i \right) \Dmat (f_i)$ of <Ref> is a sum of positive semi-definite matrices.
Therefore, if the subtraction of the two is positive semi-definite, one could apply <Ref> to prove the <Ref>.
We know from <Ref> that

\begin{equation}
\Dmat(fg) - \Dmat(f)\Dmat(g) = \sum_{p=0}^3 \Gmat^{\alpha_p \top}(f^*) \Gmat^{\alpha_p}(g) + \Jmat \left( \sum_{p=0}^3 \Gmat^{\alpha_p \top}(f) \Gmat^{\alpha_p}(g^*) \right) \Jmat.
\end{equation}

By choosing $f = f^*$, $g = f$ and with the use of <Ref>, we obtain:

\begin{align} \label{equation:widom_identity_block_topelitz}
\Dmat(f^* f) - \Dmat(f^*)\Dmat(f)
&= \Dmat(|f|^2) - \Dmat^{\top}(f)\Dmat(f) \\
&= \sum_{p=0}^3 \Gmat^{\alpha_p \top}(f)\Gmat^{\alpha_p}(f) + \Jmat \left( \sum_{p=0}^3 \Gmat^{\alpha_p \top}(f^*)\Gmat^{\alpha_p} (f^*) \right) \Jmat .
\end{align}

From <Ref>, we can see that the matrix $\Dmat(|f|^2) - \Dmat^{\top}(f)\Dmat(f)$
is positive semi-definite because it can be decomposed into a sum of positive semi-definite matrices and because positive semi-definiteness is closed under addition, we have:

\begin{align}
\sum_{i=1}^{\cin} \leftmat \Dmat \left( |f_i|^2 \right) - \Dmat^{\top} \left(f_i \right) \Dmat (f_i) \rightmat &\geq 0
\end{align}

By re-arranging and with the use <Ref>, we obtain:

\begin{align}
\sum_{i=1}^{\cin} \Dmat \left( |f_i|^2 \right) - \sum_{i=1}^{\cin} \leftmat \Dmat^{\top} \left(f_i \right) \Dmat (f_i) \rightmat &\geq 0 \\
\Dmat \left( \sum_{i=1}^{\cin} |f_i|^2 \right) - \sum_{i=1}^{\cin} \leftmat \Dmat^{\top} \left(f_i \right) \Dmat (f_i) \rightmat &\geq 0
\end{align}

We can conclude that the <Ref> is true and therefore by <Ref> we have:

\begin{align}
\lambda_1 \left( \sum_{i=1}^{\cin} \Dmat^{\top} \left(f_i \right) \Dmat (f_i) \right) &\leq \lambda_1 \left( \Dmat \left( \sum_{i=1}^{\cin} |f_i|^2 \right) \right) \\
\sigma_1^2 \left( \Mmat \right) &\leq \sup_{\omega_1, \omega_2 \in [0, 2\pi]^2} \sum_{i=1}^{\cin} |f_i(\omega_1, \omega_2)|^2 \\
\sigma_1 \left( \Mmat \right) &\leq \sup_{\omega_1, \omega_2 \in [0, 2\pi]^2} \sqrt{ \sum_{i=1}^{\cin} |f_i(\omega_1, \omega_2)|^2 }
\end{align}

which concludes the proof.

To have a bound on the full convolution operation, we extend <Ref> to take into account the number of output channels.
The matrix of a full convolution operation is a block matrix where each block is a doubly-block Toeplitz matrix.
Below, we present our main result:

[Bound on the largest singular value on the discrete convolution operation]
Let us define doubly-block Toeplitz matrices $\Dmat(f_{1, 1}), \dots, \Dmat(f_{\cin, \cout})$ where each $f_{i,j}: \Rbb^2 \rightarrow \Cbb$ is a multivariate polynomial of the same form as <Ref>.
Construct a matrix $\Mmat$ with $\cin\times n^2$ rows and $\cout\times n^2$ columns.
We can bound the largest singular value of the matrix $\Mmat$ as follows:

\begin{equation} \label{equation:lipbound_sv}
\sigma_1 \left( \Mmat \right) \leq \sqrt{ \sum_{i=1}^{\cout} \sup_{\omega_1, \omega_2 \in [0, 2\pi]^2} \sum_{j = 1}^{\cin} \left|f_{ij}(\omega_1, \omega_2) \right|^2 } .
\end{equation}

First, in order to prove <Ref>, we will need the following lemma which bounds the singular values of a matrix constructed from the concatenation of multiple matrices.

Let us define matrices $\Amat^{(1)}, \dots, \Amat^{(p)}$ with $\Amat^{(i)} \in \Rbb^{n \times n}$. Let us construct the matrix $\Mmat \in \Rbb^{n \times pn}$ as follows:

\begin{equation}
\Mmat \triangleq \leftmat \Amat^{(1)}, \dots, \Amat^{(p)} \rightmat
\end{equation}

where $\leftmat\ \cdot\ \rightmat$ define the concatenation operation. Then, we can bound the singular values of the matrix $\Mmat$ as follows:

\begin{equation}
\sigma_1(\Mmat) \leq \sqrt{\sum_{i=1}^p \sigma_1(\Amat^{(i)})^2}
\end{equation}

\begin{align}
\sigma_1\left(\Mmat\right)^2 &= \lambda_1\left(\Mmat \Mmat^\top\right) \\
&= \lambda_1\left( \sum_{i=1}^p\Amat^{(i)} \Amat^{(i)\top}  \right) \\
&\leq \sum_{i=1}^p \lambda_1\left( \Amat^{(i)} \Amat^{(i)\top}  \right) \\
&\leq \sum_{i=1}^p \sigma_1\left( \Amat^{(i)} \right)^2 \\
\Leftrightarrow \sigma_1\left(\Mmat\right) &\leq \sqrt{\sum_{i=1}^p \sigma_1(\Amat^{(i)})^2}
\end{align}

which concludes the proof.

Now, the proof of <Ref> is a combination of <Ref> and <Ref>.

Let us define the matrix $\Mmat^{(i)}$ as follows:

\begin{equation}
\Mmat^{(i)} = \leftmat \Dmat(f_{1, i})^\top, \dots, \Dmat(f_{\cin, i})^\top \rightmat^\top .
\end{equation}

We can express the matrix $\Mmat$ as the concatenation of multiple $\Mmat^{(i)}$ matrices:

\begin{equation}
\Mmat = \leftmat \Mmat^{(1)}, \dots, \Mmat^{(\cout)} \rightmat
\end{equation}

Then, we can bound the singular values of the matrix $\Mmat$ as follows:

\begin{align}
\sigma_1\left(\Mmat\right) &\leq \sqrt{\sum_{i=1}^{\cout} \sigma_1(\Mmat^{(i)})^2} \\
\sigma_1\left(\Mmat\right) &\leq \sqrt{\sum_{j=1}^{\cout} \sup_{\omega_1, \omega_2 \in [0, 2\pi]^2} \sum_{i=1}^{\cin} |f_{i,j}(\omega_1, \omega_2)|^2 }
\end{align}

where the first inequality is due to <Ref> and the second one is due to <Ref>.
This concludes the proof.

<Ref> depends on the convolution matrix $\Mmat$, however, we can easily formulate the bound with the values of a 4-dimensional kernel.
Let us define a kernel $\Kmat \in \Rbb^{\cout \times \cin \times s \times s}$, a padding $p \in \Nbb$ and $d = \lfloor s / 2 \rfloor$ the degree of the trigonometric polynomial, then:

\begin{equation}
f_{ij}(\omega_1, \omega_2) = \sum_{h_1 = -d}^d \sum_{h_2 = -d}^d k_{i, j, h_1,h_2} e^{\ci (h_1 \omega_1 + h_2 \omega_2)}.
\end{equation}

where $k_{i, j, h_1,h_2} = \leftmat \Kmat \rightmat_{i, j, a, b}$ with $a =  s - p - 1 + h_1$ and $b =  s - p - 1 + h_2$.

In the rest of the chapter, we will refer to the bound in <Ref> applied to a kernel as $\lipbound$ and we denote $\lipbound(\Kmat)$ the Lipschitz upper bound of the convolution performed by the kernel $\Kmat$.

§ COMPUTATION AND PERFORMANCE ANALYSIS OF LIPBOUND

This section aims at analyzing the bound introduced in <Ref>.
First, we present an algorithm to efficiently compute the bound, we analyze its tightness by comparing it against the true largest singular value.
Finally, we compare the efficiency and the accuracy of our bound against the state-of-the-art methods.

§.§ The Maximum Modulus of a Trigonometric Polynomial

In order to compute $\lipbound$ from <Ref>, we have to compute the maximum modulus of several trigonometric polynomials.
However, finding the maximum modulus of a trigonometric polynomial has been known to be NP-Hard <cit.>, and in practice they exhibit low convexity (see <Ref>).
We found that for 2-dimensional kernels, a simple grid search algorithm such as PolyGrid (see <Ref>), works better than more sophisticated approximation algorithms ( <cit.>).
This is because the complexity of the computation depends on the degree of the polynomial which is equal to $\lfloor s / 2 \rfloor$ where $s$ is the size of the kernel and is usually small in most practical settings ($s=3$).
Furthermore, the grid search algorithm can be parallelized effectively on CPUs or GPUs and runs within less time than alternatives with lower asymptotic complexity.

kernel $1\times3\times3$

kernel $9\times3\times3$

kernel $1\times5\times5$

kernel $9\times5\times5$

Contour plots of multivariate trigonometric polynomials where the values of the coefficient are the values of a random convolutional kernel. The red dots in the figures represent the maximum modulus of the trigonometric polynomials.

[1]
PolyGrid$f, S$polynomial $f$, number of samples $S$
$\sigma \gets 0$, $\omega_1 \gets 0$, $\epsilon \gets \frac{2\pi}{S}$
$i=0$ to $S-1$
$\omega_2 \gets 0$
$j=0$ to $S-1$
$\omega_2 \gets \omega_2 + \epsilon$
$\sigma \gets \max( \sigma, |f(\omega_1, \omega_2)|)$
$\omega_1 \gets \omega_1 + \epsilon$
return $\sigma$ approximated maximum modulus of $f$

PolyGrid Algorithm

To fix the number of samples $S$ in the grid search, we rely on the work of <cit.>, who has analyzed the quality of the approximation depending on $S$.
Following this work we first define $\Theta_S$, the set of $S$ equidistant sampling points as follows:

\begin{equation}
\Theta_S \triangleq \left\{ \omega \mid \omega = k \cdot \frac{2\pi}{S} \mbox{ with }  k = 0, \dots, S-1 \right\}.
\end{equation}

Then, for a trigonometric polynomial $f: [0, 2\pi]^2 \rightarrow \Cbb$, we have:

\begin{equation}
\max_{\omega_1, \omega_2 \in [0,2\pi]^2} \left| f(\omega_1, \omega_2) \right| \leq (1 - \alpha)^{-1} \max_{\omega_1', \omega_2' \in \Theta_S^2} \left| f(\omega_1', \omega_2') \right|,
\end{equation}

where $d$ is the degree of the polynomial and $\alpha = 2d / S$.
For a $3\times3$ kernel which gives a trigonometric polynomial of degree 1, we use $S = 10$ which gives $\alpha = 0.2$.
Using this result, we can now compute $\lipbound$ for a convolution operator with $\cout$ output channels as per <Ref>.

§.§ Analysis of the Tightness of the Bound

In this section, we study the tightness of the bound with respect to the dimensions of the doubly-block Toeplitz matrices.
For each $n \in \Nbb$, we define the matrix  $\Mmat^{(n)}$ of size $kn^2 \times n^2$ as follows:

\begin{equation}
\Mmat^{(n)} \triangleq \textstyle \leftmat \Dmat^{(n)\top}(f_1), \dots, \Dmat^{(n)\top}(f_k) \textstyle \rightmat^\top
\end{equation}

where the matrices $\Dmat^{(n)}(f_i)$ are of size $n^2 \times n^2$.
To analyze the tightness of the bound, we define the function $\Gamma$, which computes the difference between $\lipbound$ and the largest singular value of the function $\Mmat^{(n)}$:

\begin{equation} \label{equation:function_convergence}
\Gamma(n) = \lipbound(\Kmat_{\Mmat^{(n)}}) - \sigma_1(\Mmat^{(n)})
\end{equation}

where $\Kmat_{\Mmat^{(n)}}$ is the convolution kernel associated with the matrix $\Mmat^{(n)}$.

Representation of the function $\Gamma(n)$ defined for different kernel size.

To compute a very close approximation of the exact largest singular value of $\Mmat^{(n)}$ for a specific $n$, we use the Implicitly Restarted Arnoldi Method (IRAM)  <cit.> available in SciPy.
The results of this experiment are presented in <Ref>.
We observe that the difference between the bound and the actual value (approximation gap) quickly decreases as the input size increases.
For an input size of $50$, the approximation gap is as low as $0.012$ using a standard $6\times3\times3$ convolution kernel.
For a larger input size such as ImageNet ($224$), the gap is lower than $4.10^{-4}$.
Therefore $\lipbound$ gives an almost exact value of the largest singular value of the operator matrix for most realistic settings.

§.§ Comparison of LipBound with State-of-the-Art Approaches

4c1x3x3           4c32x3x3

3-6 8-11
2cRatio     2cTime (ms)
2cRatio     2cTime (ms)

0.431     0.042     1088     251          0.666     0.123     1729     399

1.293     0.126     1.90     0.48         1.441     0.188     1.90     0.46

0.973     0.006     4.30     0.64         0.972     0.004     4.93     0.67

LipBound          0.992     0.012     0.49     0.05         0.984     0.021     0.63     0.46

Comparison of the accuracy of approximation methods for computing an approximation of the largest singular value of a convolution layer.

Network         2cLipBound (ms)         2cPower Method (ms)         Ratio

AlexNet         4.75     1.10         38.75     2.52         8.14

ResNet 18         29.88     1.73         148.35     14.92         4.96

ResNet 34         54.73     3.62         266.85     25.35         4.87

ResNet 50         60.77     4.62         467.61     36.52         7.69

ResNet 101         102.72     11.53         817.06     102.87         7.95

ResNet 152         158.80     20.84         1373.57     164.37         8.64

DenseNet 121         125.55     14.59          937.35      11.52         7.46

DenseNet 161         176.11     19.13         1292.61      30.50         7.33

DenseNet 169         188.03     19.74         1372.62      21.16         7.29

DenseNet 201         281.13     23.41         1930.19     170.79         6.86

VGG 11         13.73     1.19          81.78     4.45         5.95

VGG 13         14.96     1.99         102.04     4.20         6.82

VGG 16         21.92     1.94         132.29     5.99         6.03

VGG 19         29.05     0.66         162.28     4.87         5.58

WideResNet 50-2         113.28     45.44         468.74     6.54         4.13

SqueezeNet 1-0         18.44     5.93         222.40     25.49         12.05

SqueezeNet 1-1         18.26     6.65         209.80      3.59         11.48

Efficiency of LipBound computation the Power Method with 10 iterations on full networks.

In this section we compare our PolyGrid algorithm with the values obtained using alternative approaches.
We consider the 3 alternative techniques by <cit.> which have been described in <Ref>, <Ref>.

To compare the different approaches, we extracted 20 kernels from a trained model.
For each kernel we construct the corresponding doubly-block Toeplitz matrix and compute its largest singular value.
Then, we compute the ratio between the approximation obtained with the considered approach and the approximated singular value obtained by IRAM, and average the ratios over the 20 kernels.
Thus good approximations result in approximation ratios that are close to 1.
The results of this experiment are presented in <Ref>.
The comparison has been made on a Tesla V100 GPU.
The time was computed with the PyTorch CUDA profiler and we “warmed” up the GPU before starting the timer for caching purposes.

The method introduced by <cit.> and presented in <Ref> computes the largest singular value of convolution layers based on doubly-block circulant matrices.
Doubly-block circulant matrices perform a convolution with a “wrapping around” operation which do not correspond to the more general setting.
We can see in <Ref> that the values differ by an important margin.
This technique is also computationally expensive as it requires computing the SVD of $n^2$ small matrices where $n$ is the size of inputs.
<cit.> have shown that the singular value of the reshape kernel is a bound on the largest singular value of the convolution layer.
Their approach is very efficient but the approximation is loose and overestimate the real value.
As said previously, the power method provides a good approximation at the expense of the efficiency.
We also compare our approach to the power method with 10 iterations from  <cit.> (see <Ref>).
The results show that our proposed technique: PolyGrid algorithm can get the best of both worlds.
It achieves a near perfect accuracy while being very efficient to compute.

The results of <Ref> shows the performance for the computation for only one convolution layer.
However, during the training Lipbound or the power method need to be computed for every layer of the network and the computation time is dependent on the architecture of the network, for example, the size of the activations or the size of the kernels.
In  <Ref>, we compare our approach method against the power method on the full architecture, , the time needed for the computation on all the layers of the networks.
We measure on the following convolutional architectures: AlexNet <cit.>, ResNet <cit.>, DenseNet <cit.>, VGG <cit.>, WideResNet <cit.>, SqueezeNet <cit.>.
<Ref> shows that our approach is systematically faster than the power method by a factor up to 12 when considering all the layers of the networks.
This demonstrates the scalability of our method.

§ LIPSCHITZ REGULARIZATION FOR ADVERSARIAL ROBUSTNESS

One promising application of Lipschitz regularization is in the area of adversarial robustness.
Empirical techniques to improve robustness against adversarial examples such as Adversarial Training only impact the training data,  and often show poor generalization capabilities <cit.>.
<cit.> have shown that the adversarial generalization error depends on the Lipschitz constant of the network, which suggests that the adversarial test error can be improved by applying Lipschitz regularization in addition to adversarial training.

In this section, we illustrate the usefulness of LipBound by training a Wide ResNet <cit.> with Lipschitz regularization and adversarial training.
Our regularization scheme is inspired by the one used by <cit.> but instead of using the power method, we use our PolyGrid algorithm presented in <Ref> which efficiently computes an upper bound on the largest singular value of convolution layers.

We introduce the AT+LipReg loss to combine Adversarial Training and our Lipschitz regularization scheme in which layers with a large Lipschitz constant are penalized.
We consider a neural network $N_\Omega : \Xset \rightarrow \Yset$ with $\depth$ layers $\phi^{(1)}_{\Wmat^{(1)}, \bvec^{(1)}}, \dots, \phi^{(\depth)}_{\Wmat^{(\depth)}, \bvec^{(\depth)}}$ where $\Wmat^{(1)}, \dots, \Wmat^{(\depth)}$ are the weight matrices and $\Omega$ is the union of all the parameters as defined in <Ref>.
Given a distribution $\Dset$ over $\Xset \times \Yset$, we can train the parameters of the network by minimizing the AT+LipReg loss as follows:

\begin{equation} \label{equation:ch5-obj_function}
\min_{\Omega} \Ebb_{\xvec,y \sim \Dset} \left[ L(N_\Omega(\xvec + \adv^{\text{adv}}_\Omega(\xvec)), y) + r(\Omega) \right]
\end{equation}

where $L$ is the cross-entropy loss function, $\adv^{\text{adv}}_\Omega(\xvec)$ is an adversarial perturbation following the loss maximization strategy presented in <Ref> and the regularization function $r$ is defined as follows:

\begin{equation}
r(\Omega) =  C_1 \underbrace{\sum_{(\Wmat, \bvec) \in \Omega} \left( \norm{\Wmat}_\fro + \norm{\bvec}_2 \right)}_{
\text{$\ell_2$ regularization}} + C_2 \underbrace{\vphantom{\sum_{(\Wmat, \bvec) \in \Omega}} \sum_{i=1}^{\depth-1} \log\left(\lipbound\left(\Kmat_{\Wmat^{(i)}}\right)\right)}_{\text{Lipschitz regularization}}
\end{equation}

where $C_1$, $C_2$ are two user-defined hyper-parameters.
Note that regularizing the sum of logs is equivalent to regularizing the product of all the $\lipbound$ which is an upper bound on the global Lipschitz constant.
In practice, we also include the upper bound on the Lipschitz of the batch normalization since we can compute it very efficiently (see Appendix C.4.1 of  <cit.>) but we omit the last fully connected layer.

In this section, we compare the robustness of Adversarial Training <cit.> against the combination of Adversarial Training and Lipschitz regularization.
To regularize the Lipschitz constant of the network, we use the objective function defined in Equation <ref>.
We train Lipschitz regularized neural networks with LipBound (see <Ref>) implemented with PolyGrid (see <Ref>) (AT+LipBound) with $S = 10$ or with the specific power method for convolutions introduced by <cit.> with 10 iterations (AT+PM).

Table <ref> shows the gain in robustness against strong adversarial attacks across different datasets.
We can observe that both AT+LipBound and AT+PM offer a better defense against adversarial attacks and that AT+LipBound offers a further improvement over the Power Method.
<Ref> show the Accuracy under attack with different numbers of iterations of the PGD algorithm.
Table <ref> presents our results on the ImageNet Dataset.
First, we can observe that the AT+LipReg trained networks offer a better generalization than with standalone Adversarial Training.
Secondly, we can observe the gain in robustness against strong adversarial attacks.
Network trained with Lipschitz regularization and Adversarial Training offer a consistent increase in robustness across $\ell_\infty$ and $\ell_2$ attacks with different $\epsilon$ value.