Nicolas Barroso

Amélie Fondevilla

David Vanderhaeghe

/ R.Vergne and N.Mellado. (Paper Chairs)

Keywords:

Figure 1: From a user provided set of key frames (designated by arrows) and motion field, our approach generates intermediate to create an animation.

Introduction

2D animated feature film exhibits a wide range of visual appearances. In this paper we focus on appearance that shows individual 2D marks as in sketching or painting. In this kind of style we distinguish two main workflows. The classical workflow is to create key † E-mail: nicolas.barroso@irit.fr frames to define the movement, and then create the necessary inbetweening to smooth this movement. The other workflow, known as "under the camera", is to create a first frame, then slightly modify this frame to obtain the next one, and so on. Each frame of the animation only exists when its captured and is then destroyed to produce the next one. The under the camera workflow needs specific artistic skills to anticipate the whole animation sequence. Computers algorithm can help artists in many tedious animation tasks, improving the efficiency of the pipeline in various ways [GM19]. Our long term goal is to provide computer tools to have the look and feel of the under the camera workflow with the level of control of the classical workflow. This paper focus on the automatic creation of intermediate frames given a set of key frames. We target professional animators and need to provide fine control over each frame, even those computed automatically. To this end we use stroke based rendering [VC12] where each stroke of the animation is represented independently, rather than keeping only the pixels representation of the frames. Strokes are rendered using a paint simulator to compute the final frame appearance.

The main contributions of the paper are:

• the extrapolation of 3D rendered motion field to obtain a 2D motion field suitable for stroke propagation (Sec 4), • an optimization scheme to generate guide strokes from key frames and motion field (Sec 5), • a frame generator that leverage guide stroke from up to two key frames to produce intermediate frames (Sec. 6). These contributions allows us to elaborate a first inbetween pipeline to generate intermediates frames in a stroke based rendering framework.

Related Works

Synthesis of Intermediate frames. The automatic synthesis of intermediate frames, also called inbetweening, is a well-known domain in the field of 2D animation. Most methods work in image space, by synthesizing texture depicting as-rigid-as-possible deformation [SDC09,BCK * 13,DBB * 17], or using a triangulation of the image space and applying deformation directly on this triangulated mesh, using 2D skinning techniques [BKLP16]. Texture synthesis approaches are agnostic to the tools used to produce the example. The resulting frames are raster images. Our approach generates strokes that give high level structure for post editing. Closer to our method, other approaches are working at the stroke level, deducing the deformation between two frames through stroke pairing [WNS * 10], or using directly the motion of an underlying 3D animation [WDK * 12]. While most of the stroke-based approaches are focused on sketches, a few of them deal with paintings. In this case, the inbetweening problem can be overcome by embedding the painting strokes in a 3D space [START_REF] Schmid | Overcoat: An implicit canvas for 3d painting[END_REF], where the strokes can follow the motion of a 3D surface [BBS * 13]. We propose to define all the strokes in the 2D image which corresponds to the traditional animation pipeline.

Temporal coherence. Temporal coherence plays a significant role in the perception of stylized animations. In texture approach it can be ensure by energy minimization [BCK * 13] or neural network architectures [START_REF] Ruder | Artistic style transfer for videos and spherical images[END_REF]. To get closer to the hand drawn look some works try to reduce the temporal coherence of the stylization process and inject randomness through noise control in their algorithm [FLJ * 14] [FJS * 17] or use marks approach to obtain a more finer control of the randomness [CZBB20]. As chown by Delanoy et. al. [DBH19], processing the motion rather than the image, with motion rigidifaction for instance, break temporal coherence provided by the initial optical flow and results in 2D look closer to manual methods. Our approach focus on automatic, yet plausible, intermediate frame generation, with the noise inherent to hand made animation, rather than enforcing temporal coherence of the animation.

Generation of strokes by example. Chen et. al. [CZBB20]

builds strokes of intermediate frames by rigid transformation of the strokes of key frames. Intermediate frames are composed by selecting and transforming strokes to reach a target density map. The density map is estimated by interpolating key frame's densities with the image registration method [SDC09]. Combined with transformed strokes the animation sequence shows an hand drawn look. We share the same approach, but in our case case, the strokes' transformation are guided by an underlying motion field that gives the expected motion. In addition, our approach modifies the curve shape using optimization scheme that preserve geometric properties from the key frames.

Overview

Key frames refer to frame drawn, by the artist, and generated frames are automatically rendered by our approach. The animation workflow is as follow: the artist draw one or more key frames which convey the example style and appearance to reproduce (Fig. 1). The artist also provides a motion field that encompass the animation motion. Our approach generates intermediate frames taking key frames and motion, one frame after the other (Fig. 2). At any time, the artist can decide that frame do not corresponds to his wishes and modify or redraw it, this introduce a new key frame in the animation. Intermediate frames are updated to take into account this new key frame. For the rest of the paper we consider that intermediate frames are surrounded by two key frames, one before and one after along the timeline.

We focus on the animation of one object, multiple objects and background are handled using layers. The final animation is composed of multiple layers, each layer is animated and rendered independently, using a different motion field and key frames. Layers composition can then be mixed as painted on different transparent canvas, or mixed as if they are paint one over the other. Key frames and intermediate frames are defined by an ordered list of strokes. A stroke is defined by a curve represented as a polyline, and a set of parameters, i.e. the quantity and color of paint on the virtual tool and the pressure of the tool on the canvas along the curve. The final image is obtained by rendering the strokes with a paint simulator. There is two main steps to generate an intermediate frame: the propagation of guide strokes (Sec. 5) using motion field (Sec. 4) and the generation of frame's strokes (Sec. 6).

Motion Field

Our method needs a motion field for each frame of the animation as input. This motion field can come from different sources, for instance it can be handcrafted by the artist. In practice, for the example shown in the paper, we start from a simple 3D animation capturing the motion to convey. Since the motion field is created from the render of a 3D scene, we propose to display the 3D scene as a background template to draw a key frame. The proposed UI allows the artist to choose whether to draw or not over the underlying image.

For each frame, the motion field contains two motion vectors for each pixel of the image, one vector to next frame position of this pixel and one vector to previous frame position. Since the motion field is null outside of the animated object's surface, we need to extend the motion information for each background pixels to let the artist paint over the background as well. To this end we compute biharmonic weights as describe by Baster et. al. [START_REF] Baxter | N-way morphing for 2d animation[END_REF] The main step of the algorithm is to build a triangular 2D mesh over the image plane to finally interpolate vertices motion vectors for each pixels. Vertices of this mesh are evenly distributed over the image plane. Each vertex over a surface is assigned the underlying motion vector and become a seed point to the interpolation. The motion vector of vertices falling on background pixels are computed using biharmonic weight computed for each vertices according to the seeds, as shown Figure 3. Finally, the motion vector of each pixels of the image plane is computed using barycentric interpolation from the motion field of the triangular mesh.

Moreover we store a confidence value with each motion vector. The confidence is computed by the same barycentric interpolation considering a value of one for each vertices of the 2D mesh that has valid motion and zero for null vectors.

Guide Strokes

To guide the generation of an intermediate frame we generate a set of guide strokes. We compute a guide stroke for each strokes of a key frame. A guide stroke captures both key frame content, i.e. the curve, color, and pressure variation, and the motion information, i.e. the motion fields from stroke time to intermediate frame time. Let's consider one stroke of one key frame at time t 0 , as stated before, the stroke's curve is described by a poly line, we call the curve of the stroke in the key frame the reference curve. The advected curve is computed by moving each poly-line vertex according to underlying motion field. Advection only is not sufficient, it can break the curve shape and also spread curve away from the underlying object due to the motion field interpolation. We regularize the stroke curve using an optimization.

We leverage the confidence value to optimize the guide stroke curve. The optimization process mix the reference curve, with the advected curve according to confidence value: closer to reference curve when the confidence is low ; or closer to the advected curve when confidence is high.

We use the following three energy functions. These energies are measured along the curve for N evenly distributed sample points. The first constraint ensures the optimization will converge spatially close to the advected curve.

Es = N-1 ∑ i=0 c adv i N -copt i N 2
where c(p) is the evaluation of a 2D point of the parametric curve function at x ∈ [0, 1], c adv is the advected stroke and copt is the curve of the guide stroke.

The two next energy functions are weighted by the confidence σ(x) along the curve of the advected stroke. The second constraint penalizes curvature differences between the reference curve and the curve of the guide stroke:

Ec = N-1 ∑ i=0 γ re f i N -γopt i N 2 × 1 -σ i N
where γ is the curvature computed with finite differences.

The last energy ensures that the length between two consecutive samples of the curve of the guide stroke remains close to the length of corresponding samples in the reference curve.

δ(a, b) = ∥ f (a) -f (b)∥ E l = N-1 ∑ i=1 δ re f i N , i -1 N -δopt i N , i -1 N 2 × 1 -σ i N
Since the stroke poly-line have a lot of points, the optimization would be underconstrained. So, we first fit the poly-line by piecewise cubic Bézier curves, resulting in fewer variables to optimize. We use the Levenberg Marquardt algorithm to solve the optimization. The variables are control points of Bézier curves, and the energy minimized is

E = Es + Ec + E l
The curve regularization is shown Figure 4, where reference curve, advected curve and guide curve are shown for a complex motion. This process is repeated for each stroke of the key frame. And for each intermediate frame, computing the new advected curve using the previous guide stroke as a starting position. For each intermediate frame we also compute guide strokes from the next key frame, using the same process using backward motion vector.

Stroke Generation

The stroke generation stage consists in creating the set of stroke of an intermediate frame given the guide strokes. To do so, we make three straightforward assumptions:

1. The more an intermediate frame is close in time to key frame, the more it should look like this key frame. 2. An intermediate frame should reflect the content of the key frames before and after it. 3. The intermediate frame should have been done by hand, or at least look likes a hand made drawing.

We propose to generate the strokes of an intermediate frame as follow: We select a subset of the guide strokes from the two key frames.

To go from key frame A at time ta to key frame B at time t b = ta + N we progressively select guide strokes generated from key frame B while deselecting stroke generated from key frame A.

To ensure a complete representation of the animated object, we choose to have at least the guide strokes from one of the key frame fully selected at each intermediate frame. To this end we define a selection ratio for intermediate frame at time ta +i as min 1, 2(N-i) N

Selected guide strokes generated from key frame A is 100% for intermediate frames from time ta to ta + N/2, and decrease to 0% at time t b . The ratio of guide stroke from key frame B is computed similarly, but reversed in time. We design this selection scheme to ensure correct coverage in the intermediate frames.

We tested two strategies to select the guide stroke given the ratio. The first strategy select a random subset of the guide strokes according to the ratio independently for each intermediate frame.

The second strategy randomize the list of guide stroke once for all, then selects the strokes according to the ratio following the list order. The artist can test both strategy and pick the one he prefer.

When rendering the intermediate frame, the stroke drawing order follow the order in the key frames: We assign a scalar value as draw time between 0 and 1 to each of the strokes of a key frame, according to the drawing order of the artist. The selected strokes of both key frames are sorted by this draw time before rendering.

Implementation and Results

Our C++/OpenGL prototype use Radium Engine [MRB * 21] as main rendering engine. We use our own implementation of the paint simulator presented by Baxter et. al. [CBWG10]. This paint simulator compute bi-directional paint exchanges between the brush and the canvas, on the GPU. Strokes advection is done on the GPU using compute shader. The Levemberg Marquardt optimization uses Eigen [GJ * 10] implementation. We extract motion fields through Blender AOV rendering [Com18], we use triangle lib in python [She96, R * 20] and our implementation of bi-harmonic weight interpolation to obtain the interpolated motion field, as a pre-process. The time needed to render a frame depends on the number of guide strokes generated and the number of pixel covered by the strokes. In practice, the examples shown in the paper and supplemental videos took less than one second per frame to render, including advection, optimization and paint simulation. These render time are suitable for interaction.

Conclusion and Discussion

Our method generate intermediate frames in a stroke painterly style and produce an animation with a traditional hand drawn look and feel. It takes as input key frames drawn by an artist, and a motion field to describe how the strokes should move. While we diverge from the traditional paint under the camera look and feel, we ease the creation of animation in a stroke based rendering context by automatically creating intermediate frames.

The main limitation of the approach concern the complexity of the motion we can depict. For instance, when two moving objects cross each others on the same layer, some point of the poly-line of stroke of one object will follow the other object motion. To solve this issue, we think a better registration of each stroke with the underlying motion field could be envisioned.

Several improvements are left as future work, for instance during the generation of the strokes of an intermediate frame, we can imagine to add some randomness in the strokes parameters (curve shape, pressure, brush color). In our current selection scheme, the intermediate frame at mid-time between two key frames render the full set of strokes from the two key frames. The resulting image depict thicker paint on the canvas, which could be perceived as an artifact. Also, the generation of strokes from guides strokes could have a target coverage over the image plane to remove this artifact.

We think that the confidence value could convey artistic expression, so that some strokes keep their original shapes during advection, while others better follow the underlying motion. To this end we plan to propose confidence brushes to alter the confidence of the motion field, either for some regions, or for specific selected strokes.

The creation of an animated 3D scene to compute the motion field adds a large amount of work. We see two diverging workflows here: Either stay in pure 2D animation, and explore new tools to build the motion field from artist gesture, or take more benefit from the 3D animation, such as color, surface normal, and lighting, to make working time profitable.

Figure 2 :

 2 Figure 2: We generate an intermediate frame, here at time 12, using the guide strokes computed from the key frames before (at time 10) and after (at time 15) on the timeline.

Figure 3 :

 3 Figure 3: Computed motion field of a bouncing ball. Left: motion field as computed from the 3D renderer. Middle: 2D mesh and known vectors in green. Right: computed motion vector (in red) from the known vectors using bi-harmonic interpolation.

Figure 4 :

 4 Figure 4: (a) Example of two reference strokes on a key frame painted in yellow, (b) advected strokes deformed by advection only (c) or modified by our optimization, here with confidence set to 0 everywhere for illustration purpose.

© 2021 The Author(s) Journées Françaises de l'Informatique Graphique 2021 .

Acknowledgement This work has been partially supported by the Structures JCJC ANR project ANR-19-CE38-0009-01.


Peiyuan Zhang

Antonio Orvieto
email: aorvieto@ethz.ch

Hadi Daneshmand
email: seyed.daneshmand@inria.fr

Rethinking the Variational Interpretation of Accelerated Optimization Methods

The continuous-time model of Nesterov's momentum provides a thought-provoking perspective for understanding the nature of the acceleration phenomenon in convex optimization. One of the main ideas in this line of research comes from the field of classical mechanics and proposes to link Nesterov's trajectory to the solution of a set of Euler-Lagrange equations relative to the so-called Bregman Lagrangian. In the last years, this approach led to the discovery of many new (stochastic) accelerated algorithms and provided a solid theoretical foundation for the design of structure-preserving accelerated methods. In this work, we revisit this idea and provide an in-depth analysis of the action relative to the Bregman Lagrangian from the point of view of calculus of variations. Our main finding is that, while Nesterov's method is a stationary point for the action, it is often not a minimizer but instead a saddle point for this functional in the space of differentiable curves. This finding challenges the main intuition behind the variational interpretation of Nesterov's method and provides additional insights into the intriguing geometry of accelerated paths. * Equal Contribution. 2 A differentiable function f : R d → R is said to be β-smooth if it has β-Lipschitz gradients. 3 This lower bound holds just for k < d hence it is only interesting in the high-dimensional setting.

Introduction

This paper focuses on the problem of unconstrained convex optimization, i.e. to find

x * ∈ arg min

x∈R d f (x), (P)
for some lower bounded convex L-smooth 2 loss f ∈ C 1 (R d , R).

Nesterov's acceleration. [START_REF] Nemirovskii | Problem Complexity and Method Efficiency in Optimization[END_REF] showed that no gradient-based optimizer can converge to a solution of (P) faster than O(k -2 ), where k is the number of gradient evaluations 3 . While Gradient Descent (GD) converges like O(k -1 ), the optimal rate O(k -2 ) is achieved by the celebrated Accelerated Gradient Descent (AGD) method, proposed by [START_REF] Nesterov | A method for solving the convex programming problem with convergence rate o(1/k 2 )[END_REF]:

x k+1 = y k -η∇f (y k ) , with y k = x k + k -1 k + 2 (x k -x k-1 ). (AGD)
The intuition behind Nesterov's method and the fundamental reason behind acceleration is, to this day, an active area of research [START_REF] Allen-Zhu | Linear coupling: An ultimate unification of gradient and mirror descent[END_REF][START_REF] Defazio | On the curved geometry of accelerated optimization[END_REF][START_REF] Ahn | From proximal point method to Nesterov's acceleration[END_REF].

ODE models. Towards understanding the acceleration mechanism, [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights[END_REF] made an interesting observation: the convergence rate gap between GD and AGD is retained in the continuoustime limits (as the step-size η vanishes):

Ẋ + ∇f (X) = 0 (GD-ODE), Ẍ + 3 t Ẋ + ∇f (X) = 0 (AGD-ODE)

where Ẋ := dX/dt denotes the time derivative (velocity) and Ẍ := d 2 X/dt 2 the acceleration. Namely, we have that GD-ODE converges like O(t -1 ) and AGD-ODE like O(t -2 ), where t > 0 is the time variable. This seminal paper gave researchers a new tool to understand the nature of accelerated optimizers through Bessel Functions [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights[END_REF], and led to the design of many novel fast and interpretable algorithms outside the Euclidean setting [START_REF] Wibisono | A variational perspective on accelerated methods in optimization[END_REF][START_REF] Wilson | Accelerating rescaled gradient descent: Fast optimization of smooth functions[END_REF], in the stochastic setting [START_REF] Krichene | Accelerated mirror descent in continuous and discrete time[END_REF][START_REF] Xu | Accelerated stochastic mirror descent: From continuous-time dynamics to discrete-time algorithms[END_REF] and also in the manifold setting [START_REF] Alimisis | A continuous-time perspective for modeling acceleration in Riemannian optimization[END_REF][START_REF] Duruisseaux | A variational formulation of accelerated optimization on Riemannian manifolds[END_REF].

Nesterov as solution to Euler-Lagrange equations. It is easy to see that AGD-ODE can be recovered from Euler-Lagrange equations, starting from the time-dependent Lagrangian L(X, Ẋ, t) = t 3 1 2 Ẋ 2 -f (X) .

(1)

Indeed, the Euler-Lagrange equation

d dt ∂ ∂ Ẋ L(X, Ẋ, t) = ∂ ∂X L(X, Ẋ, t) (2)
reduces in this case to t 3 Ẍ + 3t 2 Ẋ + t 3 ∇f (X) = 0, which is equivalent to AGD-ODE (assuming t > 0). In a recent influential paper, [START_REF] Wibisono | A variational perspective on accelerated methods in optimization[END_REF] generalized the derivation above to non-Euclidean spaces, where the degree of separation between points x and y is measured by means of the Bregman Divergence [START_REF] Bregman | The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[END_REF]) D ψ (x, y) = ψ(y) -ψ(x) -∇ψ(x), y -x , where ψ : R d → R is a strictly convex and continuously differentiable function (see e.g. Chapter 1.3.2 in Amari ( 2016)). Namely, they introduced the so-called Bregman Lagrangian:

L α,β,γ (X, Ẋ, t) = e α(t)+γ(t) D ψ (X + e -α(t) V, X) -e β(t) f (X) , (3)
where α, β, γ are continuously differentiable functions of time. The Euler-Lagrange equations imply

Ẍ + (e α(t) -α(t)) Ẋ + e 2α(t)+β(t) ∇ 2 ψ(X + e -α(t) Ẋ) -1 ∇f (X) = 0. (4)
The main result of [START_REF] Wibisono | A variational perspective on accelerated methods in optimization[END_REF] is that, under the ideal-scaling conditions β(t) ≤ e α(t) and γ(t) = e α(t) , any solution to Eq. ( 4) converges to a solution of (P) at the rate O(e -β(t) ). Under the choice ψ(x) = 1 2 x 2 2 , we get back to the Euclidean metric D ψ (x, y) = 1 2 x -y 2 2 . Moreover, choosing α(t) = log(2/t), β(t) = γ(t) = 2 log(t), we recover the original Lagrangian in Eq. ( 1) and O(e -β(t) ) = O(t -2 ), as derived in [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights[END_REF]. Impact of the variational formulation. The variational formulation in [START_REF] Wibisono | A variational perspective on accelerated methods in optimization[END_REF] has had a considerable impact on the recent developments in the theory of accelerated methods. Indeed, this approach can be used to design and analyze new accelerated algorithms. For instance, [START_REF] Xu | Accelerated stochastic mirror descent: From continuous-time dynamics to discrete-time algorithms[END_REF] used the Lagrangian mechanics formalism to derive a novel simplified variant of accelerated stochastic mirror descent. Similarly, [START_REF] França | On dissipative symplectic integration with applications to gradient-based optimization[END_REF], [START_REF] Muehlebach | Optimization with momentum: Dynamical, controltheoretic, and symplectic perspectives[END_REF] used the dual Hamiltonian formalism to study the link between symplectic integration of dissipative ODEs and acceleration. Due to its rising importance in the field of optimization, the topic was also presented by Prof. M. I. Jordan as a plenary lecture at the International Congress of Mathematicians in 2018 [START_REF] Jordan | Dynamical, symplectic and stochastic perspectives on gradient-based optimization[END_REF], centered around the question "what is the optimal way to optimize?". Imprecise implications of the variational formulation. While the Lagrangian formalism has been inspiring and successful for algorithm design and analysis, its precise implications for the geometry and the path of accelerated solutions have not been examined in a mathematically rigorous way (to the best of our knowledge). In [START_REF] Jordan | Dynamical, symplectic and stochastic perspectives on gradient-based optimization[END_REF] it is hinted that, since Nesterov's method solves the Euler-Lagrange equations, it minimizes the action functional t2 t1 L α,β,γ (Y, Ẏ , t)dt over the space of curves by the minimum action principle of classical mechanics [START_REF] Arnol'd | Mathematical methods of classical mechanics[END_REF]. This claim4 is inaccurate. Indeed, the term minimum action principle is misleading5 : solving Euler-Lagrange only makes the action stationary (necessary condition: vanishing first-order derivative), but does not guarantee minimality -this only holds in physics for very special cases6 , which do not include even simple mechanical systems like the pendulum (proof in Section 36.2 of [START_REF] Gelfand | Calculus of variations[END_REF]). Indeed, from a theoretical perspective, the claim requires computing the second variation along Nesterov's path. Quite surprisingly, even though many papers are dedicated to the variational formulation [START_REF] Wibisono | A variational perspective on accelerated methods in optimization[END_REF][START_REF] Jordan | Dynamical, symplectic and stochastic perspectives on gradient-based optimization[END_REF][START_REF] Casgrain | A latent variational framework for stochastic optimization[END_REF][START_REF] Duruisseaux | A variational formulation of accelerated optimization on Riemannian manifolds[END_REF], to the best of our knowledge there is no work which provides an in-depth rigorous study of the action relative to Bregman Lagrangian and that characterizes minimality of Nesterov in the space of curves. Our contributions. Intrigued by the non-trivial open question of minimality of Nesterov's path and by the enigmatic geometry of accelerated flows, in this paper, we examine the properties of accelerated gradient methods from the perspective of calculus of variations. 1. In Sec. 3 we study the minimality of classical Nesterov's ODE (damping 3/t) proposed by [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights[END_REF] on multidimensional quadratic losses. By using Jacobi's theory for the second variation (summarized in Sec. 2), we find that Nesterov's path is optimal only if the integration interval [t 1 , t 2 ] is small enough. In contrast, if t 2 -t 1 > 40/β (β is Lipschitz constant for the gradient), Nesterov's path is actually a saddle point for the action (see Fig. 1). 2. In Sec. 4 we extend the analysis to the µ-strongly convex setting and thus consider a constant damping α. We show that, for extremely overdamped Nesterov flows (α ≥ 2 √ β), i.e for highly suboptimal parameter tuning (acceleration holds only for α ≈ 2 √ µ), Nesterov's path is always a minimizer for the action. In contrast, we show that for α < 2 √ β (acceleration setting), if t 2 -t 1 > 2π/ √ 4β -α, Nesterov's path is again a saddle point for the action. 3. In Sec. 5 we discuss the implications of our results for the theory of accelerated methods and propose a few interesting directions for future research.

We start by recalling some definitions and results from calculus of variations, which we adapt from classical textbooks [START_REF] Curran Associates | [END_REF][START_REF] Arnol'd | Mathematical methods of classical mechanics[END_REF][START_REF] Gelfand | Calculus of variations[END_REF].

Background on calculus of variations

We work on the vector space of curves

C 1 ([t 1 , t 2 ], R d ) with t 1 , t 2 ∈ [0, ∞).
We equip this space with the standard norm Y = max t1≤t≤t2 Y (t) 2 + max t1≤t≤t2 Ẏ (t) 2 . Under this choice, for any regular Lagrangian L, the action functional

J[Y ] := t2 t1 L(Y, Ẏ , t)dt is continuous. First variation. Let D be the linear subspace of continuously differentiable displacements curves h such that h(t 1 ) = h(t 2 ) = 0. The corresponding increment of J at Y along h is defined as ∆J[Y ; h] := J[Y + h] -J[Y ]. Suppose that we can write ∆J[Y ; h] = ϕ[Y ; h] + h ,
where ϕ is linear in h and → 0 as h → 0. Then, J is said to be differentiable at Y and the linear functional

δJ[Y ; •] : D → R such that δJ[Y ; h] := ϕ[Y ; h] is called first variation of J at Y . It can be shown that, if J is differentiable at Y , then its first variation at Y is unique. Extrema. J is said to have an extremum at Y if ∃δ > 0 such that, ∀h ∈ D with h ≤ δ, the sign of J[Y + h] -J[Y ] is constant. A necessary condition for J to have an extremum at Y is that δJ[Y ; h] = 0, for all h ∈ D. (5)
We proceed with stating one of the most well-known results in calculus of variations, which follows by using Taylor's theorem on J

[Y ] = t2 t1 L(Y, Ẏ , t) dt. Theorem 1 (Euler-Lagrange equation). A necessary condition for the curve Y ∈ C 1 ([t 1 , t 2 ], R d ) to be an extremum for J (w.r.t. D) is that it satisfies the Euler-Lagrange equations (2).
It is crucial to note that Theorem 1 provides a necessary, but not sufficient condition for an extremum -indeed, the next paragraph is completely dedicated to this. Second Variation. Thm. 1 does not distinguish between extrema (maxima or minima) and saddles. For this purpose, we need to look at the second variation.

Suppose that the increment of J at Y can be written as

∆J[Y ; h] = ϕ 1 [Y ; h] + ϕ 2 [Y ; h] + h 2 ,
where ϕ 1 is linear in h, ϕ 2 is quadratic in h and → 0 as h → 0. Then J is said to be twice differentiable and the functional

δ 2 J[Y :, •] : D → R s.t. δ 2 J[Y, h] := ϕ 2 [Y ; h] is called the second variation of J at Y .
Uniqueness of second variation is proved in the same way as the first variation. Theorem 2. A necessary condition for the curve

Y ∈ C 1 ([t 1 , t 2 ], R d ) to be a local minimum for J (w.r.t D) is that it satisfies δ 2 J[Y ; h] ≥ 0. For local maxima, the sign is flipped. Jacobi equations. Recall that J[Y ] = t2 t1 L(Y, Ẏ , t)dt. Using the notation L Y Z = ∂ 2 L/(∂Y ∂Z), the Taylor expansion for ∆J[Y ; h] = J[Y + h] -J[Y ] if h → 0 converges to ∆J[Y ; h] = t2 t1 L Y h + L Ẏ ḣ dt + 1 2 t2 t1 L Y Y h 2 + L Ẏ Ẏ ḣ2 + 2L Y Ẏ h ḣ dt, (6)
where the equality holds coordinate-wise. Therefore, δJ[Y

; h] = t2 t1 L Y h + L Ẏ ḣ dt and δ 2 J[Y ; h] = 1 2 t2 t1 L Y Y h 2 + 2L Y Ẏ h ḣ + L Ẏ Ẏ ḣ2 dt = 1 2 t2 t1 L Y Y - d dt L Y Ẏ h 2 dt + 1 2 t2 t1 L Ẏ Ẏ ḣ2 dt = 1 2 t2 t1 P ḣ2 + Qh 2 dt, (7)
where

P = L Ẏ Ẏ , Q = L Y Y -d dt L Y Ẏ
, and the second equality follows from integration by parts since h vanishes at t 1 and t 2 . Using this expression, it is possible to derive an easy necessary (but not sufficient) condition for minimality. Theorem 3 (Legendre's necessary condition). A necessary condition for the curve Y to be a minimum of J is that L Ẏ Ẏ is positive semidefinite. Conjugate points. A crucial role in the behavior of δ 2 J[Y ; h] is played by the shape of the solutions to Jacobi's differential equation d dt (P ḣ) -Qh = 0. A point t ∈ (t 1 , t 2 ) is said to be conjugate to point t 1 (w.r.t. J) if Jacobi's equation admits a solution which vanishes at both t 1 and t but is not identically zero. We have the following crucial result. Theorem 4 (Jacobi's condition). Necessary and sufficient conditions for Y to be a local minimum for J are: (1) Y satiesfies the Euler-Lagrange Equation ;(2) P positive definite; (3) (t 1 , t 2 ) contains no points conjugate to t 1 .

Analysis of the action of Nesterov's path with vanishing damping 3/t

This section is dedicated to the analysis of the action functional relative to Eq. (3). By incorporating the tool of variational calculus, we study the optimality of Nesterov's method in minimizing the action. We start by a general abstract analysis in the convex quadratic case in Sec. 3.1, and then present an intuitive analytical computation in Sec. 3.2. The non-quadratic case is discussed in Sec. 3.4.

Solutions to Jacobi's equation for the Bregman Lagrangian in the quadratic setting

For the sake of clarity, we start by considering the Lagrangian in Eq. ( 1) for the simple onedimensional case f (x) = βx 2 /2. We have Therefore, Jacobi's equation relative to the action functional

Q = L Y Y - d dt L Y Ẏ = -βt 3 , P = L Ẏ Ẏ = t 3 . ( 8
)
t2 t1 L(Y, Ẏ , t)dt with t 1 > 0 is d dt (t 3 ḣ) -βt 3 h = 0. =⇒ t 3 ḧ + 3t 2 ḣ + βt 3 h = 0 =⇒ ḧ + 3 t ḣ + βh = 0, (9)
which is itself Nesterov's ODE. Following the procedure outlined in Theorem 2, we now study the solutions h such that h(t 1 ) = 0. Any solution to Eq. ( 9) can be written as7

h(t) = C Y 1 ( √ β t) t -C Y 1 ( √ β t 1 ) J 1 ( √ β t) J 1 ( √ β t 1 ) t , (10)
where C > 0 specifies the initial velocity (see Fig. 2), J α is the Bessel function of the first kind and Y α is the Bessel function of the second kind.

J α (x) = ∞ m=0 (-1) m m! Γ(m + α + 1) x 2 2m+α , Y α (x) = J α (x) cos(απ) -J -α (x) sin(απ) . ( 11
) Points t > t 1 conjugate to t 1 satisfy h(t) = 0, which results in the identity Y 1 ( √ β t)/Y 1 ( √ β t 1 ) = J 1 ( √ β t)/J 1 ( √ β t 1 ).
Remarkably, this condition does not depend on C, but only on t 1 and on the sharpness β. Let us now fix these parameters and name

K β,t1 = Y 1 ( √ β t 1 )/J 1 ( √ β t 1 ). Points conjugate to t 1 then satisfy Y 1 ( √ β t) = K β,t1 J 1 ( √ β t)
. We now recall the following expansions [START_REF] Watson | A treatise on the theory of Bessel functions[END_REF], also used by [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights[END_REF]:

J 1 (x) = 2 πx cos x - 3π 4 + O 1 x , Y 1 (x) = 2 πx sin x - 3π 4 + O 1 x .
(12) Since J 1 and Y 1 asymptotically oscillate around zero and are out of synch (π/2 difference in phase), for t big enough the condition Y 1 (

√ β t) = K β,t1 J 1 ( √ β t
) is going to be satisfied. Further, this condition is going to be satisfied for a smaller value for t if β is increased, as confirmed by Figure 2.

Theorem 5 (Local optimality of Nesterov with vanishing damping). Let f : R d → R be a convex quadratic, and let X : R → R d be a solution to the ODE Ẍ +

3 t Ẋ + ∇f (X) = 0. For 0 < t 1 < t 2 , consider the action functional J[Y ] = t2 t1 L(Y, Ẏ , t)dt, mapping Y ∈ C 1 ([t 1 , t 2 ], R d ), to a real number. Then, if |t 2 -t 1 |
is small enough, there are no points conjugate to t 1 and Nesterov's path minimizes J over all curves such that Y (t 1 ) = X(t 1 ) and Y (t 2 ) = X(t 2 ). The length of the optimality interval |t 2 -t 1 | shrinks as β, the maximum eigenvalue of the Hessian of f , increases.

Proof. The argument presented in this section can be lifted to the multidimensional case. Indeed, since the dynamics in phase space is linear, it's geometry is invariant to rotations and we can therefore assume the Hessian is diagonal. Next, Jacobi's equation has to be solved coordinate-wise, which leads to a logical AND between conjugacy conditions. By the arguments above, the dominating condition is the one relative to the maximum eigenvalue β.

The following corollary shows that Nesterov's path actually becomes suboptimal if the considered time interval is big enough. This is also verified numerically in Figure 1.

Corollary 6 (Nesterov with vanishing damping is not globally optimal). In the settings of Thm. 5, for |t 2 -t 1 | big enough, Nesterov's path becomes a saddle point for J.

Proof. Non-existence of conjugate points is necessary and sufficient for minimality/maximality.

In the next subsection, we provide a constructive proof for Cor. 6, which allows us to derive a concrete bound for |t 2 -t 1 |. In particular, we show in Prop. 7 that, for the case of vanishing damping 3/t, Nesterov's path is always a saddle point for the action if |t 2 -t 1 | > 40/β. Optimality for the special case X(t 2 ) = x * . Some readers might have already realized that Jacobi's equation Eq. ( 9) in the quadratic potential case is itself the solution of Nesterov's ODE. This means that, if t 1 ≈ 0, the first time conjugate to t 1 is exactly when Nesterov's path reaches the minimizer. Hence -only in the one-dimensional case -it is actually true that, if we do not consider arbitrary time intervals but only the first interval before the solution first touches the optimizer, Nesterov's path always minimizes the action. Sadly, this strong and interesting result is not valid in higher dimensions, since Nesterov's path in general never actually crosses the minimizer in finite time8 . Remark on dropping boundary conditions. The results in this section are formulated for the fixed boundaries case Y (t 1 ) = X(t 1 ) and Y (t 2 ) = X(t 2 ), where X is the solution to Nesterov's ODE. From an optimization viewpoint, this requirement seems strong. Ideally, we would want X(t 2 ) to be any point close to the minimizer (say inside an -ball). A simple reasoning proves that Nesterov's path can be a saddle point for the action also in this case. By contradiction, assume Nesterov's trajectory minimizes the action among all curves that reach any point inside a small -ball at time t 2 . Nesterov's path also minimizes the action in the (smaller) space of curves that reach exactly X(t 2 ). By Cor. 6, this leads to a contradiction if t 2 is big enough.

A constructive proof of Nesterov's suboptimality in the quadratic setting

We now present a direct computation, to shed some light on the suboptimality of Nesterov's path in the context of Thm. 5. In the setting of Sec. 3.1, the second variation of J along γ is 1 2 t2 t1 t 3 [ ḣ(t) 2βh(t) 2 ]dt, independent of γ. Consider now the finite-norm perturbation (vanishing at boundary):

h ,c (t) =    0 t ≤ c -or t ≥ c + t-c+ t ∈ (c -, c) c+ -t t ∈ (c, c + ) , c ∈ (t 1 , t 2 ), < min(c -t 1 , t 2 -c). ( 13
)
This is a triangular function with support (c -, c + ) and height one. Let h ,c be a C 1 modification of h ,c such that h ,c -h ,c is negligible9 . For any scaling factor σ > 0,

δ 2 J(σ • h ,c ) = -σ 2 3β 4 10 + (βc 2 -3) 2 -3c 2 c 3 . ( 14
)
The denominator is always positive. Hence, we just need to study the sign of the numerator, with respect to changes in > 0 and c > 0. Consider for now c fixed, then the zeros of the numerator are at (15 -5u ± √ 25u 2 -60u + 225)/(3β), with u := βc 2 . Since 25u 2 -60u + 225 > 0 for all u ≥ 0, the solution has two real roots. However, only one root 2 * (u, β) is admissible, since the smallest one is always negative 10 . As a result, for fixed c > 0, δ 2 J(σh ,c ) changes sign only at 2 * (u, β). Note that 2 * (u, β) is decreasing as a function of u and 2 * (0, β) = 10/β as well as lim u→∞ 2 * (u, β) = 3/β. Therefore, for any c, β > 0, we showed that δ 2 J(σh ,c ) changes sign when * ∈ [ 3/β, 10/β]. If choosing big is allowed by the considered interval (h ,c has to vanish at t 1 , t 2 ), then the second variation is indefinite. This happens if |t 2 -t 1 | > 2 * . By taking σ → 0, we get the following result.

Proposition 7 (Sufficient condition for saddle). The second variation of the action of Nesterov's Lagrangian (damping 3/t) on f (x) = βx 2 /2 is an indefinite quadratic form for |t 2 -t 1 | > 40/β. This result generalizes to β-smooth multidimensional convex quadratics.

We remark that the inverse dependency on the square root of β is also predicted by the general proof in Sec. 3.1, where the argument of the Bessel functions is always t √ β.

3.3 Unboundedness of the action for large integration intervals (from above and below)

In Prop. 7, we showed that for big enough integration intervals, Nesterov's method with damping 3/t on f (x) = βx 2 /2 is saddle point for the action. This suggests that the action is itself unboundedboth from above and below. It is easy to show this formally.

Proposition 8 (Unboundedness of the action). Let L be Lagrangian of Nesterov's method with damping 3/t on a β-smooth convex quadratic and

J[Y ] = t2 t1 L(Y, Ẏ , t)dt. Let a, b be two arbitrary vectors in R d . There exists a sequence of curves (Y k ) k∈N , with Y k ∈ C 1 ([t 1 , t 2 ], R d ) and Y k (t 1 ) = a, Y k (t 2 ) = b for all k ∈ N, such that J[Y k ] k → ∞. In addition, if |t 2 -t 1 | > 40/β
there exists another sequence with the same properties diverging to -∞.

Proof. The proof is based on the computation performed for Prop. 7. Crucially, note that for the quadratic loss function case we have δ 2 J = J. For the case a = b = 0, we showed that for any interval [t 1 , t 2 ], by picking small enough, we have J(h ,c ) = δ 2 J(h ,c ) > 0 (also illustrated in the figure supporting the proof). Hence, J(σ • h ,c ) → +∞ as σ → ∞. Same argument holds for -∞ in the large interval case. This proves the assertion for a = b = 0. Note that the curves corresponding to the diverging sequences can be modified to start/end at any a, b ∈ R d at the price of a bounded error in the action. This does not modify the behavior in the limit; hence, the result follows.

Optimality of Nesterov with vanishing damping if curvature vanishes (polynomial loss)

Note that the bound on |t 2 -t 1 | in Prop. 7 gets loose as β decreases. This is also predicted by the argument with Bessel functions in Sec. 3.1, and clear from the simulation in Fig. 2. As a result, as curvature vanishes, Nesterov's path becomes optimal for larger and larger time intervals. This setting is well described by polynomial losses f (x) ∝ (x -x * ) p , with p > 2. As Nesterov's path approaches the minimizer x * , the curvature vanishes; hence, for every β > 0 there exists a time interval (t 1 , ∞) where the curvature along Nesterov's path is less then β. This suggest that, for losses with vanishing curvature at the solution, there exists a time interval (t * , ∞) where Nesterov's path is actually a minimizer for the action. While this claim is intuitive, it is extremely hard to prove formally since in this case the second variation of J depends on the actual solution of Nesterov's equations -for which no closed-form formula is known in the polynomial case [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights[END_REF].

However, we also note that the vanishing sharpness setting is only interesting from a theoretical perspective. Indeed, regularized machine learning objectives do not have this property. Actually, in the deep neural network setting, it is known that the sharpness (maximum eig. of the Hessian) actually increases overtime [START_REF] Yao | Pyhessian: Neural networks through the lens of the Hessian[END_REF][START_REF] Cohen | Gradient descent on neural networks typically occurs at the edge of stability[END_REF]. Hence, it is safe to claim that in the machine learning setting Nesterov's path is only optimal for small time intervals, as shown in Thm. 5.

Analysis of the action of Nesterov's path with constant damping α

In the µ-strongly convex case 11 , it is well known that a constant damping α = 2 √ µ yields acceleration compared to gradient descent 12 . This choice completely changes the geometry of Nesterov's path and needs a separate discussion. The corresponding Lagrangian is

L α (Y, Ẏ , t) = e αt 1 2 Ẏ 2 -f (Y ) . ( 15
)
11 Hessian eigenvalues lower bounded by µ > 0.

12 The corresponding rate is linear and depends on µ/β, as opposed to µ/β (GD case).

Again, we consider the quadratic function f (x) = βx 2 /2 and examine Jacobi's ODE ḧ(t) + α ḣ(t) + βh(t) = 0. We have to determine whether there exists a non-trivial solution such that h(t 1 ) = h(t 2 ) = 0 and h(t) vanishes also at a point t ∈ (t 1 , t 2 ), the conjugate point.

For the critical damping case α = 2 √ β the general solution such that h(t 1 ) = 0 is

h(t) = Ce - √ βt (t -t 1 ). (16)
There is no non-trivial solution h that vanishes also at t ∈ (t 1 , t 2 ) -no conjugate points. The same holds for the overdamping case α > 2 √ β, where the solution that vanishes at t 1 is

h(t) = Ce -αt 2 e 1 2 √ α 2 -4βt -e 1 2 √ α 2 -4β(2t1-t) . (17)
For the underdamping case α < 2 √ β, the picture gets more similar to the vanishing damping case (Sec.3.1). The solution under h(t 1 ) = 0 is

h(t) = Ce -αt 2 sin 4β -α 2 2 t -tan 4β -α 2 2 t 1 cos 4β -α 2 2 t . (18)
Hence all points t > t 1 conjugate to t 1 satisfy t = t 1 + 2kπ/ 4β -α 2 for k ∈ N. Therefore for any t 2 > t 1 + 2π/ 4β -α 2 there exists a conjugate point t ∈ (t 1 , t 2 ).

Theorem 9 (Global optimality of overdamped Nesterov, suboptimality of accelerated Nesterov).

Let f : R d → R be a strongly convex quadratic, and let X : R → R d be a solution to the ODE

Ẍ + α Ẋ + ∇f (X) = 0. For 0 ≤ t 1 < t 2 , consider the action J[Y ] = t2 t1 L α (Y, Ẏ , t)dt. If α ≥ 2 √ β,
where β is the max. eigenvalue of the Hessian of f , then Nesterov's path minimizes J over all curves s.t. Y (t 1 ) = X(t 1 ) and Y (t 2 ) = X(t 2 ). Else (e.g. acceleration setting α ≈ 2 √ µ),

Nesterov's path is optimal only for |t 2 -t 1 | ≤ 2π/ 4β -α 2 and otherwise is a saddle point.

Proof. As for the proof of Thm. 5, the condition on conjugate points has to hold for each eigendirection separately. We conclude by noting that eigenvalues are in the range [µ, β].

For the underdamping case, we give a concrete example for α = β = 1, to show the saddle point nature. Consider the finite-norm perturbation h(t) = sin(kπ(t -t 1 )/(t 1 -t 2 )), where k ∈ N. Then,

δ 2 J[γ](σh) = σ 2 e 2t1 k 2 e -(t2-t1) π 2 (e t2-t1 -1)(2k 2 π 2 -(t 2 -t 1 ) 2 ) (t 2 -t 1 ) 2 (4k 2 π 2 + (t 2 -t 1 ) 2 ) . (19)
Hence, for any t 2 -t 1 > √ 2kπ, it holds that δ 2 J[γ](σh) < 0.

Extending the optimality claims to t 2 = ∞ with Γ-convergence. From an optimization viewpoint, the most interesting setting is to study the action over the complete trajectory, i.e. to consider

Y ∈ C 1 ([t 1 , ∞), R d ) such that Y (t 1 ) = X(t 1 )
and Y (∞) = x * , a minimizer. Prop. 7 and Thm. 9 show that the question of optimality in this case deserves a discussion only in the extremely overdamped case α ≥ 2 √ β, where minimality is guaranteed for any time interval. A careful study of the infinite-time setting would require the theory of Γ-convergence [START_REF] Braides | Gamma-convergence for Beginners[END_REF]. The usual pipeline consists in defining a sequence of problems J k , on intervals [t 1 , t k 2 ], with t k 2 → ∞ as k → ∞. Under the assumption that each J k admits a global minimizer (only true for the overdamped case), one can study convergence of

J * k = min{J k [Y ] : Y ∈ C 1 ([t 1 , t k 2 ], R d )} to J * ∞ = min{J ∞ (Y ) : Y ∈ C 1 ([t 1 , ∞), R d )}. While existence of J *
∞ and J ∞ is not trivial in general, for our setting the pipeline directly yields minimality of overdamped Nesterov's path until ∞.

Discussion of the main findings and directions for future research

In this section, we summarize the results of Sec. 3 & 4 and discuss some implications that our analysis delivers on the geometry of accelerated flows in the convex and strongly convex setting.

We summarize below the main high-level findings of our theoretical analysis:

1. The optimality of Nesterov's path for minimization of the action corresponding to the Bregman Lagrangian is strictly linked to the curvature around the minimizer reached by the flow.

2. As the maximal curvature β increases, it gets increasingly difficult for accelerated flows to minimize the action over long integration intervals: both the accelerated ODEs Ẍ + 3/t Ẋ + ∇f (X) = 0 and Ẍ + 2 √ µ Ẋ + ∇f (X) = 0 are optimal only for intervals of length ∝ 1/ √ β. 3. If Nesterov's path does not minimize the action, there does not exist a "better" (through the eyes of the action) algorithm, as the functional gets unbounded from below (Sec. 3.3). 4. This suboptimality is due precisely to the oscillations in the accelerated paths. In contrast, as long as each coordinate in parameter space decreases monotonically (as it is the case for gradient descent), Nesterov's path is optimal. See also O'donoghue and Candes (2015). 5. Hence, Nesterov's method with very high damping α > 2 √ β -which does not oscillate and hence does not lead to acceleration (see Fig. 3) -minimizes the action. In a nutshell, locally Nesterov's method does indeed optimize a functional over curves. However, this property breaks down precisely when the geometry gets interesting -i.e. when the loss evolution is non-monotonic. Since acceleration is a global phenomenon, i.e. is the cumulative result of many consecutive oscillations (see Fig. 3), our results suggest that the essence of acceleration cannot be possibly captured by the minimization of the action relative to the Bregman Lagrangian. Non-uniqueness of the Lagrangian. A possible reason for the non-optimality of Nesterov's path is, simply put -that we are not looking at the right action functional. Indeed, there are many Lagrangians that can generate Nesterov ODE. Let F (Y, t) be any function which does not involve the velocity, then it is easy to see that the Lagrangian L is equivalent to L(Y, Ẏ , t) = L(Y, Ẏ , t) + Ẏ , ∂F ∂Y (X, t) + ∂F ∂t (X, t).

(20) This simple fact opens up new possibilities for analyzing and interpreting Nesterov's method using different functionals -which perhaps have both a more intuitive form and better properties. Higher order ODEs. On a similar note, it could be possible to convexify the action functional by considering a logical OR between symmetric ODEs, e.g ( d 2 dt 2 + α d dt + β)( d 2 dt 2 -α d dt + β)X = 0. Such tricks are often used in the literature on dissipative systems [START_REF] Szegleti | Dissipation in Lagrangian formalism[END_REF]. Noether Theorem. In physics, the variational framework is actually never used to claim the minimality of the solution to the equations of motion. Its power relies almost completely in the celebrated Noether's Theorem [START_REF] Noether | Invariante variationsprobleme[END_REF], which laid the foundations for modern quantum field theory by linking the symmetries in the Lagrangian (and of the Hamiltonian) to the invariances in the dynamics. Crucially, for the application of Noether's Theorem, one only needs the ODE to yield a stationary point for the action (also saddle points work). Coincidentally, while finalizing this manuscript, two preprints [START_REF] Tanaka | Noether's learning dynamics: The role of kinetic symmetry breaking in deep learning[END_REF]Głuch and Urbanke, 2021) came out on some implications of Noether's Theorem for optimization. However, we note that these works do not discuss the direct link between Noether's Theorem and acceleration, but instead study the interaction between the symmetries in neural network landscapes and optimizers. While some preliminary implications of Noether's theorem for time-rescaling of accelerated flows are discussed in [START_REF] Wibisono | A variational perspective on accelerated methods in optimization[END_REF], we suspect that a more in-depth study could lead, in combination with recent work on the Hamiltonian formalism [START_REF] Diakonikolas | Generalized momentum-based methods: A Hamiltonian perspective[END_REF], to substantial insights on the hidden invariances of accelerated paths. We note that finding these invariances might not be an easy task, and requires a dedicated work: indeed, even for simple linear damped harmonic oscillators (constant damping), invariance in the dynamics can be quite complex [START_REF] Choudhuri | Symmetries and conservation laws of the damped harmonic oscillator[END_REF].

Conclusion

We provided an in-depth theoretical analysis of Nesterov's method from the perspective of calculus of variations, and showed that accelerated paths only minimize the action of the Bregman Lagrangian locally. This suggests that further research is needed for understanding the underlying principle of generation for acceleration -which has been an open question in optimization for almost 40 years [START_REF] Nesterov | A method for solving the convex programming problem with convergence rate o(1/k 2 )[END_REF]. To this end, we proposed a few concrete directions for future work.

Figure 1 :

 1 Figure1: Optimization of f (x) = x 2 /2 using AGD-ODE. Peturbations (vanishing at extrema) are added to the AGD-ODE solution: depending on the perturbation kind (i.e. direction in space of curves), the local behavior is either a max or a min. Hence, Nesterov's path can be a saddle point for the action (formally shown in Sec. 3.1).

Figure 2 :

 2 Figure 2: First conjugate point to t1 = 1, 4 for quadratics βx 2 /2, under the settings of Section 3.1. For each value of β, six solutions h(t) to the Jacobi equation (each one has different velocity) are shown.

Figure 3 :

 3 Figure 3: Optimization of potential f (x, y) = 0.02x 2 + 0.0004y 2 , where 2e -2 = β µ = 3e -4. Nonmonotonic trajectories (i.e. the accelerated curves) minimize the action only for short time intervals. Simulation with Runge-Kutta 4 integration.

Paragraph before Eq. (9) in[START_REF] Jordan | Dynamical, symplectic and stochastic perspectives on gradient-based optimization[END_REF]: "[...] we use standard calculus of variations to obtain a differential equation whose solution is the path that optimizes the time-integrated Bregman Lagrangian".

From Section 36.2 in[START_REF] Gelfand | Calculus of variations[END_REF]: "The principle of least action is widely used[...]. However, in a certain sense the principle is not quite true [...]. We shall henceforth replace the principle of least action by the principle of stationary action. In other words, the actual trajectory of a given mechanical system will not be required to minimize the action but only to cause its first variation to vanish."

e.g. free particle in vanishing potentials, or t1 ≈ t2, see Remark 2 at the end of Section 21 and Section 36.2.

Symbolic computations are checked in Maple/Mathematica, numerical simulations are performed in Matlab.

The first crossing time in each direction depends on the sharpness in each direction, hence by the time each coordinate reaches zero, we already have a conjugate point.

Standard technique in calculus of variation, see e.g. proof of Legendre's Thm[START_REF] Gelfand | Calculus of variations[END_REF].

If u > 0, then 15 -5u -√ 25u 2 -60u + 225 < 0.

Acknowledgments and Disclosure of Funding

The authors thank Prof. Luca Martinazzi for his helpful discussions. This work was partially funded in part by the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001(PRAIRIE 3IA Institute). Hadi Daneshmand acknowledges support from the European Research Council (grant SEQUOIA 724063).


Hadi Daneshmand
email: seyed.daneshmand@inria.fr

Jonas Kohler
email: jonas.kohler@inf.ethz.ch

Francis Bach
email: francis.bach@inria.fr

Thomas Hofmann
email: thomas.hofmann@inf.ethz.ch

Aurelien Lucchi
email: aurelien.lucchi@inf.ethz.ch

Batch Normalization Provably Avoids Rank Collapse for Randomly Initialised Deep Networks

Randomly initialized neural networks are known to become harder to train with increasing depth, unless architectural enhancements like residual connections and batch normalization are used. We here investigate this phenomenon by revisiting the connection between random initialization in deep networks and spectral instabilities in products of random matrices. Given the rich literature on random matrices, it is not surprising to find that the rank of the intermediate representations in unnormalized networks collapses quickly with depth. In this work we highlight the fact that batch normalization is an effective strategy to avoid rank collapse for both linear and ReLU networks. Leveraging tools from Markov chain theory, we derive a meaningful lower rank bound in deep linear networks. Empirically, we also demonstrate that this rank robustness generalizes to ReLU nets. Finally, we conduct an extensive set of experiments on real-world data sets, which confirm that rank stability is indeed a crucial condition for training modern-day deep neural architectures.

Introduction and related work

Depth is known to play an important role in the expressive power of neural networks [START_REF] Telgarsky | Benefits of depth in neural networks[END_REF]. Yet, increased depth typically leads to a drastic slow down of learning with gradient-based methods, which is commonly attributed to unstable gradient norms in deep networks [START_REF] Hochreiter | The vanishing gradient problem during learning recurrent neural nets and problem solutions[END_REF]. One key aspect of the training process concerns the way the layer weights are initialized. When training contemporary neural networks, both practitioners and theoreticians advocate the use of randomly initialized layer weights with i.i.d. entries from a zero mean (Gaussian or uniform) distribution. This initialization strategy is commonly scaled such that the variance of the layer activation stays constant across layers [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF][START_REF] He | Delving deep into rectifiers: Surpassing human-level performance on imagenet classification[END_REF]. However, this approach can not avoid spectral instabilities as the depth of the network increases. For example, [START_REF] Saxe | Exact solutions to the nonlinear dynamics of learning in deep linear neural networks[END_REF] observes that for linear neural networks, such initialization lets all but one singular values of the last layers activation collapse towards zero as the depth increases.

Nevertheless, recent advances in neural architectures have allowed the training of very deep neural networks with standard i.i.d. initialization schemes despite the above mentioned shortcomings. Among these, both residual connections and normalization layers have proven particularly effective and are thus in widespread use (see [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF][START_REF] Salimans | Weight normalization: A simple reparameterization to accelerate training of deep neural networks[END_REF][START_REF] He | Delving deep into rectifiers: Surpassing human-level performance on imagenet classification[END_REF] to name just a few). Our goal here is to bridge the explanatory gap between these two observations by studying the effect of architectural enhancements on the spectral properties of randomly initialized neural networks. We also provide evidence for a strong link of the latter with the performance of gradient-based optimization algorithms.

One particularly interesting architectural component of modern day neural networks is Batch Normalization (BN) [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF]. This simple heuristics that normalizes the pre-activation of hidden units across a mini-batch, has proven tremendously effective when training deep neural networks with gradient-based methods. Yet, despite of its ubiquitous use and strong empirical benefits, the research community has not yet reached a broad consensus, when it comes to a theoretical explanation for its practical success. Recently, several alternatives to the original "internal covariate shift" hypothesis [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF] have appeared in the literature: decoupling optimization of direction and length of the parameters [START_REF] Kohler | Exponential convergence rates for batch normalization: The power of length-direction decoupling in non-convex optimization[END_REF], auto-tuning of the learning rate for stochastic gradient descent [START_REF] Arora | Theoretical analysis of auto rate-tuning by batch normalization[END_REF], widening the learning rate range [START_REF] Bjorck | Understanding batch normalization[END_REF], alleviating sharpness of the Fisher information matrix [START_REF] Karakida | The normalization method for alleviating pathological sharpness in wide neural networks[END_REF], and smoothing the optimization landscape [START_REF] Santurkar | How does batch normalization help optimization?[END_REF]. Yet, most of these candidate justifications are still actively debated within the community. For example, [START_REF] Santurkar | How does batch normalization help optimization?[END_REF] first made a strong empirical case against the original internal covariate shift hypothesis. Secondly, they argued that batch normalization simplifies optimization by smoothing the loss landscape. However, their analysis is on a per-layer basis and treats only the largest eigenvalue. Furthermore, even more recent empirical studies again dispute these findings, by observing the exact opposite behaviour of BN on a ResNet20 network [START_REF] Yao | PyHessian: Neural networks through the lens of the Hessian[END_REF].

On random initialization and gradient based training

In light of the above discussion, we take a step back -namely to the beginning of training -to find an interesting property that is provably present in batch normalized networks and can serve as a solid basis for a more complete theoretical understanding.

The difficulty of training randomly initialized, un-normalized deep networks with gradient methods is a long-known fact, that is commonly attributed to the so-called vanishing gradient effect, i.e., a decreasing gradient norm as the networks grow in depth (see, e.g., [START_REF] Samuel | Deep information propagation[END_REF]). A more recent line of research tries to explain this effect by the condition number of the input-output Jacobian (see, e.g., [START_REF] Yang | Mean field residual networks: On the edge of chaos[END_REF][START_REF] Yang | A mean field theory of batch normalization[END_REF][START_REF] Pennington | The emergence of spectral universality in deep networks[END_REF][START_REF] Bjorck | Understanding batch normalization[END_REF]). Here, we study the spectral properties of the above introduced initialization with a particular focus on the rank of the hidden layer activations over a batch of samples. The question at hand is whether or not the network preserves a diverse data representation which is necessary to disentangle the input in the final classification layer.

As a motivation, consider the results of Fig. 1, which plots accuracy and output rank when training batch-normalized and un-normalized neural networks of growing depth on the Fashion-MNIST dataset [START_REF] Xiao | Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms[END_REF]. As can be seen, the rank in the last hidden layer of the vanilla networks collapses with depth and they are essentially unable to learn (in a limited number of epochs) as soon as the number of layers is above 10. The rank collapse indicates that the direction of the output vector has become independent of the actual input. In other words, the randomly initialized network no longer preserves information about the input. Batch-normalized networks, however, preserve a high rank across all network sizes and their training accuracy drops only very mildly as the networks reach depth 32.

The above example shows that both rank and optimization of even moderately-sized, unnormalized networks scale poorly with depth. Batch-normalization, however, stabilizes the rank in this setting and the obvious question is whether this effect is just a slow-down or even simply a numerical phenomenon, or whether it actually generalizes to networks of infinite depth.

In this work we make a strong case for the latter option by showing a remarkable stationarity aspect of BN. Consider for example the case of passing N samples x i ∈ R d arranged column-wise in an input matrix X ∈ R d×N through a very deep network with fully-connected layers. Ideally, from an information propagation perspective, the network should be able to differentiate between individual samples, regardless of its depth [START_REF] Samuel | Deep information propagation[END_REF]. However, as can be seen in Fig. 2, the hidden representation of X collapses to a rank one matrix in vanilla networks, thus mapping all x i to the same line in R d . Hence, the hidden layer activations and along with it the individual gradient directions become perceptrons (MLPs) of depth 1-32 and width 128 hidden units. Left: Rank 3 after random initialization as in PyTorch [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF]. Right: Training accuracy after training 75 epochs with SGD, batch size 128 and grid-searched learning rate. Mean and 95% confidence interval of 5 independent runs. independent from the input x i as depth goes to infinity. We call this effect "directional" gradient vanishing (see Section 3 for a more thorough explanation).

Interestingly this effect does not happen in batch-normalized networks, which yield -as we shall prove in Theorem 2 -a stable rank for any depth, thereby preserving a disentangled representation of the input and hence allowing the training of very deep networks. These results substantiate earlier empirical observations made by [START_REF] Bjorck | Understanding batch normalization[END_REF] for random BN-nets, and also validates the claim that BN helps with deep information propagation [START_REF] Samuel | Deep information propagation[END_REF].

Contributions

In summary, the work at hand makes the following two key contributions:

(i) We theoretically prove that BN indeed avoids rank collapse for deep linear neural nets under standard initialization and for any depth. In particular, we show that BN can be seen as a computationally cheap rank preservation operator, which may not yield hidden matrices with full rank but still preserves sufficient modes of variation in the data to achieve a scaling of the rank with Ω( √ d), where d is the width of the network. Subsequently, we leverage existing results from random matrix theory [START_REF] Bougerol | Products of Random Matrices with Applications to Schrödinger Operators[END_REF] to complete the picture with a simple proof of the above observed rank collapse for linear vanilla networks, which interestingly holds regardless of the presence of residual connections (Lemma 3). Finally, we connect the rank to difficulties in gradient based training of deep nets by showing that rank collapse makes the directional component of the gradients independent of the input.

(ii) We empirically show that the rank is indeed a crucial quantity for gradient-based learning. In particular, we show that both the rank and the final training accuracy quickly diminish in depth unless BN layers are incorporated in both simple feed-forward and convolutional neural nets. To take this reasoning beyond mere correlations, we actively intervene with the rank of networks before training and show that (a) one can break the training stability of BN by initializing in a way that reduces its rank-preserving properties, and (b) a rank-increasing pre-training procedure for vanilla networks can recover their training ability even for large depth. Interestingly, our pre-training method allows vanilla SGD to outperform BN on very deep MLPs. In all of our experiments, we find that SGD updates preserve the order of the initial rank throughout optimization, which underscores the importance of the rank at initialization for the entire convergence behavior.

Background and Preliminaries

Network description. We consider a given input X ∈ R d×N containing N samples in R d . Let 1 k ∈ R k denote the k-dimensional all one vector and H (γ) denote the hidden representation of X in layer of a BN-network with residual connections. The following recurrence summarizes the network mapping

H (γ) +1 = BN 0,1 d (H (γ) + γW H (γ) ), H (γ) 0 = X, (1)
where W ∈ R d×d and γ regulates the skip connection strength (in the limit, γ = ∞ recovers a network without skip connection) 4 . Throughout this work, we consider the network weights W to be initialized as follows.

Definition 1 (Standard weight initialization). The elements of weight matrices W are i.i.d. samples from a distribution P that has zero-mean, unit-variance, and its density is symmetric around zero 5 . We use the notation µ for the probability distribution of the weight matrices.

We define the BN operator BN α,β as in the original paper [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF], namely

BN α,β (H) = β • (diag (M (H))) -1 /2 H + α1 N , M (H) := 1 N HH , (2)
where • is a row-wise product. Both α ∈ R d and β ∈ R d are trainable parameters. Throughout this work we assume the initialization α = 0 and β = 1 d , and also omit corrections of the mean activity. As demonstrated empirically in Fig. 5, and theoretically in App. C this simplification does not change the performance of BN in our settings.

Rank notions. To circumvent numerical issues involved in rank computations we introduce a soft notion of the rank denoted by rank τ (H) (soft rank). Specifically, let σ 1 , . . . , σ d be the singular values of H. Then, given a τ > 0, we define rank τ (H) as

rank τ (H) = d i=1 1(σ 2 i /N ≥ τ ). (3)
Intuitively, rank τ (H) indicates the number of singular values whose absolute values are greater than √ N τ . It is clear that rank τ (H) is less or equal to rank(H) for all matrices H. For analysis purposes, we need an analytic measure of the collinearity of the columns and rows of H. Inspired by the so-called stable rank (see, e.g., [START_REF] Tropp | An introduction to matrix concentration inequalities[END_REF]), we thus introduce the following quantity

r(H) = Tr(M (H)) 2 / M (H) 2 F , M (H) = HH /N. (4)
In contrast to the algebraic rank, r(H) is differentiable with respect to H. Furthermore, the next lemma proves that the above quantity lower-bounds the soft-rank for the hidden representations.

Lemma 1. For an arbitrary matrix H ∈ R d×d , rank(H) ≥ r(H). For the sequence {H

(γ) } ∞ =1 defined in Eq. (2), rank τ (H (γ) ) ≥ (1 -τ ) 2 r(H (γ) ) holds for τ ∈ [0, 1].

Batch normalization provably prevents rank collapse

Since our empirical observations hold equally for both non-linear and linear networks, we here focus on improving the theoretical understanding in the linear case, which constitutes a growing area of research [START_REF] Saxe | Exact solutions to the nonlinear dynamics of learning in deep linear neural networks[END_REF][START_REF] Kawaguchi | Deep learning without poor local minima[END_REF][START_REF] Peter | Gradient descent with identity initialization efficiently learns positive-definite linear transformations by deep residual networks[END_REF][START_REF] Arora | A convergence analysis of gradient descent for deep linear neural networks[END_REF]. First, inspired by [START_REF] Dieuleveut | Bridging the gap between constant step size stochastic gradient descent and Markov chains[END_REF] and leveraging tools from Markov Chain theory, our main result proves that the rank of linear batch-normalized networks scales with their width as Ω( √ width). Secondly, we leverage results from random matrix theory [START_REF] Boucheron | Concentration inequalities: A nonasymptotic theory of independence[END_REF] to contrast our main result to unnormalized linear networks which we show to provably collapse to rank one, even in the presence of residual connections.

Main result

In the following we state our main result which proves that batch normalization indeed prevents the rank of all hidden layer activations from collapsing to one. Please see Appendix E for the more formal version of this theorem statement. Theorem 2.

[Informal] Suppose that the rank(X) = d and that the weights W are initialized in a standard i.i.d. zero-mean fashion (see Def. 1). Then, the following limits exist such that

lim L→∞ 1 L L =1 rank τ (H (γ) ) ≥ lim L→∞ (1 -τ ) 2 L L =1 r(H (γ) ) = Ω((1 -τ ) 2 √ d) (5)
holds almost surely for a sufficiently small γ (independent of ) and any τ ∈ [0, 1), under some additional technical assumptions. Please see Theorem 14 in the Appendix for the formal statement.

Theorem 2 yields a non trivial width-dependency. Namely, by setting for example τ := 1/2, the result states that the average number of singular values with absolute value greater than N/2 is at least Ω( √ d) on average. To put this into context: If one were to replace diag(M ) -1 /2 by the full inverse (M ) -1 /2 in Eq. ( 2), then BN would effectively constitute a classical whitening operation such that all {H (γ) } L =1 would be full rank (equal to d). However, as noted in the original BN paper [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF], whitening is obviously expensive to compute and furthermore prohibitively costly to incorporate in back-propagation. As such, BN can be seen as a computationally inexpensive approximation of whitening, which does not yield full rank hidden matrices but still preserves sufficient variation in the data to provide a rank scaling as Ω(

√ d).

Although the lower-bound in Thm. 2 is established on the average over infinite depth (i.e., L → ∞), Corollary 15 (in App. E) proves that the same bound holds for all rank(H ) and rank τ (H ).

Necessary assumptions. The above result relies on two key assumptions: (i) First, the input X needs to be full rank. (ii) Second, the weights have to be drawn according to the standard initialization scheme. We believe that both assumptions are indeed necessary for BN to yield a robust rank.

Regarding (i), we consider a high input rank a natural condition since linear neural nets cannot possibly increase the rank when propagating information through their layers. Of course, full rank is easily achieved by an appropriate data pre-processing. Yet, even when the matrix is close to low rank we find that BN is actually able to amplify small variations in the data (see Fig. 3.b). 6 Notably, we observed that hidden representations remain full rank if

H (γ) 0 is full-rank and N = O( √ d).
Regarding (ii), we derive -based on our theoretical insights -an adversarial initialization strategy that corrupts both the rank robustness and optimization performance of batch-normalized networks, thus suggesting that the success of BN indeed relies heavily on the standard i.i.d. zero-mean initialization.

Experimental validation. In order to underline the validity of Theorem 2 we run multiple simulations by feeding Gaussian data of dimensionality d = N into networks of growing size and with different residual strengths. For each network, we compute the mean and standard deviation of the soft rank rank τ with τ = 0.5. As depicted in Fig. 3, the curves clearly indicate a Ω( √ d) dependency for lim L→∞ L =1 rank τ (H )/L, just as predicted in the Theorem. Although the established guarantee requires the weight on the parametric branch (i.e., γ) to be small, the results of Fig. 3 indicate that the established lower bound holds for a much wider range including the case where no residual connections are used at all (γ = ∞). (γ) for increasing values of given nearly collinear inputs. As can be seen, BN quickly amplifies smaller variations in the data while reducing the largest one.

Comparison with unnormalized networks

In order to stress the importance of the above result, we now compare the predicted rank of H with the rank of unnormalized linear networks, which essentially constitute a linear mapping in the form of a product of random matrices. The spectral distribution of products of random matrices with i.i.d. standard Gaussian elements has been studied extensively [START_REF] Bjorck | Understanding batch normalization[END_REF][START_REF] Peter | Lyapunov exponents for products of complex Gaussian random matrices[END_REF][START_REF] Liu | Bulk and soft-edge universality for singular values of products of ginibre random matrices[END_REF]. Interestingly, one can show that the gap between the top and the second largest singular value increases with the number of products (i.e., ) at an exponential rate7 [START_REF] Peter | Lyapunov exponents for products of complex Gaussian random matrices[END_REF][START_REF] Liu | Bulk and soft-edge universality for singular values of products of ginibre random matrices[END_REF]. Hence, the matrix converges to a rank one matrix after normalizing by the norm. In the following, we extend this result to products of random matrices with a residual branch that is obtained by adding the identity matrices. Particularly, we consider the hidden states H of the following linear residual network:

H = B X, B := k=1 (I + γW k ). ( 6
)
Since the norm of H is not necessarily bounded, we normalize as H = B X/ B . The next lemma characterizes the limit behaviour of { H }. Lemma 3. Suppose that γ ∈ (0, 1) and assume the weights W to be initialized as in Def. 1 with element-wise distribution P. Then we have for linear networks, which follow recursion [START_REF] Peter | Gradient descent with identity initialization efficiently learns positive-definite linear transformations by deep residual networks[END_REF], that:

a. If P is standard Gaussian, then the sequence { H } converges to a rank one matrix.

b.

If P is uniform[- √ 3, √ 3
], then there exists a monotonically increasing sequence of integers 1 < 2 , . . . such that the sequence { H k } converges to a rank one matrix.

This results stands in striking contrast to the result of Theorem 2 established for batch-normalized networks. 8 Interestingly, even residual skip connections cannot avoid rank collapse for very deep neural networks, unless one is willing to incorporate a depth dependent down-scaling of the parametric branch as for example done in [START_REF] Allen-Zhu | A convergence theory for deep learning via over-parameterization[END_REF], who set γ = O( 1 L ) . Remarkably, Theorem 2 shows that BN layers provably avoid rank collapse without requiring the networks to become closer and closer to identity. Remarkably, the remaining direction after rank collapse depends exclusively on the random weights and it is independent of the input.

Implications of rank collapse on gradient based learning. In order to explain the severe consequence of rank collapse on optimization performance reported in Fig. 1, we study the effect of rank one hidden-layer representations on the gradient of the training loss for distinct input samples. Let L i denote the training loss for datapoint i on a vanilla network as in Eq. ( 6). Furthermore, let the final classification layer be parametrized by W L+1 ∈ R dout×d . Then, given that the hidden presentation at the last hidden layer L is rank one, the normalized gradients of the loss with respect to weights of individual neurons k ∈ 1, ..., d out in the classification layer (denoted by ∇ W L+1,k L i , where ∇ W L+1,k L i = 1) are collinear for any two datapoints i and j, i.e. ∇ W L+1,k L i = ∓∇ W L+1,k L j . A formal statement is presented in Prop. 19 in the Appendix alongside empirical validations on a VGG19 network (Fig. 10). This result implies that the commonly accepted vanishing gradient norm hypothesis is not descriptive enough since SGD does not take small steps into the right direction, but into a random direction that is independent from the input. In other words, deep neural networks are prone to directional gradient vanishing after initialization, which is caused by the collapse of the last hidden layer activations to a very small subspace (one line in R d in the extreme case of rank one activations).

The important role of the rank

The preceding sections highlight that the rank of the hidden representations is a key difference between random vanilla and BN networks. We now provide three experimental findings that substantiate the particular importance of the rank at the beginning of training: First, we find that an unsupervised, rank-increasing pre-training allows SGD on vanilla networks to outperform BN networks. Second, we show that the performance of BN-networks is closely tied to a high rank at initialization. Third, we report that SGD updates preserve the initial rank magnitude throughout the optimization process.

Outperforming BN using a pre-training step. As discussed above, batch normalization layers are very effective at avoiding rank collapse. Yet, this is of course not the only way to preserve rank. Based upon our theoretical insights, we leverage the lower bound established in Eq. ( 4) to design a pre-training step that not only avoids rank collapse but also accelerates the convergence of SGD. Our proposed procedure is both simple and computationally cheap. Specifically, we maximize the lower-bound r(H ) (in Eq. ( 4)) on the rank of the hidden presentation H in each layer . Since this function is differentiable with respect to its input, it can be optimized sufficiently by just a few steps of (stochastic) gradient ascent (see Section G in the Appendix for more details). As can be seen, the slow down in depth is much less severe for the pre-trained networks. This improvement is, also, reflected both in terms of training accuracy and test loss (see Fig. 11 in Appendix). Interestingly, the pre-training is not only faster than BN on deep networks, but it is also straight-forward to use in settings where the application of BN is rather cumbersome such as for very small batch sizes or on unseen data [START_REF] Ioffe | Batch renormalization: Towards reducing minibatch dependence in batchnormalized models[END_REF][START_REF] Wu | Group normalization[END_REF].

Breaking batch normalization. Some scholars hypothesize that the effectiveness of BN stems from a global landscape smoothing [START_REF] Santurkar | How does batch normalization help optimization?[END_REF] or a certain learning rate tuning [START_REF] Arora | Theoretical analysis of auto rate-tuning by batch normalization[END_REF], that are thought to be induced by the normalization. Under these hypotheses, one would expect that SGD converges fast on BN-nets regardless of the initialization. Yet, we here show that the way that networks are initialized does play a crucial role for the subsequent optimization performance of BN-nets. Particularly, we train two MLPs with batchnorm, but change the initialization for the second net from the standard PyTorch way W l,i,j ∼ uniform -1/ √ d l , 1/ √ d l [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF][START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF] to W l,i,j ∼ uniform 0, +2/ √ d l , where d l is the layer size. As can be seen to the right, this small change reduces the rank preserving quality of BN significantly, which is reflected in much slower learning behaviour. Even sophisticated modern day architectures such as VGG and ResNet networks are unable to fit the CIFAR-10 dataset after changing the initialization in this way (see Fig. 5).

Rank through the optimization process. The theoretical result of Theorem 2 considers the rank at random initialization. To conclude, we perform two further experiments which confirm that the initial rank strongly influences the speed of SGD throughout the entire optimization process. In this regard, Fig. 6 reports that SGD preserves the initial magnitude of the rank to a large extent, regardless of the specific network type. This is particularly obvious when comparing the two BN initializations. A further noteworthy aspect is the clear correlation between the level of pre-training and optimization performance on vanilla nets. Interestingly, this result does again not only hold on simple MLPs but also generalizes to modern day networks such as the VGG-19 (see Fig. 5) and ResNet50 architecture (see Appendix I).

Discussions

In summary, our work highlights a key difference between random vanilla-and BN networks. While the rank of the hidden representations quickly collapses to one as the depth of vanilla networks increases, BN is robust against such rank collapse. This intriguing property arises due to the standard initialization of weights and also it is preserved through the optimization process. Notably, our theoretical analysis proves this striking difference for linear MLPs and holds empirically across a wide range of data sets and network architectures. Our experiments further highlight the determining role of the rank quantity in the training speed. Inspired by these observations, we develop a novel pretraining method that allows previously un-trainable very deep vanilla networks to learn, sometimes even faster than BN-MLPs of the same size. Thereby our study reveals a key requirement for a proper initialization of deep neural networks, opening doors to the development of effective initialization schemes for modern-day architectures.

We thus consider our work a relevant step towards a better understanding of optimization for deep neural networks. Furthermore, our findings give rise to several interesting follow-up questions: (i) Can one generalize the analysis of Theorem 2 to ReLU and other non-linear nets to prove the observed rank robustness (e.g. Fig. 2)? (ii) is it possible to rigorously prove that SGD updates preserve the rank magnitude throughout optimization, as observed in Fig. 6)? (iii) Is it possible to use the develop a similarly effective pre-training for convolution and recurrent networks? (iv) How can one theoretically characterize the connection between the convergence of SGD and the rank quantity (a follow-up on directional gradient vanishing)? (v) Does rank robustness explain the success of related architectures such as layer normalization [START_REF] Lei Ba | Layer normalization[END_REF], weight normalization [START_REF] Salimans | Weight normalization: A simple reparameterization to accelerate training of deep neural networks[END_REF]) and modern initialization techniques such as fix up initialization [START_REF] Zhang | Fixup initialization: Residual learning without normalization[END_REF]? We believe that these questions will spark an interesting line of future research towards the goal of fully understanding optimization in deep neural networks.

Broader impact

As we only contribute to a better understanding of neural network training in general, we consider our work fundamental research without any specific application. Hence a broader impact discussion is not applicable.

Figure 1 :

 1 Figure 1: Effect of depth on rank and learning, on the Fashion-MNIST dataset with ReLU multilayer

Figure 2 :

 2 Figure 2: Rank comparison of last hidden activation: Log(rank) of the last hidden layer's activation over total number of layers (blue for BN-and orange for vanilla-networks) for Gaussian inputs. Networks are MLPs of width d = 32. (Left) Linear activations, (Right) ReLU activations. Mean and 95% confidence interval of 10 independent runs. While the rank quickly drops in depth for both networks, BN stabilizes the rank above √ d.

Figure 3 :

 3 Figure 3: a) Result of Theorem 2 for different values of γ, where γ = ∞ depicts networks without skip connections. Each point is the average rank1 /2 over depth (L = 10 6 ) of nets of width d ∈ {8, 16, .., 256} an on x-axis. b) Top 10 singular values of H(γ) for increasing values of given nearly collinear inputs. As can be seen, BN quickly amplifies smaller variations in the data while reducing the largest one.

Figure 4 :

 4 Figure 4: Pre-training versus BN: Loss over epochs on CIFAR-10 for MLPs of increasing depth with 128 hidden units and ReLU activation. Trained with SGD (batchsize 64) and grid-searched stepsize. See Fig. 11 for the corresponding test loss and accuracy as well as Fig. 12 for FashionMNIST results.

Fig. 4

 4 Fig.4compares the convergence rate of SGD on pre-trained vanilla networks and BN-networks. As can be seen, the slow down in depth is much less severe for the pre-trained networks. This improvement is, also, reflected both in terms of training accuracy and test loss (see Fig.11in Appendix). Interestingly, the pre-training is not only faster than BN on deep networks, but it is also straight-forward to use in settings where the application of BN is rather cumbersome such as for very small batch sizes or on unseen data[START_REF] Ioffe | Batch renormalization: Towards reducing minibatch dependence in batchnormalized models[END_REF][START_REF] Wu | Group normalization[END_REF].

  in last hidden layer W U[ a, a] W U[0, 2a] W U[ a, a] no centering

Figure 5 :

 5 Figure 5: Breaking Batchnorm: CIFAR-10 on VGG19 with standard PyTorch initialization as well as a uniform initialization of same variance. (Left) training accuracy, (Right) Rank of last hidden layer computed using torch.matrix_rank(). Plot also shows results for standard initialization and BN without mean deduction. Avg. and 95% CI of 5 independent runs. (See Fig. 13 in Appendix for similar results on ResNet-50).

Figure 6 :

 6 Figure 6: Pretraining: Fashion-MNIST on MLPs of depth 32 and width 128. (Left) Training accuracy, (Right) Lower bound on rank. Blue line is a ReLU network with standard initialization. Other solid lines are pre-trained layer-wise with 25 (orange) and 75 (green) iterations to increase the rank. Dashed lines are batchnorm networks with standard and asymmetric initialization. Average and 95% confidence interval of 5 independent runs.

Computed using torch.matrix_rank(), which regards singular values below σmax × d × 10 -7 as zero. This is consistent with both Matlab and Numpy.

For the sake of simplicity, we here assume that the numbers of hidden units is equal across layers. In App. E we show how our results extend to nets with varying numbers of hidden units.

Two popular choices for P are the Gaussian distribution N (0, 1) and the uniform distribution U([-1, 1]). The variance can be scaled with the choice of γ to match the prominent initializations from[START_REF] He | Delving deep into rectifiers: Surpassing human-level performance on imagenet classification[END_REF] and[START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF]. Note that the symmetry implies that the law of each element [W ]ij equates the law of -[W ]ij.

Intuitively this means that even if two data points are very close to each other in the input space, their hidden presentation can still be disentangled in batch-normalized networks (see Appendix E for more details)

The growth-rate of the i-th singular value is determined by the i-th Lyapunov exponent of the product of random matrices. We refer the reader to[START_REF] Peter | Lyapunov exponents for products of complex Gaussian random matrices[END_REF] for more details on Lyapunov exponents.

According to the observations in Fig.2, the result of part b holds for the usual sequence of indices { k = k}, which indicates that { H k } converges to a rank one matrix even in the case of uniform initialization.

Acknowledgement

We thank Dr. Aran Raoufi, Dr. Olivier Ledoit, Gary Becigneul, and Olivier Ledoit for their helpful discussions.

Funding

This project is fully supported by ETH-Zurich fellowships.