Introduction
Consider the problem of minimizing a smooth function L : R d → R L: \mathbb{R}^d \to\mathbb{R} L : R d → R x 0 ∈ R d x_0 \in\mathbb{R}^d x 0 ∈ R d h > 0 h > 0 h > 0
x n + 1 = x n − h ∇ L ( x n ) , n ≥ 0. (1) x_{n+1}=x_n-h\nabla L(x_n), \qquad n\geq0. \tag*{(1)} x n + 1 = x n − h ∇ L ( x n ) , n ≥ 0. (1) Since (1) is an Euler iteration for the ordinary differential equation x ′ = − ∇ L ( x ) x'=-\nabla L(x) x ′ = − ∇ L ( x )
d d t L ( x ( t ) ) = − ∥ ∇ L ( x ( t ) ) ∥ 2 , \frac{d}{dt}L(x(t))=-\lVert\nabla L(x(t))\rVert^2, d t d L ( x ( t )) = − ∥ ∇ L ( x ( t )) ∥ 2 , a Taylor expansion shows that, at any noncritical point, there exists a sufficiently small h h h L L L ∇ L \nabla L ∇ L { x : L ( x ) ≤ L ( x 0 ) } \{x:L(x)\leq L(x_0)\} { x : L ( x ) ≤ L ( x 0 )} h > 0 h>0 h > 0 L L L { x n } n ≥ 0 \{x_n\}_{n\geq0} { x n } n ≥ 0 L L L h h h § 2.1 below).
The algorithm is simple, involving only first partial derivatives and elementary vector operations. It scales well with dimension as, apart from the cost of evaluating the gradient, memory usage and arithmetic operations are linear in the number of variables. Moreover, convergence of the algorithm in the small h h h
This being said, gradient descent has well-known shortcomings. Applying a single learning rate to all components of the gradient leads to poor performance in problems where L L L
These issues lead naturally to interest in adaptive step modifications like AdaGrad [15 ] , RMSProp [44 ] , Adam [30 ] , and Muon [25 ] . RMSProp is a commonly used modification of gradient descent and is the main focus of the present work. The algorithm is defined formally in § 3.1 , but the idea is to adapt to the geometry of L L L x ∈ R d x \in\mathbb{R}^{d} x ∈ R d ϵ \epsilon ϵ § 10 ). See § 2.3 for some additional context.
From a dynamical systems perspective, the starting point for the discussion of RMSProp is the convergence result of Dimitrieski et al. [14 Theorem 1] . The result is restated precisely as Theorem 4.4 below but, briefly, it says that if L L L 0 < h ∗ 0 < h_{*} 0 < h ∗ 0 < h < h ∗ 0 < h < h_{*} 0 < h < h ∗ L L L h ∗ > 0 h_{*} > 0 h ∗ > 0 L L L
The present work studies the dynamics of RMSProp for the quadratic loss functional L ( x ) = 1 2 x T A x L(x) = \frac{1}{2}x^{\mathsf{T}}Ax L ( x ) = 2 1 x T A x A A A § 2.4 . However, we note that the stability theorem of [14 ] gives explicitly a critical learning rate of
h ∗ = 2 ϵ λ max ( A ) (2) h_{*} = \frac{2\epsilon}{\lambda_{\max}(A)} \tag*{(2)} h ∗ = λ m a x ( A ) 2 ϵ (2) for the quadratic case. Here λ max ( A ) \lambda_{\max}(A) λ m a x ( A ) A A A ϵ \epsilon ϵ § 3.1 for the precise definition of the map.) For quadratic loss functions, we establish two global stability/instability results which complement the existing small h h h [14 ] .
First, we show that for all h > h ∗ h > h_{*} h > h ∗ § 3.1 , the version of RMSProp we work with is not C 1 C^{1} C 1
Nevertheless, we show that for every h > h ∗ h > h_{*} h > h ∗ h h h h h h h 2 h^{2} h 2
Taken together, these results lead to a practical observation about quadratic RMSProp, and constitute one of the main results of the present work. That is, while for h > h ∗ h > h_{*} h > h ∗ L L L h h h This has practical applications, for example in scheduling adjustments to the size of the learning rate h h h § 9.1 .
A natural next question is “what kinds of dynamics actually occur in the compact trapping region?” We prove that, for generic A A A h > h ∗ h > h_{*} h > h ∗ h ∗ h_{*} h ∗ h = h ∗ h = h_{*} h = h ∗ C 1 C^{1} C 1 2 d − 1 2^{d-1} 2 d − 1
The results described above hold for a generic symmetric positive definite matrix A A A d d d A A A 2 × 2 2 \times2 2 × 2 5 × 5 5 \times5 5 × 5 § 6 and § 7 , but a brief summary is as follows: we often find Neimark–Sacker bifurcations of the family of h > h ∗ h > h^{*} h > h ∗ § 2.5 ) leads sometimes to further bifurcations into invariant tori, and sometimes to period doubling and the appearance of strange chaotic attractors. We also perform a limited parameter space search, classifying parameters by the kind of attractor they lead to.
Finally, the relevance of this study is not limited to RMSProp. For example, with its time-dependent bias correction removed, Adam reduces to RMSProp when β 1 = 0 \beta_{1} = 0 β 1 = 0 [23 ] , we continue some of our results into Adam, suggesting that it may be amenable to the same analysis. We discuss this continuation further in § 10 .
The remainder of the paper is organized as follows. § 2 discusses related work and motivation. § 3 defines the RMSProp map, § 4 contains the stability and trapping results, and § 5 treats period-two orbits. § 6 and § 7 describe the numerical methods and the numerical bifurcation analysis, § 8 surveys the other parameters, and § 10 treats autonomous Adam.
Related work and motivation
Machine learning context
Gradient descent and its variants are enjoying a period of substantial renewed interest, due at least in part to applications in machine learning (ML) like supervised and unsupervised training of neural networks and large language models. A thorough review of this literature is a task far beyond the scope of the present work, and we refer the interested reader to the review article of LeCun, Bengio, and Hinton [35 ] , and to the notes of Ruder [42 ] for much more exhaustive discussion. We note that referring to L L L h h h
Dynamics of discretized gradient descent
Treating the gradient descent algorithm as a dynamical system with the learning rate as a bifurcation parameter is a classic topic in numerical analysis and dynamics of algorithms. That increasing the learning rate to accelerate convergence leads to oscillations is well known (see for example the book of Hairer and Wanner [18 ] ). It is also common to study gradient descent for quadratic loss functions, in which case the algorithm involves iterating a linear system, and all dynamics are understood by studying the eigenvalues of the resulting matrix. In this case it’s easy to see that convergence to the minimum, oscillations, and unbounded orbits are all possible. Beyond quadratic minimization, the gradient descent algorithm admits nonlinearities and the situation is much more complicated. Kong and Tao have proved that, for multiscale objectives, the algorithm exhibits chaotic dynamics when the learning rate is not small enough, providing a mechanism for escaping local minima as an alternative to noisy escape due to stochastic gradients [31 ] . We refer also to the work of Liang and Montúfar [36 ] , who prove the existence of chaotic behavior and fractal convergence boundaries when gradient descent is applied to certain matrix factorization problems.
Dynamics of RMSProp and the adaptive edge of stability
Violent but bounded oscillations occur in adaptive methods as well as in gradient descent. Their behavior, however, depends on more than the learning rate and the local curvature. In RMSProp, we also divide by a quantity determined by the memory variables. When gradients become large, the memory can grow and reduce subsequent updates. Thus, fixing h h h ϵ \epsilon ϵ [8 , 9 ] study this behavior in neural-network training. They account for the optimizer’s coordinate rescaling when measuring the largest eigenvalue of the Hessian, obtaining what they call the preconditioned sharpness. In their experiments, this quantity often stays near the stability boundary of the update with the rescaling held fixed. They call this the adaptive edge of stability. Moreover, oscillations at the edge of stability can be chaotic [8 ] .
A natural next question is “what kinds of dynamics actually occur when the iterates remain bounded?” For our quadratic loss, the Hessian is constant, so we focus on the effects of adaptation without also having to account for changes in the loss curvature. We note that the adaptive edge of stability concerns the changing preconditioned curvature along a training trajectory, while h ∗ h_* h ∗ (2) corresponds to a fixed stability boundary of the quadratic minimizer. Bock and Weiß [3 ] give local stability conditions for several versions of Adam and RMSProp, while Dimitrieski et al. [14 ] prove global convergence for epsilon-outside RMSProp below 2 ϵ / L 2\epsilon/L 2 ϵ / L L L L
Quadratic loss functionals
Restricting to a quadratic loss corresponds to isolating the quadratic term in the Taylor expansion of L L L L L L L L L
Beyond this, Meterez et al. [37 ] have argued in defense of the quadratic model, showing that local quadratic models predict well the training results of a language model with millions of parameters. The case of quadratic loss also illuminates the boundary between gradient descent and RMSProp. This is because gradient descent for a quadratic functional results in a linear dynamical system, while RMSProp remains nonlinear due to the fact that the step sizes depend on past gradients. Complicated behavior of RMSProp in a quadratic example is then a feature of the optimizer, and not the chosen problem. Beyond this, the quadratic case admits detailed analytical study, serving as a base for further analysis of nonconvex, noisy, or multiscale losses.
The parameterization method
The parameterization method is a functional analytic framework for studying invariant manifolds, first proposed by Cabré, Fontich, and de la Llave for studying spectral submanifolds associated with fixed points of maps between Banach spaces [5 , 6 , 7 ] , and later extended to quasi-periodic orbits of maps, and their attached invariant stable/unstable manifolds by Haro and de la Llave [20 , 21 , 17 ] . By now the parameterization method supports a thriving sub-industry, and a thorough review of the literature is a task beyond the scope of the present work. We refer the interested reader to the recent book of Haro et al. [19 ] . We describe our computation and continuation of invariant circles and tori in § 6.2 .
RMSProp as a dynamical system
Throughout, d d d ⊘ \oslash ⊘ ⊙ \odot ⊙ x ∘ r x^{\circ r} x ∘ r x ∘ \sqrt[\circ]{x} ∘ x x x x I d I_d I d d × d d \times d d × d 1 ∈ R d \mathbf{1} \in\mathbb{R}^d 1 ∈ R d d d d A A A d × d d \times d d × d a i j a_{ij} a ij ∣ A ∣ |A| ∣ A ∣ ∣ a i j ∣ |a_{ij}| ∣ a ij ∣ ∣ x ∣ = ( ∣ x 1 ∣ , … , ∣ x d ∣ ) T |x| = (|x_1|,\ldots,|x_d|)^\mathsf{T} ∣ x ∣ = ( ∣ x 1 ∣ , … , ∣ x d ∣ ) T
Algorithm definition
Given a learning rate h > 0 h > 0 h > 0 β ∈ [ 0 , 1 ) \beta\in[0,1) β ∈ [ 0 , 1 ) [41 ] .
We note that the definition is not universal, as different optimizer implementations place the numerical stability constant ϵ \epsilon ϵ ϵ \epsilon ϵ [41 ] . The current Keras/TensorFlow RMSProp implementation instead adds ϵ \epsilon ϵ v k + 1 + ϵ ∘ \sqrt[\circ]{v_{k+1}+\epsilon} ∘ v k + 1 + ϵ [29 ] . Optax, an optimization library for JAX, supports both conventions and defaults to epsilon inside the square root [39 ] .
Both epsilon conventions preserve the phase space R d × R ≥ 0 d \mathbb{R}^d \times\mathbb{R}^d_{\ge0} R d × R ≥ 0 d F F F (5) Algorithm 1 RMSProp Update Rule (Epsilon-Outside Variant): Tieleman and Hinton [44 ]
Require: Initial spatial vector x 0 ∈ R d x_0 \in\mathbb{R}^d x 0 ∈ R d v 0 ∈ R ≥ 0 d v_0 \in\mathbb{R}_{\geq0}^d v 0 ∈ R ≥ 0 d v 0 = 0 v_0 = 0 v 0 = 0 ϵ ∈ R > 0 \epsilon\in\mathbb{R}_{>0} ϵ ∈ R > 0
Require: Objective function L : R d → R L : \mathbb{R}^d \to\mathbb{R} L : R d → R h > 0 h > 0 h > 0 β ∈ [ 0 , 1 ) \beta\in[0,1) β ∈ [ 0 , 1 )
Require: Number of iterations N N N
for k = 0 , 1 , … , N − 1 k = 0,1,\ldots,N-1 k = 0 , 1 , … , N − 1
g k ← ∇ L ( x k ) g_k \gets\nabla L(x_k) g k ← ∇ L ( x k )
v k + 1 ← β v k + ( 1 − β ) ( g k ) ⊙ 2 v_{k+1} \gets\beta v_k + (1-\beta)(g_k)^{\odot2} v k + 1 ← β v k + ( 1 − β ) ( g k ) ⊙ 2
x k + 1 ← x k − h ( g k ) ⊘ v k + 1 + ϵ 1 ⊙ x_{k+1} \gets x_k - h(g_k) \oslash\sqrt[\odot]{v_{k+1}+\epsilon1} x k + 1 ← x k − h ( g k ) ⊘ ⊙ v k + 1 + ϵ 1
end for
return x N , v N x_N,v_N x N , v N
below, in the phase space R d × R ≥ 0 d \mathbb{R}^d \times\mathbb{R}_{\geq0}^d R d × R ≥ 0 d v ≥ 0 v \geq0 v ≥ 0 C 1 C^1 C 1
The epsilon-inside map additionally extends smoothly to an open neighborhood of the origin that includes negative memory coordinates. We focus on epsilon outside because the loss of a C 1 C^1 C 1
As already mentioned in the introduction, we focus on the performance and properties of RMSProp with L L L
L ( x ) = 1 2 x T A x , (3) L(x) = \frac{1}{2}x^{\mathsf{T}}Ax, \tag*{(3)} L ( x ) = 2 1 x T A x , (3) where A = A T ≻ 0 A = A^{\mathsf{T}} \succ0 A = A T ≻ 0 h , ϵ > 0 h,\epsilon> 0 h , ϵ > 0 0 < β < 1 0 < \beta< 1 0 < β < 1 v 0 ≥ 0 v_0 \geq0 v 0 ≥ 0
The map
Combining the memory and spatial updates defines a map on R d × R ≥ 0 d \mathbb{R}^d \times\mathbb{R}_{\geq0}^d R d × R ≥ 0 d (3) , set
F V ( x , v ) = β v + ( 1 − β ) ( A x ) ⊙ 2 . (4) F_V(x,v) = \beta v + (1-\beta)(Ax)^{\odot2}. \tag*{(4)} F V ( x , v ) = β v + ( 1 − β ) ( A x ) ⊙ 2 . (4) The epsilon-outside RMSProp map F : R d × R ≥ 0 d → R d × R ≥ 0 d F : \mathbb{R}^d \times\mathbb{R}_{\geq0}^d \to\mathbb{R}^d \times\mathbb{R}_{\geq0}^d F : R d × R ≥ 0 d → R d × R ≥ 0 d
F ( x , v ) = [ x − h A x ⊘ β v + ( 1 − β ) ( A x ) ⊙ 2 + ϵ 1 ⊙ β v + ( 1 − β ) ( A x ) ⊙ 2 ] = [ x − h A x ⊘ F V ( x , v ) + ϵ 1 ⊙ F V ( x , v ) ] (5) F(x,v) =
\begin{bmatrix}
x-hAx \oslash\sqrt[\odot]{\beta v+(1-\beta)(Ax)^{\odot2}+\epsilon1}\\
\beta v+(1-\beta)(Ax)^{\odot2}
\end{bmatrix}
=
\begin{bmatrix}
x-hAx \oslash\sqrt[\odot]{F_V(x,v)+\epsilon1}\\
F_V(x,v)
\end{bmatrix} \tag*{(5)} F ( x , v ) = [ x − h A x ⊘ ⊙ β v + ( 1 − β ) ( A x ) ⊙ 2 + ϵ 1 β v + ( 1 − β ) ( A x ) ⊙ 2 ] = [ x − h A x ⊘ ⊙ F V ( x , v ) + ϵ 1 F V ( x , v ) ] (5) For fixed parameters, write F ( x , v ) = ( F 1 ( x , v ) , F 2 ( x , v ) ) F(x,v) = (F_1(x,v),F_2(x,v)) F ( x , v ) = ( F 1 ( x , v ) , F 2 ( x , v )) F 1 F_1 F 1 F 2 = F V F_2 = F_V F 2 = F V F I : R d × R ≥ 0 d → R d × R ≥ 0 d F_I : \mathbb{R}^d \times\mathbb{R}_{\geq0}^d \to\mathbb{R}^d \times\mathbb{R}_{\geq0}^d F I : R d × R ≥ 0 d → R d × R ≥ 0 d
F I ( x , v ) = [ x − h A x ⊘ F V ( x , v ) + ϵ 1 ⊙ β v + ( 1 − β ) ( A x ) ⊙ 2 ] F_I(x,v) =
\begin{bmatrix}
x-hAx \oslash\sqrt[\odot]{F_V(x,v)+\epsilon1}\\
\beta v+(1-\beta)(Ax)^{\odot2}
\end{bmatrix} F I ( x , v ) = [ x − h A x ⊘ ⊙ F V ( x , v ) + ϵ 1 β v + ( 1 − β ) ( A x ) ⊙ 2 ] Origin stability and global trapping
Dimitrieski et al. [14 ] establish global asymptotic stability for epsilon-outside RMSProp, as defined in (5), below 2 ϵ / L 2\epsilon/L 2 ϵ / L L L L
Definition 4.1 (Lyapunov stability and instability [16 ] ). Let F : R d × R ≥ 0 d → R d × R ≥ 0 d F : \mathbb{R}^d \times\mathbb{R}_{\geq0}^d \to\mathbb{R}^d \times\mathbb{R}_{\geq0}^d F : R d × R ≥ 0 d → R d × R ≥ 0 d z ∗ z_\ast z ∗ z ∗ z_\ast z ∗ stable in the sense of Lyapunov if, for every ρ > 0 \rho> 0 ρ > 0 δ > 0 \delta> 0 δ > 0
∥ z 0 − z ∗ ∥ < δ ⟹ ∥ F k ( z 0 ) − z ∗ ∥ < ρ for every k ≥ 0 \lVert z_0-z_\ast\rVert< \delta\quad\Longrightarrow\quad\lVert F^k(z_0)-z_\ast\rVert< \rho\quad\text{for every } k \geq0 ∥ z 0 − z ∗ ∥ < δ ⟹ ∥ F k ( z 0 ) − z ∗ ∥ < ρ for every k ≥ 0 for all z 0 ∈ R d × R ≥ 0 d z_0 \in\mathbb{R}^d \times\mathbb{R}_{\geq0}^d z 0 ∈ R d × R ≥ 0 d Lyapunov unstable if it is not Lyapunov stable; equivalently, there exists ρ 0 > 0 \rho_0 > 0 ρ 0 > 0 δ > 0 \delta> 0 δ > 0 z 0 ∈ R d × R ≥ 0 d z_0 \in\mathbb{R}^d \times\mathbb{R}_{\geq0}^d z 0 ∈ R d × R ≥ 0 d k ≥ 0 k \geq0 k ≥ 0
∥ z 0 − z ∗ ∥ < δ , ∥ F k ( z 0 ) − z ∗ ∥ ≥ ρ 0 . \lVert z_0-z_\ast\rVert< \delta,\qquad\lVert F^k(z_0)-z_\ast\rVert\geq\rho_0. ∥ z 0 − z ∗ ∥ < δ , ∥ F k ( z 0 ) − z ∗ ∥ ≥ ρ 0 . Global asymptotic stability below the threshold
Definition 4.2 (μ \mu μ f : R d → R f:\mathbb{R}^d \to\mathbb{R} f : R d → R μ \mu μ μ \mu μ μ > 0 \mu> 0 μ > 0
f ( y ) ≥ f ( x ) + ∇ f ( x ) T ( y − x ) + μ 2 ∥ y − x ∥ 2 2 f(y) \geq f(x) + \nabla f(x)^{\mathsf{T}}(y-x) + \frac{\mu}{2}\lVert y-x\rVert_2^2 f ( y ) ≥ f ( x ) + ∇ f ( x ) T ( y − x ) + 2 μ ∥ y − x ∥ 2 2 for all x , y ∈ R d x,y \in\mathbb{R}^d x , y ∈ R d
Definition 4.3 (L L L f : R d → R f:\mathbb{R}^d \to\mathbb{R} f : R d → R L L L L > 0 L > 0 L > 0 L L L
∥ ∇ f ( y ) − ∇ f ( x ) ∥ 2 ≤ L ∥ y − x ∥ 2 \lVert\nabla f(y)-\nabla f(x)\rVert_2 \leq L\lVert y-x\rVert_2 ∥ ∇ f ( y ) − ∇ f ( x ) ∥ 2 ≤ L ∥ y − x ∥ 2 for all x , y ∈ R d x,y \in\mathbb{R}^d x , y ∈ R d
We recall, for reference, the constant-step stability result, whose statement and proof are found in [14 Thm. 1] .
Theorem 4.4 (Global stability of epsilon-outside RMSProp). Let f : R d → R f:\mathbb{R}^d \to\mathbb{R} f : R d → R μ \mu μ μ > 0 \mu> 0 μ > 0 L L L L > 0 L > 0 L > 0 x ∗ x_\ast x ∗
v k + 1 = β v k + ( 1 − β ) ( ∇ f ( x k ) ) ⊙ 2 , x k + 1 = x k − h ∇ f ( x k ) ⊘ ( v k + 1 + ϵ 1 ) , v_{k+1}=\beta v_k+(1-\beta)(\nabla f(x_k))^{\odot2},\qquad x_{k+1}=x_k-h\nabla f(x_k)\oslash\left(\sqrt{\smash[b]{v_{k+1}}}+\epsilon\mathbf{1}\right), v k + 1 = β v k + ( 1 − β ) ( ∇ f ( x k ) ) ⊙ 2 , x k + 1 = x k − h ∇ f ( x k ) ⊘ ( v k + 1 + ϵ 1 ) , with x 0 ∈ R d x_0 \in\mathbb{R}^d x 0 ∈ R d v 0 ∈ R ≥ 0 d v_0 \in\mathbb{R}_{\geq0}^d v 0 ∈ R ≥ 0 d 0 < β < 1 0 < \beta< 1 0 < β < 1 ϵ > 0 \epsilon> 0 ϵ > 0
0 < h < 2 ϵ L . 0<h<\frac{2\epsilon}{L}. 0 < h < L 2 ϵ . Then ( x ∗ , 0 ) (x_\ast,0) ( x ∗ , 0 ) R d × R ≥ 0 d \mathbb{R}^d \times\mathbb{R}_{\geq0}^d R d × R ≥ 0 d
For the quadratic objective L ( x ) = 1 2 x T A x L(x)=\frac{1}{2}x^{\mathsf{T}}Ax L ( x ) = 2 1 x T A x x ∗ = 0 x_\ast=0 x ∗ = 0 μ = λ min ( A ) \mu=\lambda_{\min}(A) μ = λ m i n ( A ) L = L A L=L_A L = L A L A = λ max ( A ) = ∥ A ∥ 2 L_A=\lambda_{\max}(A)=\lVert A\rVert_2 L A = λ m a x ( A ) = ∥ A ∥ 2
Corollary 4.5 (Global asymptotic stability). If 0 < h < h ∗ 0<h<h_\ast 0 < h < h ∗ L ( x ) = 1 2 x T A x L(x)=\frac{1}{2}x^{\mathsf{T}}Ax L ( x ) = 2 1 x T A x R d × R ≥ 0 d \mathbb{R}^d \times\mathbb{R}_{\geq0}^d R d × R ≥ 0 d F F F
In the instability proof below, we use the component maps F 1 ( x , v ) F_1(x,v) F 1 ( x , v ) F 2 ( x , v ) F_2(x,v) F 2 ( x , v )
Instability above the threshold
To prove Lyapunov instability for h > h ∗ h>h_\ast h > h ∗
We begin by observing that the RMSProp mapping is Fréchet differentiable at the origin, relative to the phase space v ≥ 0 v\geq0 v ≥ 0
Definition 4.6 (Relative derivative). Let z ∗ ∈ R d × R ≥ 0 d z_\ast\in\mathbb{R}^d \times\mathbb{R}_{\geq0}^d z ∗ ∈ R d × R ≥ 0 d F : R d × R ≥ 0 d → R 2 d F:\mathbb{R}^d \times\mathbb{R}_{\geq0}^d \to\mathbb{R}^{2d} F : R d × R ≥ 0 d → R 2 d z ∗ z_\ast z ∗ R d × R ≥ 0 d \mathbb{R}^d \times\mathbb{R}_{\geq0}^d R d × R ≥ 0 d J : R 2 d → R 2 d J:\mathbb{R}^{2d} \to\mathbb{R}^{2d} J : R 2 d → R 2 d
lim z → z ∗ , z ≠ z ∗ z ∈ R d × R ≥ 0 d ∥ F ( z ) − F ( z ∗ ) − J ( z − z ∗ ) ∥ 2 ∥ z − z ∗ ∥ 2 = 0. \lim_{\substack{z\to z_\ast,\ z\ne z_\ast\\ z\in\mathbb{R}^d\times\mathbb{R}_{\geq0}^d}}
\frac{\lVert F(z)-F(z_\ast)-J(z-z_\ast)\rVert_2}{\lVert z-z_\ast\rVert_2}=0. z → z ∗ , z = z ∗ z ∈ R d × R ≥ 0 d lim ∥ z − z ∗ ∥ 2 ∥ F ( z ) − F ( z ∗ ) − J ( z − z ∗ ) ∥ 2 = 0. We call J J J relative derivative and write D R d × R ≥ 0 d F ( z ∗ ) = J D_{\mathbb{R}^d\times\mathbb{R}_{\geq0}^d}F(z_\ast)=J D R d × R ≥ 0 d F ( z ∗ ) = J
Lemma 4.7 (Relative derivative at the origin). For every h > 0 h>0 h > 0 F F F R d × R ≥ 0 d \mathbb{R}^d\times\mathbb{R}_{\geq0}^d R d × R ≥ 0 d
D R d × R ≥ 0 d F ( 0 ) = J 0 = [ I d − h ϵ A 0 0 β I d ] . (6) D_{\mathbb{R}^d\times\mathbb{R}_{\geq0}^d}F(0)=J_0=\begin{bmatrix} I_d-\frac{h}{\epsilon}A & 0 \\ 0 & \beta I_d \end{bmatrix}. \tag*{(6)} D R d × R ≥ 0 d F ( 0 ) = J 0 = [ I d − ϵ h A 0 0 β I d ] . (6) More precisely, F ( z ) = J 0 z + R ( z ) F(z)=J_0z+R(z) F ( z ) = J 0 z + R ( z ) ∥ R x ( z ) ∥ 2 = O ( ∥ z ∥ 2 3 / 2 ) \lVert R_x(z)\rVert_2=O(\lVert z\rVert_2^{3/2}) ∥ R x ( z ) ∥ 2 = O (∥ z ∥ 2 3/2 ) ∥ R v ( z ) ∥ 2 = O ( ∥ z ∥ 2 2 ) \lVert R_v(z)\rVert_2=O(\lVert z\rVert_2^2) ∥ R v ( z ) ∥ 2 = O (∥ z ∥ 2 2 ) z → 0 z\to0 z → 0 R d × R ≥ 0 d \mathbb{R}^d\times\mathbb{R}_{\geq0}^d R d × R ≥ 0 d
Proof. Write z = ( x , v ) ∈ R d × R ≥ 0 d z=(x,v)\in\mathbb{R}^d\times\mathbb{R}_{\geq0}^d z = ( x , v ) ∈ R d × R ≥ 0 d F 2 = F V F_2=F_V F 2 = F V ( h / ϵ ) A x (h/\epsilon)Ax ( h / ϵ ) A x F ( z ) = J 0 z + R ( z ) F(z)=J_0z+R(z) F ( z ) = J 0 z + R ( z ) J 0 J_0 J 0
R ( z ) = [ h A x ⊙ ( 1 ϵ 1 − 1 F 2 ( x , v ) + ϵ 1 ) ( 1 − β ) ( A x ) ⊙ 2 ] . R(z)=\begin{bmatrix}
hAx\odot\left(\frac{1}{\epsilon}\mathbf{1}-\frac{1}{\sqrt{F_2(x,v)}+\epsilon\mathbf{1}}\right)\\
(1-\beta)(Ax)^{\odot2}
\end{bmatrix}. R ( z ) = h A x ⊙ ( ϵ 1 1 − F 2 ( x , v ) + ϵ 1 1 ) ( 1 − β ) ( A x ) ⊙ 2 . Here and below the fractions in the first component are componentwise. Thus the two blocks of the remainder are R x ( z ) = h A x ⊙ ( F 2 ( x , v ) ⊘ ( ϵ ( F 2 ( x , v ) + ϵ 1 ) ) ) R_x(z)=hAx\odot\left(\sqrt{F_2(x,v)}\oslash\left(\epsilon\left(\sqrt{F_2(x,v)}+\epsilon\mathbf{1}\right)\right)\right) R x ( z ) = h A x ⊙ ( F 2 ( x , v ) ⊘ ( ϵ ( F 2 ( x , v ) + ϵ 1 ) ) ) R v ( z ) = ( 1 − β ) ( A x ) ⊙ 2 R_v(z)=(1-\beta)(Ax)^{\odot2} R v ( z ) = ( 1 − β ) ( A x ) ⊙ 2
To make the remainder estimate explicit, use the Euclidean norm and L A = ∥ A ∥ 2 L_A=\lVert A\rVert_2 L A = ∥ A ∥ 2 [ F 2 ( x , v ) ] i + ϵ ≥ ϵ \sqrt{[F_2(x,v)]_i}+\epsilon\geq\epsilon [ F 2 ( x , v ) ] i + ϵ ≥ ϵ ∣ ( A x ) i ∣ ≤ L A ∥ x ∥ 2 |(Ax)_i|\leq L_A\lVert x\rVert_2 ∣ ( A x ) i ∣ ≤ L A ∥ x ∥ 2 v i ≤ ∥ v ∥ 2 v_i\leq\lVert v\rVert_2 v i ≤ ∥ v ∥ 2 ( A x ) i 2 ≤ L A 2 ∥ x ∥ 2 2 (Ax)_i^2\leq L_A^2\lVert x\rVert_2^2 ( A x ) i 2 ≤ L A 2 ∥ x ∥ 2 2 ∣ ( R x ( z ) ) i ∣ ≤ ( h / ϵ 2 ) ∣ ( A x ) i ∣ [ F 2 ( x , v ) ] i |(R_x(z))_i|\leq(h/\epsilon^2)|(Ax)_i|\sqrt{[F_2(x,v)]_i} ∣ ( R x ( z ) ) i ∣ ≤ ( h / ϵ 2 ) ∣ ( A x ) i ∣ [ F 2 ( x , v ) ] i
∣ ( R x ( z ) ) i ∣ ≤ h L A ϵ 2 ∥ x ∥ 2 β ∥ v ∥ 2 + ( 1 − β ) L A 2 ∥ x ∥ 2 2 , ∥ R v ( z ) ∥ 2 ≤ ( 1 − β ) ∥ A x ∥ 2 2 ≤ ( 1 − β ) L A 2 ∥ x ∥ 2 2 . (7) \lvert(R_x(z))_i\rvert\leq\frac{hL_A}{\epsilon^2}\lVert x\rVert_2\sqrt{\beta\lVert v\rVert_2+(1-\beta)L_A^2\lVert x\rVert_2^2},\qquad\lVert R_v(z)\rVert_2\leq(1-\beta)\lVert Ax\rVert_2^2\leq(1-\beta)L_A^2\lVert x\rVert_2^2. \tag*{(7)} ∣( R x ( z ) ) i ∣ ≤ ϵ 2 h L A ∥ x ∥ 2 β ∥ v ∥ 2 + ( 1 − β ) L A 2 ∥ x ∥ 2 2 , ∥ R v ( z ) ∥ 2 ≤ ( 1 − β ) ∥ A x ∥ 2 2 ≤ ( 1 − β ) L A 2 ∥ x ∥ 2 2 . (7) If ϱ = ∥ z ∥ 2 ≤ 1 \varrho=\lVert z\rVert_2\leq1 ϱ = ∥ z ∥ 2 ≤ 1 ∥ x ∥ 2 , ∥ v ∥ 2 ≤ ϱ \lVert x\rVert_2,\lVert v\rVert_2\leq\varrho ∥ x ∥ 2 , ∥ v ∥ 2 ≤ ϱ ϱ 2 ≤ ϱ \varrho^2\leq\varrho ϱ 2 ≤ ϱ
∥ R x ( z ) ∥ 2 ≤ h d L A ϵ 2 β + ( 1 − β ) L A 2 ϱ 3 / 2 , ∥ R v ( z ) ∥ 2 ≤ ( 1 − β ) L A 2 ϱ 2 . \lVert R_x(z)\rVert_2\leq\frac{h\sqrt{d}L_A}{\epsilon^2}\sqrt{\beta+(1-\beta)L_A^2}\,\varrho^{3/2},\qquad\lVert R_v(z)\rVert_2\leq(1-\beta)L_A^2\varrho^2. ∥ R x ( z ) ∥ 2 ≤ ϵ 2 h d L A β + ( 1 − β ) L A 2 ϱ 3/2 , ∥ R v ( z ) ∥ 2 ≤ ( 1 − β ) L A 2 ϱ 2 . Consequently,
∥ R ( z ) ∥ 2 ∥ z ∥ 2 ≤ h d L A ϵ 2 β + ( 1 − β ) L A 2 ϱ 1 / 2 + ( 1 − β ) L A 2 ϱ ⟶ 0 , \frac{\lVert R(z)\rVert_2}{\lVert z\rVert_2}\leq\frac{h\sqrt{d}L_A}{\epsilon^2}\sqrt{\beta+(1-\beta)L_A^2}\,\varrho^{1/2}+(1-\beta)L_A^2\varrho\longrightarrow0, ∥ z ∥ 2 ∥ R ( z ) ∥ 2 ≤ ϵ 2 h d L A β + ( 1 − β ) L A 2 ϱ 1/2 + ( 1 − β ) L A 2 ϱ ⟶ 0 , as z → 0 z\to0 z → 0 R d × R ≥ 0 d \mathbb{R}^d\times\mathbb{R}_{\geq0}^d R d × R ≥ 0 d R ( z ) = o ( ∥ z ∥ 2 ) R(z)=o(\lVert z\rVert_2) R ( z ) = o (∥ z ∥ 2 )
Note that the classic Hartman–Grobman/stable manifold theory cannot be applied off the shelf, as the map is not differentiable in a neighborhood of the origin. This said, we find that a Lipschitz extension of the RMSProp map permits application of the Lipschitz regularity unstable manifold construction of [22 ] (see also [23 ] ), whose proof is based on the Lipschitz Hadamard–Perron theorem to obtain the instability of the origin, and even the existence of a low regularity unstable manifold.
We do not present this argument here, giving instead a direct instability proof based on the construction of a simple invariant cone. The argument is fairly short, makes the present paper more self-contained, and provides some insight into the geometry of the instability which will be useful later. The construction and preparation of the Lipschitz extension, on the other hand, are somewhat technical, and more opaque, and not genuinely shorter than the explicit construction given below.
Theorem 4.8 (Local Lyapunov instability at the origin). If
h > h ∗ : = 2 ϵ L A , h>h_*:=\frac{2\epsilon}{L_A}, h > h ∗ := L A 2 ϵ , then the origin is Lyapunov unstable under F F F
Proof. By Lemma 4.7, F ( z ) = J 0 z + R ( z ) F(z)=J_0z+R(z) F ( z ) = J 0 z + R ( z ) R ( z ) = o ( ∥ z ∥ 2 ) R(z)=o(\lVert z\rVert_2) R ( z ) = o (∥ z ∥ 2 ) z → 0 z\to0 z → 0 z ∈ R d × R ≥ 0 d z\in\mathbb{R}^d\times\mathbb{R}_{\geq0}^d z ∈ R d × R ≥ 0 d U : = E L A × { 0 } U:=E_{L_A}\times\{0\} U := E L A × { 0 } E L A E_{L_A} E L A A A A L A L_A L A U ‾ \overline{U} U J 0 ∣ U = − a I J_0|_U=-aI J 0 ∣ U = − a I a : = h L A ϵ − 1 > 1 a:=\frac{hL_A}{\epsilon}-1>1 a := ϵ h L A − 1 > 1 A A A U ⊥ U^\perp U ⊥ J 0 J_0 J 0 b : = ∥ J 0 ∣ U ⊥ ∥ 2 < a b := \lVert J_0|_{U^\perp}\rVert_2 < a b := ∥ J 0 ∣ U ⊥ ∥ 2 < a x x x J 0 J_0 J 0 1 − h ϵ λ j 1-\frac{h}{\epsilon}\lambda_j 1 − ϵ h λ j 0 ≤ λ j ≤ L A 0 \le\lambda_j \le L_A 0 ≤ λ j ≤ L A v v v β \beta β h L A ϵ > 2 \frac{hL_A}{\epsilon} > 2 ϵ h L A > 2 λ j < L A \lambda_j < L_A λ j < L A β \beta β a a a
Choose η > 0 \eta> 0 η > 0
a − 2 η > max { 1 , b + 2 η } . (8) a - 2\eta> \max\{1,b+2\eta\}. \tag*{(8)} a − 2 η > max { 1 , b + 2 η } . (8) Since R ( z ) = o ( ∥ z ∥ 2 ) R(z)=o(\lVert z\rVert_2) R ( z ) = o (∥ z ∥ 2 ) ϱ 0 > 0 \varrho_0 > 0 ϱ 0 > 0 ∥ R ( z ) ∥ 2 ≤ η ∥ z ∥ 2 \lVert R(z)\rVert_2 \le\eta\lVert z\rVert_2 ∥ R ( z ) ∥ 2 ≤ η ∥ z ∥ 2 z ∈ B ϱ 0 ( 0 ) ∩ ( R d × R ≥ 0 d ) z \in B_{\varrho_0}(0) \cap(\mathbb{R}^d \times\mathbb{R}_{\ge0}^d) z ∈ B ϱ 0 ( 0 ) ∩ ( R d × R ≥ 0 d )
Let Π U \Pi_U Π U Π U ⊥ \Pi_{U^\perp} Π U ⊥ U U U U ⊥ U^\perp U ⊥
C : = { z = u + w : u ∈ U , w ∈ U ⊥ , ∥ w ∥ 2 ≤ ∥ u ∥ 2 } . \mathcal{C}:=\{z=u+w:u\in U,w\in U^\perp,\ \lVert w\rVert_2\le\lVert u\rVert_2\}. C := { z = u + w : u ∈ U , w ∈ U ⊥ , ∥ w ∥ 2 ≤ ∥ u ∥ 2 } . If z = u + w ∈ C ∩ B ϱ 0 ( 0 ) ∩ ( R d × R ≥ 0 d ) z=u+w\in\mathcal{C}\cap B_{\varrho_0}(0)\cap(\mathbb{R}^d\times\mathbb{R}_{\ge0}^d) z = u + w ∈ C ∩ B ϱ 0 ( 0 ) ∩ ( R d × R ≥ 0 d ) ∥ z ∥ 2 ≤ ∥ u ∥ 2 + ∥ w ∥ 2 ≤ 2 ∥ u ∥ 2 \lVert z\rVert_2\le\lVert u\rVert_2+\lVert w\rVert_2\le2\lVert u\rVert_2 ∥ z ∥ 2 ≤ ∥ u ∥ 2 + ∥ w ∥ 2 ≤ 2 ∥ u ∥ 2
∥ Π U F ( z ) ∥ 2 ≥ a ∥ u ∥ 2 − ∥ R ( z ) ∥ 2 ≥ ( a − 2 η ) ∥ u ∥ 2 , \lVert\Pi_U F(z)\rVert_2\ge a\lVert u\rVert_2-\lVert R(z)\rVert_2\ge(a-2\eta)\lVert u\rVert_2, ∥ Π U F ( z ) ∥ 2 ≥ a ∥ u ∥ 2 − ∥ R ( z ) ∥ 2 ≥ ( a − 2 η ) ∥ u ∥ 2 , while
∥ Π U ⊥ F ( z ) ∥ 2 ≤ b ∥ w ∥ 2 + ∥ R ( z ) ∥ 2 ≤ ( b + 2 η ) ∥ u ∥ 2 . \lVert\Pi_{U^\perp}F(z)\rVert_2\le b\lVert w\rVert_2+\lVert R(z)\rVert_2\le(b+2\eta)\lVert u\rVert_2. ∥ Π U ⊥ F ( z ) ∥ 2 ≤ b ∥ w ∥ 2 + ∥ R ( z ) ∥ 2 ≤ ( b + 2 η ) ∥ u ∥ 2 . By (8), ∥ Π U ⊥ F ( z ) ∥ 2 < ∥ Π U F ( z ) ∥ 2 \lVert\Pi_{U^\perp}F(z)\rVert_2<\lVert\Pi_UF(z)\rVert_2 ∥ Π U ⊥ F ( z ) ∥ 2 < ∥ Π U F ( z ) ∥ 2 ∥ Π U F ( z ) ∥ 2 ≥ ( a − 2 η ) ∥ u ∥ 2 > ∥ u ∥ 2 \lVert\Pi_UF(z)\rVert_2\ge(a-2\eta)\lVert u\rVert_2>\lVert u\rVert_2 ∥ Π U F ( z ) ∥ 2 ≥ ( a − 2 η ) ∥ u ∥ 2 > ∥ u ∥ 2 B ϱ 0 ( 0 ) B_{\varrho_0}(0) B ϱ 0 ( 0 ) U U U q : = a − 2 η > 1 q:=a-2\eta>1 q := a − 2 η > 1
Let e max ∈ E L A e_{\max}\in E_{L_A} e m a x ∈ E L A δ > 0 \delta>0 δ > 0 0 < t < min { δ , ϱ 0 } 0<t<\min\{\delta,\varrho_0\} 0 < t < min { δ , ϱ 0 } z 0 = ( t e max , 0 ) z_0=(te_{\max},0) z 0 = ( t e m a x , 0 ) z 0 ∈ C ∩ ( R d × R ≥ 0 d ) z_0\in\mathcal{C}\cap(\mathbb{R}^d\times\mathbb{R}_{\ge0}^d) z 0 ∈ C ∩ ( R d × R ≥ 0 d ) B ϱ 0 ( 0 ) B_{\varrho_0}(0) B ϱ 0 ( 0 ) k ≥ 0 k\ge0 k ≥ 0 ∥ Π U F k ( z 0 ) ∥ 2 ≥ q k t \lVert\Pi_UF^k(z_0)\rVert_2\ge q^kt ∥ Π U F k ( z 0 ) ∥ 2 ≥ q k t q > 1 q>1 q > 1 B ϱ 0 ( 0 ) B_{\varrho_0}(0) B ϱ 0 ( 0 )
Global trapping
We now show that every RMSProp orbit eventually enters a compact forward-invariant set. The argument uses a Lyapunov or energy-like function that decreases outside a bounded region. Indeed, the argument is similar to that found in the proof of Theorem 4.4, but adapted to the quadratic loss function and taking into account the loss of stability of the origin. We first state a general trapping criterion, then establish the required decrease for normalized RMSProp and obtain bounds in the original coordinates.
Theorem 4.9. Let F : R d × R ≥ 0 d → R d × R ≥ 0 d F:\mathbb{R}^d\times\mathbb{R}_{\ge0}^d\to\mathbb{R}^d\times\mathbb{R}_{\ge0}^d F : R d × R ≥ 0 d → R d × R ≥ 0 d V : R d × R ≥ 0 d → R V:\mathbb{R}^d\times\mathbb{R}_{\ge0}^d\to\mathbb{R} V : R d × R ≥ 0 d → R V V V
V ( z ) ⟶ + ∞ as ∥ z ∥ ⟶ ∞ , z ∈ R d × R ≥ 0 d . (9) V(z)\longrightarrow+\infty\quad\text{as }\lVert z\rVert\longrightarrow\infty,\qquad z\in\mathbb{R}^d\times\mathbb{R}_{\ge0}^d. \tag*{(9)} V ( z ) ⟶ + ∞ as ∥ z ∥ ⟶ ∞ , z ∈ R d × R ≥ 0 d . (9) Let R ⊂ R d × R ≥ 0 d R\subset\mathbb{R}^d\times\mathbb{R}_{\ge0}^d R ⊂ R d × R ≥ 0 d
V ( F ( z ) ) < V ( z ) for every z ∈ ( R d × R ≥ 0 d ) ∖ R . (10) V(F(z))<V(z)\qquad\text{for every }z\in(\mathbb{R}^d\times\mathbb{R}_{\ge0}^d)\setminus R. \tag*{(10)} V ( F ( z )) < V ( z ) for every z ∈ ( R d × R ≥ 0 d ) ∖ R . (10) Then there exists C 0 > 0 C_0>0 C 0 > 0 C ≥ C 0 C\ge C_0 C ≥ C 0 Ω C = { z ∈ R d × R ≥ 0 d : V ( z ) ≤ C } \Omega_C=\{z\in\mathbb{R}^d\times\mathbb{R}_{\ge0}^d:V(z)\le C\} Ω C = { z ∈ R d × R ≥ 0 d : V ( z ) ≤ C }
(i) Compactness. Ω C \Omega_C Ω C
(ii) Containment of the region. R ⊂ int R d × R ≥ 0 d ( Ω C ) R\subset\operatorname{int}_{\mathbb{R}^d\times\mathbb{R}_{\ge0}^d}(\Omega_C) R ⊂ int R d × R ≥ 0 d ( Ω C )
(iii) Strict trapping. F ( Ω C ) ⊂ int R d × R ≥ 0 d ( Ω C ) F(\Omega_C)\subset\operatorname{int}_{\mathbb{R}^d\times\mathbb{R}_{\ge0}^d}(\Omega_C) F ( Ω C ) ⊂ int R d × R ≥ 0 d ( Ω C ) Ω C \Omega_C Ω C
(iv) Finite-time entry. For every z 0 ∈ R d × R ≥ 0 d z_0\in\mathbb{R}^d\times\mathbb{R}_{\ge0}^d z 0 ∈ R d × R ≥ 0 d N = N ( z 0 , C ) ≥ 0 N=N(z_0,C)\ge0 N = N ( z 0 , C ) ≥ 0 F n ( z 0 ) ∈ Ω C F^n(z_0)\in\Omega_C F n ( z 0 ) ∈ Ω C n ≥ N n\ge N n ≥ N
The interiors in parts (ii) and (iii) are relative to R d × R ≥ 0 d \mathbb{R}^d\times\mathbb{R}_{\ge0}^d R d × R ≥ 0 d
The theorem is closely related to standard results on Lyapunov functions for discrete time dynamical systems. The paper of Kalman and Bertram [26 ] as well as the work of LaSalle [33 ] are excellent resources. A proof of the theorem allowing application on the cone R d × R ≥ 0 d \mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d} R d × R ≥ 0 d
The rest of this section is devoted to showing that the epsilon-outside RMSProp map admits an energy satisfying conditions (9) and (10), which imply that every orbit eventually enters and remains in a compact trapping set.
Normalization and energy estimates
First, we introduce the notation and normalization used below. Let A = A T ≻ 0 A = A^{\mathsf{T}} \succ0 A = A T ≻ 0 d × d d \times d d × d 0 < β < 1 0 < \beta< 1 0 < β < 1 ϵ > 0 \epsilon> 0 ϵ > 0 m = λ min ( A ) m = \lambda_{\min}(A) m = λ m i n ( A ) L A = λ max ( A ) L_A = \lambda_{\max}(A) L A = λ m a x ( A ) h ∗ = 2 ϵ / L A h_* = 2\epsilon/L_A h ∗ = 2 ϵ / L A h > 0 h > 0 h > 0
X = x / h , W = v ⊙ / h , and c = ϵ / h . (11) X = x/h,\qquad W = \sqrt[\odot]{v}/h,\qquad\text{and}\qquad c = \epsilon/h. \tag*{(11)} X = x / h , W = ⊙ v / h , and c = ϵ / h . (11) The normalized map G = ( G 1 , G 2 ) G = (G_1,G_2) G = ( G 1 , G 2 )
G 2 ( X , W ) = β W ⊙ 2 + ( 1 − β ) ( A X ) ⊙ 2 ⊙ , G_2(X,W) = \sqrt[\odot]{\beta W^{\odot2} + (1-\beta)(AX)^{\odot2}}, G 2 ( X , W ) = ⊙ β W ⊙ 2 + ( 1 − β ) ( A X ) ⊙ 2 , G 1 ( X , W ) = X − A X ⊙ ( G 2 ( X , W ) + c 1 ) . (12) G_1(X,W) = X - AX \odot\bigl(G_2(X,W) + c\mathbf{1}\bigr). \tag*{(12)} G 1 ( X , W ) = X − A X ⊙ ( G 2 ( X , W ) + c 1 ) . (12) For each fixed h > 0 h > 0 h > 0 R d × R ≥ 0 d \mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d} R d × R ≥ 0 d c > 0 c > 0 c > 0
For the following energy estimates, assume h ≥ h ∗ h \geq h_* h ≥ h ∗
V ( X , W ) = 1 2 X T A X + α ∥ W ∥ 1 , (13) \mathcal{V}(X,W) = \frac{1}{2}X^{\mathsf{T}}AX + \alpha\lVert W\rVert_1, \tag*{(13)} V ( X , W ) = 2 1 X T A X + α ∥ W ∥ 1 , (13) where α > 0 \alpha> 0 α > 0
V ( X , W ) ≥ m 2 ∥ X ∥ 2 2 + α ∥ W ∥ 1 ⟶ + ∞ as ∥ ( X , W ) ∥ ⟶ ∞ , \mathcal{V}(X,W) \geq\frac{m}{2}\lVert X\rVert_2^2 + \alpha\lVert W\rVert_1 \longrightarrow+\infty\quad\text{as}\quad\lVert(X,W)\rVert\longrightarrow\infty, V ( X , W ) ≥ 2 m ∥ X ∥ 2 2 + α ∥ W ∥ 1 ⟶ + ∞ as ∥( X , W )∥ ⟶ ∞ , with ( X , W ) ∈ R d × R ≥ 0 d (X,W) \in\mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d} ( X , W ) ∈ R d × R ≥ 0 d
To establish condition (10), we first bound the changes in the spatial and memory terms separately. We then estimate their combined contribution coordinate by coordinate, distinguishing whether the memory or the gradient is large. This shows that the energy decreases outside a bounded set. The trapping criterion then applies, and an explicit sublevel bound gives the corresponding size estimates in the original coordinates.
Remark 4.10 (Choice of energy for the scaling bounds). We choose the energy (13) with the bounds on the size of the trapping region in mind. After the normalization (11), the learning rate appears in the map only through c = ϵ / h c = \epsilon/h c = ϵ / h h ≥ h ∗ h \geq h_* h ≥ h ∗ h h h h h h x x x h 2 h^2 h 2 v v v [14 ] use a Lyapunov function to prove global convergence and boundedness under bounded time-varying step sizes. Their function includes a 3 / 2 3/2 3/2
Lemma 4.11 (Normalized Lyapunov increments). Under the preceding hypotheses and notation, set u = X − G 1 ( X , W ) u = X - G_1(X,W) u = X − G 1 ( X , W ) c ∗ = L A / 2 c_* = L_A/2 c ∗ = L A /2 c ≤ c ∗ c \leq c_* c ≤ c ∗
K = max { 1 , 4 1 − β 1 − β } , M = K β + 1 − β , α = 1 8 M 1 − β , K = \max\left\{1,\frac{4\sqrt{1-\beta}}{1-\sqrt{\beta}}\right\},\qquad M = K\sqrt{\beta}+\sqrt{1-\beta},\qquad\alpha= \frac{1}{8M\sqrt{1-\beta}}, K = max { 1 , 1 − β 4 1 − β } , M = K β + 1 − β , α = 8 M 1 − β 1 , and use V \mathcal{V} V
V ( G ( X , W ) ) − V ( X , W ) ≤ ∑ j = 1 d B j + L A d 2 ( 1 − β ) , \mathcal{V}(G(X,W))-\mathcal{V}(X,W) \leq\sum_{j=1}^{d} B_j + \frac{L_A d}{2(1-\beta)}, V ( G ( X , W )) − V ( X , W ) ≤ j = 1 ∑ d B j + 2 ( 1 − β ) L A d , where
B j = − ∣ ( A X ) j ∣ 2 β W j + 1 − β ∣ ( A X ) j ∣ + c ∗ + α 1 − β ∣ ( A X ) j ∣ − α ( 1 − β ) W j . (14) B_j=-\frac{|(AX)_j|^2}{\sqrt{\beta}W_j+\sqrt{1-\beta}|(AX)_j|+c_*}+\alpha\sqrt{1-\beta}|(AX)_j|-\alpha(1-\sqrt{\beta})W_j. \tag*{(14)} B j = − β W j + 1 − β ∣ ( A X ) j ∣ + c ∗ ∣ ( A X ) j ∣ 2 + α 1 − β ∣ ( A X ) j ∣ − α ( 1 − β ) W j . (14) Moreover, ∥ u ∥ 2 2 ≤ d / ( 1 − β ) \lVert u\rVert_2^2\le d/(1-\beta) ∥ u ∥ 2 2 ≤ d / ( 1 − β )
Proof. Fix ( X , W ) ∈ R d × R ≥ 0 d (X,W)\in\mathbb{R}^d\times\mathbb{R}_{\geq0}^d ( X , W ) ∈ R d × R ≥ 0 d
V ( G ( X , W ) ) − V ( X , W ) = 1 2 ( G 1 ( X , W ) T A G 1 ( X , W ) − X T A X ) + α ( ∥ G 2 ( X , W ) ∥ 1 − ∥ W ∥ 1 ) . \mathcal{V}(G(X,W))-\mathcal{V}(X,W)=\frac{1}{2}\left(G_1(X,W)^{\mathsf{T}}AG_1(X,W)-X^{\mathsf{T}}AX\right)+\alpha\left(\lVert G_2(X,W)\rVert_1-\lVert W\rVert_1\right). V ( G ( X , W )) − V ( X , W ) = 2 1 ( G 1 ( X , W ) T A G 1 ( X , W ) − X T A X ) + α ( ∥ G 2 ( X , W ) ∥ 1 − ∥ W ∥ 1 ) . We bound the first term and the second term individually, then add their upper bounds. To obtain the estimates needed for both terms, the memory update gives
[ G 2 ( X , W ) ] j 2 = β W j 2 + ( 1 − β ) ( A X ) j 2 ≥ ( 1 − β ) ( A X ) j 2 . [G_2(X,W)]_j^2=\beta W_j^2+(1-\beta)(AX)_j^2\geq(1-\beta)(AX)_j^2. [ G 2 ( X , W ) ] j 2 = β W j 2 + ( 1 − β ) ( A X ) j 2 ≥ ( 1 − β ) ( A X ) j 2 . Taking nonnegative square roots yields the lower bound below. For the upper bound, apply a + b ≤ a + b \sqrt{a+b}\leq\sqrt{a}+\sqrt{b} a + b ≤ a + b a , b ≥ 0 a,b\geq0 a , b ≥ 0 W j ≥ 0 W_j\geq0 W j ≥ 0
1 − β ∣ ( A X ) j ∣ ≤ [ G 2 ( X , W ) ] j ≤ β W j + 1 − β ∣ ( A X ) j ∣ . (15) \sqrt{1-\beta}|(AX)_j|\leq[G_2(X,W)]_j\leq\sqrt{\beta}W_j+\sqrt{1-\beta}|(AX)_j|. \tag*{(15)} 1 − β ∣ ( A X ) j ∣ ≤ [ G 2 ( X , W ) ] j ≤ β W j + 1 − β ∣ ( A X ) j ∣. (15) Spatial term. Substituting G 1 ( X , W ) = X − u G_1(X,W)=X-u G 1 ( X , W ) = X − u A = A T A=A^{\mathsf{T}} A = A T
1 2 ( G 1 ( X , W ) T A G 1 ( X , W ) − X T A X ) = − u T A X + 1 2 u T A u = − ∑ j ( A X ) j 2 [ G 2 ( X , W ) ] j + c + 1 2 u T A u . \frac{1}{2}\left(G_1(X,W)^{\mathsf{T}}AG_1(X,W)-X^{\mathsf{T}}AX\right)=-u^{\mathsf{T}}AX+\frac{1}{2}u^{\mathsf{T}}Au=-\sum_j\frac{(AX)_j^2}{[G_2(X,W)]_j+c}+\frac{1}{2}u^{\mathsf{T}}Au. 2 1 ( G 1 ( X , W ) T A G 1 ( X , W ) − X T A X ) = − u T A X + 2 1 u T A u = − j ∑ [ G 2 ( X , W ) ] j + c ( A X ) j 2 + 2 1 u T A u . The lower memory bound implies
∣ u j ∣ = ∣ ( A X ) j ∣ [ G 2 ( X , W ) ] j + c ≤ ∣ ( A X ) j ∣ 1 − β ∣ ( A X ) j ∣ + c ≤ 1 1 − β . |u_j|=\frac{|(AX)_j|}{[G_2(X,W)]_j+c}\leq\frac{|(AX)_j|}{\sqrt{1-\beta}|(AX)_j|+c}\leq\frac{1}{\sqrt{1-\beta}}. ∣ u j ∣ = [ G 2 ( X , W ) ] j + c ∣ ( A X ) j ∣ ≤ 1 − β ∣ ( A X ) j ∣ + c ∣ ( A X ) j ∣ ≤ 1 − β 1 . This also holds when ( A X ) j = 0 (AX)_j=0 ( A X ) j = 0 c > 0 c>0 c > 0 ∥ u ∥ 2 2 ≤ d / ( 1 − β ) \lVert u\rVert_2^2\leq d/(1-\beta) ∥ u ∥ 2 2 ≤ d / ( 1 − β ) u T A u ≤ L A ∥ u ∥ 2 2 u^{\mathsf{T}}Au\leq L_A\lVert u\rVert_2^2 u T A u ≤ L A ∥ u ∥ 2 2 u T A u / 2 ≤ L A d / [ 2 ( 1 − β ) ] u^{\mathsf{T}}Au/2\leq L_A d/[2(1-\beta)] u T A u /2 ≤ L A d / [ 2 ( 1 − β )] c ≤ c ∗ c\leq c_* c ≤ c ∗ [ G 2 ( X , W ) ] j + c ≤ β W j + 1 − β ∣ ( A X ) j ∣ + c ∗ [G_2(X,W)]_j+c\leq\sqrt{\beta}W_j+\sqrt{1-\beta}|(AX)_j|+c_* [ G 2 ( X , W ) ] j + c ≤ β W j + 1 − β ∣ ( A X ) j ∣ + c ∗ t ↦ − g 2 / t t\mapsto-g^2/t t ↦ − g 2 / t t > 0 t>0 t > 0
1 2 G 1 ( X , W ) T A G 1 ( X , W ) − 1 2 X T A X ≤ − ∑ j ∣ ( A X ) j ∣ 2 β W j + 1 − β ∣ ( A X ) j ∣ + c ∗ + L A d 2 ( 1 − β ) . (16) \frac{1}{2}G_1(X,W)^{\mathsf{T}}AG_1(X,W)-\frac{1}{2}X^{\mathsf{T}}AX\leq-\sum_j\frac{|(AX)_j|^2}{\sqrt{\beta}W_j+\sqrt{1-\beta}|(AX)_j|+c_*}+\frac{L_A d}{2(1-\beta)}. \tag*{(16)} 2 1 G 1 ( X , W ) T A G 1 ( X , W ) − 2 1 X T A X ≤ − j ∑ β W j + 1 − β ∣ ( A X ) j ∣ + c ∗ ∣ ( A X ) j ∣ 2 + 2 ( 1 − β ) L A d . (16) Memory term. Nonnegativity of W W W G 2 ( X , W ) G_2(X,W) G 2 ( X , W ) 1 1 1 W j W_j W j α \alpha α
α ( ∥ G 2 ( X , W ) ∥ 1 − ∥ W ∥ 1 ) = α ∑ j ( [ G 2 ( X , W ) ] j − W j ) ≤ α ∑ j [ 1 − β ∣ ( A X ) j ∣ − ( 1 − β ) W j ] . (17) \alpha\left(\lVert G_2(X,W)\rVert_1-\lVert W\rVert_1\right)=\alpha\sum_j\left([G_2(X,W)]_j-W_j\right)\leq\alpha\sum_j\left[\sqrt{1-\beta}|(AX)_j|-(1-\sqrt{\beta})W_j\right]. \tag*{(17)} α ( ∥ G 2 ( X , W ) ∥ 1 − ∥ W ∥ 1 ) = α j ∑ ( [ G 2 ( X , W ) ] j − W j ) ≤ α j ∑ [ 1 − β ∣ ( A X ) j ∣ − ( 1 − β ) W j ] . (17) Sum of the two terms. Adding (16) and (17) gives
V ( G ( X , W ) ) − V ( X , W ) ≤ ∑ j = 1 d B j + L A d 2 ( 1 − β ) , \mathcal{V}(G(X,W))-\mathcal{V}(X,W)\leq\sum_{j=1}^{d}B_j+\frac{L_A d}{2(1-\beta)}, V ( G ( X , W )) − V ( X , W ) ≤ j = 1 ∑ d B j + 2 ( 1 − β ) L A d , where B j B_j B j
Lemma 4.12 (Uniform bound on the energy change). Under the preceding hypotheses and notation, define
γ = min { α ( 1 − β ) 8 , 1 8 M ( 1 + K ) } , H 0 = c ∗ M ( α 1 − β + γ ( 1 + K ) ) , H = L A d 2 ( 1 − β ) + d H 0 . \gamma=\min\left\{\frac{\alpha(1-\sqrt{\beta})}{8},\frac{1}{8M(1+K)}\right\},\qquad H_0=\frac{c_*}{M}\left(\alpha\sqrt{1-\beta}+\gamma(1+K)\right),\qquad H=\frac{L_A d}{2(1-\beta)}+dH_0. γ = min { 8 α ( 1 − β ) , 8 M ( 1 + K ) 1 } , H 0 = M c ∗ ( α 1 − β + γ ( 1 + K ) ) , H = 2 ( 1 − β ) L A d + d H 0 . These constants are positive and depend only on ( A , β , d ) (A,\beta,d) ( A , β , d ) ( X , W ) (X,W) ( X , W ) j j j
B j ≤ − γ ( ∣ ( A X ) j ∣ + W j ) + H 0 , (18) B_j\leq-\gamma\left(|(AX)_j|+W_j\right)+H_0, \tag*{(18)} B j ≤ − γ ( ∣ ( A X ) j ∣ + W j ) + H 0 , (18) V ( G ( X , W ) ) − V ( X , W ) ≤ − γ ( ∥ A X ∥ 1 + ∥ W ∥ 1 ) + H . (19) \mathcal{V}(G(X,W))-\mathcal{V}(X,W)\leq-\gamma\left(\lVert AX\rVert_1+\lVert W\rVert_1\right)+H. \tag*{(19)} V ( G ( X , W )) − V ( X , W ) ≤ − γ ( ∥ A X ∥ 1 + ∥ W ∥ 1 ) + H . (19) Proof. Fix ( X , W ) (X,W) ( X , W ) j j j g = ∣ ( A X ) j ∣ g = |(AX)_j| g = ∣ ( A X ) j ∣ w = W j w = W_j w = W j B = B j B = B_j B = B j
Case 1: w ≥ K g w \ge Kg w ≥ K g Dropping the nonpositive fractional term and using K ≥ 4 1 − β / ( 1 − β ) K \ge4\sqrt{1-\beta}/(1-\sqrt{\beta}) K ≥ 4 1 − β / ( 1 − β )
B ≤ α 1 − β g − α ( 1 − β ) w ≤ − 3 4 α ( 1 − β ) w ≤ − γ ( g + w ) , B \le\alpha\sqrt{1-\beta}g-\alpha(1-\sqrt{\beta})w \le-\frac{3}{4}\alpha(1-\sqrt{\beta})w \le-\gamma(g+w), B ≤ α 1 − β g − α ( 1 − β ) w ≤ − 4 3 α ( 1 − β ) w ≤ − γ ( g + w ) , since K ≥ 1 K \ge1 K ≥ 1 g + w ≤ 2 w g+w \le2w g + w ≤ 2 w γ ≤ α ( 1 − β ) / 8 \gamma\le\alpha(1-\sqrt{\beta})/8 γ ≤ α ( 1 − β ) /8
Case 2: w < K g w < Kg w < K g g ≥ c ∗ / M g \ge c_*/M g ≥ c ∗ / M The denominator in B B B M g + c ∗ ≤ 2 M g Mg+c_* \le2Mg M g + c ∗ ≤ 2 M g α 1 − β = 1 / ( 8 M ) \alpha\sqrt{1-\beta}=1/(8M) α 1 − β = 1/ ( 8 M )
B ≤ − g 2 M + α 1 − β g = − 3 g 8 M ≤ − γ ( g + w ) , B \le-\frac{g}{2M}+\alpha\sqrt{1-\beta}g=-\frac{3g}{8M}\le-\gamma(g+w), B ≤ − 2 M g + α 1 − β g = − 8 M 3 g ≤ − γ ( g + w ) , because g + w < ( 1 + K ) g g+w < (1+K)g g + w < ( 1 + K ) g γ ≤ 1 / [ 8 M ( 1 + K ) ] \gamma\le1/[8M(1+K)] γ ≤ 1/ [ 8 M ( 1 + K )]
Case 3: w < K g w < Kg w < K g g < c ∗ / M g < c_*/M g < c ∗ / M Dropping both nonpositive terms gives
B + γ ( g + w ) ≤ α 1 − β g + γ ( g + w ) ≤ c ∗ M ( α 1 − β + γ ( 1 + K ) ) = H 0 . B+\gamma(g+w)\le\alpha\sqrt{1-\beta}g+\gamma(g+w)\le\frac{c_*}{M}\bigl(\alpha\sqrt{1-\beta}+\gamma(1+K)\bigr)=H_0. B + γ ( g + w ) ≤ α 1 − β g + γ ( g + w ) ≤ M c ∗ ( α 1 − β + γ ( 1 + K ) ) = H 0 . Since H 0 > 0 H_0>0 H 0 > 0 (18) . Summing over j j j Lemma 4.11 gives (19) . □
Define the bounded set
R ∗ = 2 H γ , B = { ( X , W ) ∈ R d × R ≥ 0 d : ∥ A X ∥ 1 + ∥ W ∥ 1 ≤ R ∗ } . (20) R_*=\frac{2H}{\gamma},\qquad\mathcal{B}=\{(X,W)\in\mathbb{R}^d\times\mathbb{R}_{\ge0}^d:\|AX\|_1+\|W\|_1\le R_*\}. \tag*{(20)} R ∗ = γ 2 H , B = {( X , W ) ∈ R d × R ≥ 0 d : ∥ A X ∥ 1 + ∥ W ∥ 1 ≤ R ∗ } . (20) On B \mathcal{B} B m ∥ X ∥ 2 ≤ ∥ A X ∥ 2 ≤ R ∗ m\|X\|_2\le\|AX\|_2\le R_* m ∥ X ∥ 2 ≤ ∥ A X ∥ 2 ≤ R ∗ ∥ W ∥ 1 ≤ R ∗ \|W\|_1\le R_* ∥ W ∥ 1 ≤ R ∗ B \mathcal{B} B (19) gives
V ( G ( X , W ) ) − V ( X , W ) < − γ R ∗ + H = − H < 0. \mathcal{V}(G(X,W))-\mathcal{V}(X,W)<-\gamma R_*+H=-H<0. V ( G ( X , W )) − V ( X , W ) < − γ R ∗ + H = − H < 0. Thus the hypotheses of Theorem 4.9 hold. We now obtain a trapping set of the form { V ≤ C } \{\mathcal{V}\le C\} { V ≤ C } C C C h h h
Theorem 4.13 (RMSProp trapping region and scaling bounds). Under the preceding hypotheses and notation, define
B 0 = L A R ∗ 2 2 m 2 + α R ∗ , C = B 0 + 2 H , C x = 2 C / m , C v = ( C / α ) 2 . B_0=\frac{L_A R_*^2}{2m^2}+\alpha R_*,\qquad C=B_0+2H,\qquad C_x=\sqrt{2C/m},\qquad C_v=(C/\alpha)^2. B 0 = 2 m 2 L A R ∗ 2 + α R ∗ , C = B 0 + 2 H , C x = 2 C / m , C v = ( C / α ) 2 . These positive constants depend only on ( A , β , d ) (A,\beta,d) ( A , β , d ) h ≥ h ∗ h\ge h_* h ≥ h ∗ S = { ( X , W ) ∈ R d × R ≥ 0 d : V ( X , W ) ≤ C } S=\{(X,W)\in\mathbb{R}^d\times\mathbb{R}_{\ge0}^d:\mathcal{V}(X,W)\le C\} S = {( X , W ) ∈ R d × R ≥ 0 d : V ( X , W ) ≤ C }
B ⊂ int R d × R ≥ 0 d ( S ) , G ( S ) ⊂ int R d × R ≥ 0 d ( S ) . \mathcal{B}\subset\operatorname{int}_{\mathbb{R}^d\times\mathbb{R}_{\ge0}^d}(S),\qquad G(S)\subset\operatorname{int}_{\mathbb{R}^d\times\mathbb{R}_{\ge0}^d}(S). B ⊂ int R d × R ≥ 0 d ( S ) , G ( S ) ⊂ int R d × R ≥ 0 d ( S ) . Every normalized orbit enters S S S v 0 ≥ 0 v_0\ge0 v 0 ≥ 0
K h = { ( x , v ) ∈ R d × R ≥ 0 d : V ( x / h , ⊙ v / h ) ≤ C } ⊂ [ − C x h , C x h ] d × [ 0 , C v h 2 ] d , \mathcal{K}_h=\{(x,v)\in\mathbb{R}^d\times\mathbb{R}_{\ge0}^d:\mathcal{V}(x/h,\sqrt{\mathord\odot v}/h)\le C\}
\subset[-C_xh,C_xh]^d\times[0,C_vh^2]^d, K h = {( x , v ) ∈ R d × R ≥ 0 d : V ( x / h , ⊙ v / h ) ≤ C } ⊂ [ − C x h , C x h ] d × [ 0 , C v h 2 ] d , with F ( K h ) ⊂ int R d × R ≥ 0 d ( K h ) F(\mathcal{K}_h)\subset\operatorname{int}_{\mathbb{R}^d\times\mathbb{R}_{\ge0}^d}(\mathcal{K}_h) F ( K h ) ⊂ int R d × R ≥ 0 d ( K h ) h ≥ h ∗ h\ge h_* h ≥ h ∗
lim sup n → ∞ ∥ x n ∥ ∞ ≤ C x h and lim sup n → ∞ ∥ v n ∥ ∞ ≤ C v h 2 . \limsup_{n\to\infty}\|x_n\|_\infty\le C_xh\qquad\text{and}\qquad\limsup_{n\to\infty}\|v_n\|_\infty\le C_vh^2. n → ∞ lim sup ∥ x n ∥ ∞ ≤ C x h and n → ∞ lim sup ∥ v n ∥ ∞ ≤ C v h 2 . Proof. On B \mathcal{B} B (20) imply V ≤ B 0 \mathcal{V}\le B_0 V ≤ B 0 (19) also gives V ∘ G ≤ V + H \mathcal{V}\circ G\le\mathcal{V}+H V ∘ G ≤ V + H R d × R ≥ 0 d \mathbb{R}^d\times\mathbb{R}_{\ge0}^d R d × R ≥ 0 d
V ( z ) ≤ B 0 < C , V ( G ( z ) ) ≤ B 0 + H < C ( z ∈ B ) . \mathcal{V}(z)\le B_0<C,\qquad\mathcal{V}(G(z))\le B_0+H<C\qquad(z\in\mathcal{B}). V ( z ) ≤ B 0 < C , V ( G ( z )) ≤ B 0 + H < C ( z ∈ B ) . These bounds show that the level C C C § 14 .
On S S S ( m / 2 ) ∥ X ∥ ∞ 2 ≤ V ≤ C (m/2)\|X\|_\infty^2\le\mathcal{V}\le C ( m /2 ) ∥ X ∥ ∞ 2 ≤ V ≤ C α ∥ W ∥ ∞ 2 ≤ V ≤ C \alpha\|W\|_\infty^2\le\mathcal{V}\le C α ∥ W ∥ ∞ 2 ≤ V ≤ C ∥ X ∥ ∞ ≤ C x \|X\|_\infty\le C_x ∥ X ∥ ∞ ≤ C x ∥ W ∥ ∞ 2 ≤ C v \|W\|_\infty^2\le C_v ∥ W ∥ ∞ 2 ≤ C v ( X , W ) ↦ ( h X , h 2 W ⊙ 2 ) (X,W)\mapsto(hX,h^2W^{\mathord\odot2}) ( X , W ) ↦ ( h X , h 2 W ⊙ 2 ) G G G F F F S S S K h \mathcal{K}_h K h K h \mathcal{K}_h K h ∥ x ∥ ∞ ≤ C x h \|x\|_\infty\le C_xh ∥ x ∥ ∞ ≤ C x h 0 ≤ v j ≤ C v h 2 0\le v_j\le C_vh^2 0 ≤ v j ≤ C v h 2
Normalized coordinates
Theorem 4.13 gives an eventual bound of order O ( h ) O(h) O ( h ) x x x O ( h 2 ) O(h^2) O ( h 2 ) v v v h > h ∗ h > h_* h > h ∗ Figure 1 exhibit this same scaling. And for period-two orbits near h ∗ h_* h ∗ x = O ( h − h ∗ ) x = O(h-h_*) x = O ( h − h ∗ ) v = O ( ( h − h ∗ ) 2 ) v = O((h-h_*)^2) v = O (( h − h ∗ ) 2 ) x x x v v v
Figure 1. For the matrix and parameters in (54) and (50), specified in § 6.1 below, large-h h h h h h h 2 h^2 h 2 Y ∈ { X , W } Y \in\{X,W\} Y ∈ { X , W } R Y R_Y R Y ∥ Y k − Y ‾ ∥ 2 \lVert Y_k-\overline{Y}\rVert_2 ∥ Y k − Y ∥ 2 Y ‾ \overline{Y} Y x = h X x=hX x = h X v = h 2 ( W ⊙ W ) v=h^2(W\odot W) v = h 2 ( W ⊙ W )
For the numerical experiments, we use the normalized coordinates X = x / h X=x/h X = x / h W = ⊙ v / h W=\sqrt{\odot}v/h W = ⊙ v / h X X X W W W h h h h ≥ h ∗ h\ge h_* h ≥ h ∗ G : R d × R ≥ 0 d → R d × R ≥ 0 d G:\mathbb{R}^d\times\mathbb{R}_{\ge0}^d\to\mathbb{R}^d\times\mathbb{R}_{\ge0}^d G : R d × R ≥ 0 d → R d × R ≥ 0 d (12) , with G ( X k , W k ) = ( X k + 1 , W k + 1 ) G(X_k,W_k)=(X_{k+1},W_{k+1}) G ( X k , W k ) = ( X k + 1 , W k + 1 ) F I F_I F I G I : R d × R ≥ 0 d → R d × R ≥ 0 d G_I:\mathbb{R}^d\times\mathbb{R}_{\ge0}^d\to\mathbb{R}^d\times\mathbb{R}_{\ge0}^d G I : R d × R ≥ 0 d → R d × R ≥ 0 d G I ( X k , W k ) = ( X k + 1 , W k + 1 ) G_I(X_k,W_k)=(X_{k+1},W_{k+1}) G I ( X k , W k ) = ( X k + 1 , W k + 1 )
W k + 1 = G 2 ( X k , W k ) , X k + 1 = X k − A X k ⊙ ( W k + 1 ⊙ 2 + ϵ h 2 1 . W_{k+1}=G_2(X_k,W_k),\qquad X_{k+1}=X_k-AX_k\mathbin{\odot}\sqrt{\mathstrut W_{k+1}^{\odot2}+\frac{\epsilon}{h^2}\mathbf{1}}. W k + 1 = G 2 ( X k , W k ) , X k + 1 = X k − A X k ⊙ ( W k + 1 ⊙ 2 + h 2 ϵ 1 . The normalization is a change of coordinates for both epsilon conventions, where the phase space is R d × R ≥ 0 d \mathbb{R}^d\times\mathbb{R}_{\ge0}^d R d × R ≥ 0 d
Proposition 4.14 (Conjugacy under normalization). Fix h > 0 h>0 h > 0
H h W : R d × R ≥ 0 d → R d × R ≥ 0 d , H h W ( x , v ) = ( x h , ⊙ v h ) , H_h^W:\mathbb{R}^d\times\mathbb{R}_{\ge0}^d\to\mathbb{R}^d\times\mathbb{R}_{\ge0}^d,\qquad H_h^W(x,v)=\left(\frac{x}{h},\frac{\sqrt{\odot}v}{h}\right), H h W : R d × R ≥ 0 d → R d × R ≥ 0 d , H h W ( x , v ) = ( h x , h ⊙ v ) , is a homeomorphism with inverse ( H h W ) − 1 ( X , W ) = ( h X , h 2 W ⊙ 2 ) (H_h^W)^{-1}(X,W)=(hX,h^2W^{\odot2}) ( H h W ) − 1 ( X , W ) = ( h X , h 2 W ⊙ 2 )
H h W ∘ F = G ∘ H h W , H h W ∘ F I = G I ∘ H h W . H_h^W\circ F=G\circ H_h^W,\qquad H_h^W\circ F_I=G_I\circ H_h^W. H h W ∘ F = G ∘ H h W , H h W ∘ F I = G I ∘ H h W . Equivalently, the epsilon-outside maps satisfy
G = H h W ∘ F ∘ ( H h W ) − 1 . G=H_h^W\circ F\circ(H_h^W)^{-1}. G = H h W ∘ F ∘ ( H h W ) − 1 . Consequently, each original-coordinate map and its normalized representation are topologically conjugate for every fixed h > 0 h>0 h > 0
Proof. The componentwise square-root map is a homeomorphism of R ≥ 0 d \mathbb{R}_{\ge0}^d R ≥ 0 d W ↦ W ⊙ 2 W\mapsto W^{\odot2} W ↦ W ⊙ 2 x = h X x=hX x = h X v = h 2 W ⊙ 2 v=h^2W^{\odot2} v = h 2 W ⊙ 2
Thus the normalization preserves invariant sets, periodic orbits, Lyapunov stability, and asymptotic stability. The same conclusion holds for the mapping:
H h w ( x , v ) = ( x h , v h 2 ) , ( H h w ) − 1 ( X , w ) = ( h X , h 2 w ) , H_h^w(x,v)=\left(\frac{x}{h},\frac{v}{h^2}\right),\qquad(H_h^w)^{-1}(X,w)=(hX,h^2w), H h w ( x , v ) = ( h x , h 2 v ) , ( H h w ) − 1 ( X , w ) = ( h X , h 2 w ) , for which the conjugacy is an invertible linear scaling.
Additionally, on compact invariant sets separated from v i = 0 v_i=0 v i = 0
We also compare the size of the numerically observed attractors in the original and normalized coordinates. For each sampled learning rate, we discard an initial transient and estimate a radius about the sample mean from the remaining iterates, as described in the caption. Figure 1 shows that, for large h h h h h h h 2 h^2 h 2
Period-two orbits
In the previous section, we show that every RMSProp orbit eventually enters a compact forward-invariant region. A natural next question is what kinds of dynamics occur there. We begin with period-two orbits for h > h ∗ h>h_* h > h ∗
We continue with the epsilon-outside map in normalized coordinates (c = ϵ / h c=\epsilon/h c = ϵ / h (12) . For a nontrivial period-two orbit P 1 = ( X 1 , W 1 ) ⟷ P 2 = ( X 2 , W 2 ) P_1=(X_1,W_1)\longleftrightarrow P_2=(X_2,W_2) P 1 = ( X 1 , W 1 ) ⟷ P 2 = ( X 2 , W 2 )
Lemma 5.1 (Memory invariance and antipodal symmetry). Every nontrivial period-two orbit satisfies
W 1 = W 2 = ∣ A X 1 ∣ = ∣ A X 2 ∣ and X 2 = − X 1 . (21) W_1=W_2=\lvert AX_1\rvert=\lvert AX_2\rvert\quad\text{and}\quad X_2=-X_1. \tag*{(21)} W 1 = W 2 = ∣ A X 1 ∣ = ∣ A X 2 ∣ and X 2 = − X 1 . (21) Proof. Fix a coordinate and abbreviate u = ( W 1 ) i u=(W_1)_i u = ( W 1 ) i v = ( W 2 ) i v=(W_2)_i v = ( W 2 ) i a = ∣ ( A X 1 ) i ∣ a=\lvert(AX_1)_i\rvert a = ∣( A X 1 ) i ∣ b = ∣ ( A X 2 ) i ∣ b=\lvert(AX_2)_i\rvert b = ∣( A X 2 ) i ∣
( 1 + β ) u 2 = β a 2 + b 2 , ( 1 + β ) v 2 = a 2 + β b 2 , a ( u + c ) = b ( v + c ) . (1+\beta)u^2=\beta a^2+b^2,\qquad(1+\beta)v^2=a^2+\beta b^2,\qquad a(u+c)=b(v+c). ( 1 + β ) u 2 = β a 2 + b 2 , ( 1 + β ) v 2 = a 2 + β b 2 , a ( u + c ) = b ( v + c ) . If a > b a>b a > b a u < b v au<bv a u < b v β a 4 < β b 4 \beta a^4<\beta b^4 β a 4 < β b 4 a , b a,b a , b b > a b>a b > a a = b a=b a = b u = v = a u=v=a u = v = a
Thus W 1 = W 2 = ∣ A X 1 ∣ = ∣ A X 2 ∣ W_1=W_2=\lvert AX_1\rvert=\lvert AX_2\rvert W 1 = W 2 = ∣ A X 1 ∣ = ∣ A X 2 ∣ A ( X 1 + X 2 ) ⊘ ( W 1 + c 1 ) = 0 A(X_1+X_2)\mathbin{\oslash}(W_1+c\mathbf{1})=0 A ( X 1 + X 2 ) ⊘ ( W 1 + c 1 ) = 0 A A A X 2 = − X 1 X_2=-X_1 X 2 = − X 1
Theorem 5.2 (Reduced period-two equation). Period-two cycles modulo X ∼ − X X\sim-X X ∼ − X
2 X = A X ⊘ ( ∣ A X ∣ + c 1 ) ⟺ A X = 2 c X ⊘ ( 1 − 2 ∣ X ∣ ) , c = ϵ / h . (22) 2X=AX\mathbin{\oslash}(\lvert AX\rvert+c\mathbf{1})
\quad\Longleftrightarrow\quad
AX=2cX\mathbin{\oslash}(1-2\lvert X\rvert),\qquad c=\epsilon/h. \tag*{(22)} 2 X = A X ⊘ (∣ A X ∣ + c 1 ) ⟺ A X = 2 c X ⊘ ( 1 − 2 ∣ X ∣) , c = ϵ / h . (22) Every solution of the left equation satisfies ∥ X ∥ ∞ < 1 / 2 \lVert X\rVert_\infty<1/2 ∥ X ∥ ∞ < 1/2
Proof. The preceding lemma gives P 1 = ( X , ∣ A X ∣ ) P_1=(X,\lvert AX\rvert) P 1 = ( X , ∣ A X ∣) P 2 = ( − X , ∣ A X ∣ ) P_2=(-X,\lvert AX\rvert) P 2 = ( − X , ∣ A X ∣) ∣ A X ∣ \lvert AX\rvert ∣ A X ∣ P 1 P_1 P 1 P 2 P_2 P 2 X X X − X -X − X 2 ∣ X i ∣ = ∣ ( A X ) i ∣ / ( ∣ ( A X ) i ∣ + c ) < 1 2\lvert X_i\rvert=\lvert(AX)_i\rvert/(\lvert(AX)_i\rvert+c)<1 2 ∣ X i ∣ = ∣( A X ) i ∣ / (∣( A X ) i ∣ + c ) < 1 X i X_i X i ( A X ) i (AX)_i ( A X ) i ∣ ( A X ) i ∣ = 2 c ∣ X i ∣ / ( 1 − 2 ∣ X i ∣ ) \lvert(AX)_i\rvert=2c\lvert X_i\rvert/(1-2\lvert X_i\rvert) ∣( A X ) i ∣ = 2 c ∣ X i ∣ / ( 1 − 2 ∣ X i ∣) ( − 1 / 2 , 1 / 2 ) d (-1/2,1/2)^d ( − 1/2 , 1/2 ) d X i = 0 X_i=0 X i = 0
Orbits near the stability threshold
Remark. To find period-two orbits, we use the functional
Φ c ( X ) : = 1 2 X T A X − c ∑ i = 1 d p ( X i ) , p ( t ) : = − 1 2 log ( 1 − 2 ∣ t ∣ ) − ∣ t ∣ , (23) \Phi_c(X) := \frac{1}{2}X^\mathsf{T}AX-c\sum_{i=1}^{d}p(X_i), \qquad p(t):=-\frac{1}{2}\log(1-2|t|)-|t|, \tag*{(23)} Φ c ( X ) := 2 1 X T A X − c i = 1 ∑ d p ( X i ) , p ( t ) := − 2 1 log ( 1 − 2∣ t ∣ ) − ∣ t ∣ , (23) We define it on ( − 1 / 2 , 1 / 2 ) d (-1/2,1/2)^d ( − 1/2 , 1/2 ) d
∇ Φ c ( X ) = A X − 2 c X ⊘ ( 1 − 2 ∣ X ∣ ) = 0. (24) \nabla\Phi_c(X)=AX-2cX\oslash(1-2|X|)=0. \tag*{(24)} ∇ Φ c ( X ) = A X − 2 c X ⊘ ( 1 − 2∣ X ∣ ) = 0. (24) By Theorem 5.2 , we obtain a period-two orbit from any nonzero critical point. To establish existence, it is therefore enough to find a nonzero maximizer of Φ c \Phi_c Φ c
Since Φ c \Phi_c Φ c − ∞ -\infty − ∞ h > h ∗ h>h_* h > h ∗ Φ c \Phi_c Φ c Φ c ( 0 ) = 0 \Phi_c(0)=0 Φ c ( 0 ) = 0
Theorem 5.3 (Existence threshold). A nontrivial period-two orbit exists if and only if h > h ∗ : = 2 ϵ / L A h>h_*:=2\epsilon/L_A h > h ∗ := 2 ϵ / L A 0 < c < L A / 2 0<c<L_A/2 0 < c < L A /2
Proof. If X ≠ 0 X\ne0 X = 0 X T A X = 2 c ∑ i X i 2 / ( 1 − 2 ∣ X i ∣ ) > 2 c ∥ X ∥ 2 2 X^\mathsf{T}AX=2c\sum_i X_i^2/(1-2|X_i|)>2c\lVert X\rVert_2^2 X T A X = 2 c ∑ i X i 2 / ( 1 − 2∣ X i ∣ ) > 2 c ∥ X ∥ 2 2 X T A X ≤ L A ∥ X ∥ 2 2 X^\mathsf{T}AX\le L_A\lVert X\rVert_2^2 X T A X ≤ L A ∥ X ∥ 2 2 2 c < L A 2c<L_A 2 c < L A
For sufficiency, use Φ c \Phi_c Φ c p ( t ) → ∞ p(t)\to\infty p ( t ) → ∞ ∣ t ∣ ↑ 1 / 2 |t|\mathrel{\uparrow}1/2 ∣ t ∣ ↑ 1/2 X T A X X^\mathsf{T}AX X T A X Φ c → − ∞ \Phi_c\to-\infty Φ c → − ∞ p ( t ) = t 2 + 4 3 ∣ t ∣ 3 + O ( t 4 ) p(t)=t^2+\frac{4}{3}|t|^3+O(t^4) p ( t ) = t 2 + 3 4 ∣ t ∣ 3 + O ( t 4 ) e max e_{\max} e m a x L A L_A L A A e max = L A e max Ae_{\max}=L_Ae_{\max} A e m a x = L A e m a x t e max te_{\max} t e m a x Φ c ( t e max ) = ( L A 2 − c ) t 2 − 4 c 3 ∣ t ∣ 3 ∑ i ∣ e max , i ∣ 3 + O ( t 4 ) \Phi_c(te_{\max})=(\frac{L_A}{2}-c)t^2-\frac{4c}{3}|t|^3\sum_i|e_{\max,i}|^3+O(t^4) Φ c ( t e m a x ) = ( 2 L A − c ) t 2 − 3 4 c ∣ t ∣ 3 ∑ i ∣ e m a x , i ∣ 3 + O ( t 4 ) c < L A / 2 c<L_A/2 c < L A /2 Φ c ( t e max ) > 0 = Φ c ( 0 ) \Phi_c(te_{\max})>0=\Phi_c(0) Φ c ( t e m a x ) > 0 = Φ c ( 0 ) t ≠ 0 t\ne0 t = 0 Theorem 5.2 .
Remark. To find period-two orbits numerically, we maximize Φ c \Phi_c Φ c Theorem 5.3 . We apply this method for h > h ∗ h>h_* h > h ∗ § 6.4.1 .
Theorem 5.4 (Collision with the origin). Every nonzero solution X X X ∥ X ∥ ∞ ≤ d ( L A − 2 c ) 4 c = d 2 ( h L A 2 ϵ − 1 ) \lVert X\rVert_\infty\le\frac{d(L_A-2c)}{4c}=\frac{d}{2}\left(\frac{hL_A}{2\epsilon}-1\right) ∥ X ∥ ∞ ≤ 4 c d ( L A − 2 c ) = 2 d ( 2 ϵ h L A − 1 ) h ↓ h ∗ h\downarrow h_* h ↓ h ∗ ( 0 , 0 ) (0,0) ( 0 , 0 )
Proof. Taking the inner product of (22) with X X X 4 c ∑ i ∣ X i ∣ 3 / ( 1 − 2 ∣ X i ∣ ) ≤ ( L A − 2 c ) ∥ X ∥ 2 2 4c\sum_i|X_i|^3/(1-2|X_i|)\le(L_A-2c)\lVert X\rVert_2^2 4 c ∑ i ∣ X i ∣ 3 / ( 1 − 2∣ X i ∣ ) ≤ ( L A − 2 c ) ∥ X ∥ 2 2 r = ∥ X ∥ ∞ r=\lVert X\rVert_\infty r = ∥ X ∥ ∞ r 3 r^3 r 3 ∥ X ∥ 2 2 ≤ d r 2 \lVert X\rVert_2^2\le dr^2 ∥ X ∥ 2 2 ≤ d r 2 4 c r 3 ≤ d ( L A − 2 c ) r 2 4cr^3\le d(L_A-2c)r^2 4 c r 3 ≤ d ( L A − 2 c ) r 2 c ↑ L A / 2 c\uparrow L_A/2 c ↑ L A /2 X → 0 X\to0 X → 0 W = ∣ A X ∣ → 0 W=|AX|\to0 W = ∣ A X ∣ → 0
Remark. For any sequence of period-two orbits with learning rates decreasing to h ∗ h_* h ∗ Theorem 5.4 . Whether an individual branch can be continued all the way to h ∗ h_* h ∗
Under the additional generic assumptions that L A L_A L A e max e_{\max} e m a x h > h ∗ h>h_* h > h ∗ h ∗ h_* h ∗ Φ c \Phi_c Φ c
Lemma 5.5 (Construction of the period-two branch). Let A = A T ≻ 0 A=A^\mathsf{T}\succ0 A = A T ≻ 0 L A = λ max ( A ) L_A=\lambda_{\max}(A) L A = λ m a x ( A ) e = e max e=e_{\max} e = e m a x e i ≠ 0 e_i\ne0 e i = 0 i i i ϵ > 0 \epsilon>0 ϵ > 0 β ∈ ( 0 , 1 ) \beta\in(0,1) β ∈ ( 0 , 1 )
c = ϵ / h , c ∗ = L A / 2 , h ∗ = 2 ϵ / L A , δ = L A − 2 c , σ 3 = ∑ i ∣ e i ∣ 3 . c=\epsilon/h,\qquad c_*=L_A/2,\qquad h_*=2\epsilon/L_A,\qquad\delta=L_A-2c,\qquad\sigma_3=\sum_i|e_i|^3. c = ϵ / h , c ∗ = L A /2 , h ∗ = 2 ϵ / L A , δ = L A − 2 c , σ 3 = i ∑ ∣ e i ∣ 3 . For c < c ∗ c<c_{\ast} c < c ∗ c ∗ c_{\ast} c ∗ h > h ∗ h>h_{\ast} h > h ∗ h ∗ h_{\ast} h ∗
P 1 ( c ) = ( X ( c ) , W ( c ) ) , P 2 ( c ) = ( − X ( c ) , W ( c ) ) , W ( c ) = ∣ A X ( c ) ∣ . (25) P_1(c)=(X(c),W(c)),\qquad P_2(c)=(-X(c),W(c)),\qquad W(c)=|AX(c)|. \tag*{(25)} P 1 ( c ) = ( X ( c ) , W ( c )) , P 2 ( c ) = ( − X ( c ) , W ( c )) , W ( c ) = ∣ A X ( c ) ∣. (25) with
X ( c ) = δ e 2 L A σ 3 + O ( δ 2 ) , W ( c ) = δ ∣ e ∣ 2 σ 3 + O ( δ 2 ) . X(c)=\frac{\delta e}{2L_A\sigma_3}+O(\delta^2),\qquad W(c)=\frac{\delta|e|}{2\sigma_3}+O(\delta^2). X ( c ) = 2 L A σ 3 δ e + O ( δ 2 ) , W ( c ) = 2 σ 3 δ ∣ e ∣ + O ( δ 2 ) . Moreover, W ( c ) > 0 W(c)>0 W ( c ) > 0 sign X ( c ) = sign ( A X ( c ) ) = sign e \operatorname{sign}X(c)=\operatorname{sign}(AX(c))=\operatorname{sign}e sign X ( c ) = sign ( A X ( c )) = sign e ϵ \epsilon ϵ c ↑ c ∗ c \uparrow c_{\ast} c ↑ c ∗ h ↓ h ∗ h\downarrow h_{\ast} h ↓ h ∗
Proof. Use the parametrization
X = t ( e + z ) , z ∈ e ⊥ , t > 0. (26) X=t(e+z),\qquad z\in e^{\perp},\qquad t>0. \tag*{(26)} X = t ( e + z ) , z ∈ e ⊥ , t > 0. (26) Substituting (26) into (22) and dividing by t t t
Ψ ( t , z , c ) : = A ( e + z ) − 2 c ( e + z ) ⊙ ( 1 − 2 t ∣ e + z ∣ ) = 0. (27) \Psi(t,z,c):=A(e+z)-2c(e+z)\odot(1-2t|e+z|)=0. \tag*{(27)} Ψ ( t , z , c ) := A ( e + z ) − 2 c ( e + z ) ⊙ ( 1 − 2 t ∣ e + z ∣ ) = 0. (27) Every coordinate of e e e
∥ z ∥ ∞ < 1 2 min i ∣ e i ∣ , sign ( e i + z i ) = sign e i , ∣ e + z ∣ = diag ( sign e ) ( e + z ) . \|z\|_{\infty}<\frac{1}{2}\min_i|e_i|,\qquad\operatorname{sign}(e_i+z_i)=\operatorname{sign}e_i,\qquad|e+z|=\operatorname{diag}(\operatorname{sign}e)(e+z). ∥ z ∥ ∞ < 2 1 i min ∣ e i ∣ , sign ( e i + z i ) = sign e i , ∣ e + z ∣ = diag ( sign e ) ( e + z ) . The absolute value is therefore a fixed linear map here, so it does not obstruct smoothness. For small ∣ t ∣ |t| ∣ t ∣ Ψ \Psi Ψ ( 0 , 0 , c ∗ ) (0,0,c_{\ast}) ( 0 , 0 , c ∗ )
To solve for ( z , c ) (z,c) ( z , c ) t t t D ( z , c ) Ψ ( 0 , 0 , c ∗ ) D_{(z,c)}\Psi(0,0,c_{\ast}) D ( z , c ) Ψ ( 0 , 0 , c ∗ )
Ψ ( 0 , 0 , c ∗ ) = 0 , L : = D ( z , c ) Ψ ( 0 , 0 , c ∗ ) : e ⊥ × R → R d . \Psi(0,0,c_{\ast})=0,\qquad L:=D_{(z,c)}\Psi(0,0,c_{\ast}):e^{\perp}\times\mathbb{R}\to\mathbb{R}^{d}. Ψ ( 0 , 0 , c ∗ ) = 0 , L := D ( z , c ) Ψ ( 0 , 0 , c ∗ ) : e ⊥ × R → R d . Here the total derivative in ( z , c ) (z,c) ( z , c ) t t t
D ( z , c ) Ψ ( t , z , c ) ( u , s ) = D z Ψ ( t , z , c ) u + ∂ c Ψ ( t , z , c ) s = d d τ Ψ ( t , z + τ u , c + τ s ) ∣ τ = 0 . (28) D_{(z,c)}\Psi(t,z,c)(u,s)=D_z\Psi(t,z,c)u+\partial_c\Psi(t,z,c)s=\left.\frac{d}{d\tau}\Psi(t,z+\tau u,c+\tau s)\right|_{\tau=0}. \tag*{(28)} D ( z , c ) Ψ ( t , z , c ) ( u , s ) = D z Ψ ( t , z , c ) u + ∂ c Ψ ( t , z , c ) s = d τ d Ψ ( t , z + τ u , c + τ s ) τ = 0 . (28) At t = 0 t=0 t = 0 (27) is 1 1 1 (28) gives
L ( u , s ) = d d τ [ A ( e + τ u ) − 2 ( c ∗ + τ s ) ( e + τ u ) ] ∣ τ = 0 = A u − 2 c ∗ u − 2 s e = ( A − L A I ) u − 2 s e . (29) L(u,s)=\left.\frac{d}{d\tau}\left[A(e+\tau u)-2(c_{\ast}+\tau s)(e+\tau u)\right]\right|_{\tau=0}=Au-2c_{\ast}u-2se=(A-L_A I)u-2se. \tag*{(29)} L ( u , s ) = d τ d [ A ( e + τ u ) − 2 ( c ∗ + τ s ) ( e + τ u ) ] τ = 0 = A u − 2 c ∗ u − 2 se = ( A − L A I ) u − 2 se . (29) Since A A A L A L_A L A ( A − L A I ) ∣ e ⊥ (A-L_A I)|_{e^{\perp}} ( A − L A I ) ∣ e ⊥ y ∈ R d y\in\mathbb{R}^{d} y ∈ R d
y = ( e T y ) e + ( I − e e T ) y . (30) y=(e^{\mathsf T}y)e+(I-ee^{\mathsf T})y. \tag*{(30)} y = ( e T y ) e + ( I − e e T ) y . (30) Matching the components in (30) with those in (29) gives
L ( u , s ) = y ⟺ s = − 1 2 e T y , u = [ ( A − L A I ) ∣ e ⊥ ] − 1 ( I − e e T ) y . L(u,s)=y\quad\Longleftrightarrow\quad s=-\frac{1}{2}e^{\mathsf T}y,\qquad u=\left[(A-L_A I)|_{e^{\perp}}\right]^{-1}(I-ee^{\mathsf T})y. L ( u , s ) = y ⟺ s = − 2 1 e T y , u = [ ( A − L A I ) ∣ e ⊥ ] − 1 ( I − e e T ) y . Thus the implicit-function theorem [32 Chapter 3] gives smooth z ( t ) , c ( t ) z(t),c(t) z ( t ) , c ( t ) z ( 0 ) = 0 z(0)=0 z ( 0 ) = 0 c ( 0 ) = c ∗ c(0)=c_{\ast} c ( 0 ) = c ∗ (27) along ( t , z ( t ) , c ( t ) ) (t,z(t),c(t)) ( t , z ( t ) , c ( t )) t = 0 t=0 t = 0 e T e^{\mathsf T} e T
( A − L A I ) z ′ ( 0 ) − 2 c ′ ( 0 ) e − 2 L A e ⊙ ∣ e ∣ = 0 ⟹ − 2 c ′ ( 0 ) − 2 L A ∑ i ∣ e i ∣ 3 = 0. (31) (A-L_A I)z'(0)-2c'(0)e-2L_Ae\odot|e|=0\quad\Longrightarrow\quad-2c'(0)-2L_A\sum_i|e_i|^3=0. \tag*{(31)} ( A − L A I ) z ′ ( 0 ) − 2 c ′ ( 0 ) e − 2 L A e ⊙ ∣ e ∣ = 0 ⟹ − 2 c ′ ( 0 ) − 2 L A i ∑ ∣ e i ∣ 3 = 0. (31) By (31) , c ′ ( 0 ) = − L A σ 3 < 0 c'(0)=-L_A\sigma_3<0 c ′ ( 0 ) = − L A σ 3 < 0 c ( t ) c(t) c ( t ) c ( 0 ) = L A / 2 c(0)=L_A/2 c ( 0 ) = L A /2
c ( t ) = L A 2 − L A σ 3 t + O ( t 2 ) , δ ( c ( t ) ) = L A − 2 c ( t ) = 2 L A σ 3 t + O ( t 2 ) . (32) c(t)=\frac{L_A}{2}-L_A\sigma_3t+O(t^2),\qquad\delta(c(t))=L_A-2c(t)=2L_A\sigma_3t+O(t^2). \tag*{(32)} c ( t ) = 2 L A − L A σ 3 t + O ( t 2 ) , δ ( c ( t )) = L A − 2 c ( t ) = 2 L A σ 3 t + O ( t 2 ) . (32) Inverting (32) and substituting into (26) and (25) , using z ( t ) = O ( t ) z(t)=O(t) z ( t ) = O ( t ) A e = L A e Ae=L_Ae A e = L A e
t = δ 2 L A σ 3 + O ( δ 2 ) , t=\frac{\delta}{2L_A\sigma_3}+O(\delta^2), t = 2 L A σ 3 δ + O ( δ 2 ) , X ( c ) = t ( e + z ( t ) ) = t e + O ( t 2 ) = δ e 2 L A σ 3 + O ( δ 2 ) , (33) X(c)=t(e+z(t))=te+O(t^2)=\frac{\delta e}{2L_A\sigma_3}+O(\delta^2), \tag*{(33)} X ( c ) = t ( e + z ( t )) = t e + O ( t 2 ) = 2 L A σ 3 δ e + O ( δ 2 ) , (33) W ( c ) = ∣ A X ( c ) ∣ = ∣ δ e 2 σ 3 + O ( δ 2 ) ∣ = δ ∣ e ∣ 2 σ 3 + O ( δ 2 ) . (34) W(c)=|AX(c)|=\left|\frac{\delta e}{2\sigma_3}+O(\delta^2)\right|=\frac{\delta|e|}{2\sigma_3}+O(\delta^2). \tag*{(34)} W ( c ) = ∣ A X ( c ) ∣ = 2 σ 3 δ e + O ( δ 2 ) = 2 σ 3 δ ∣ e ∣ + O ( δ 2 ) . (34) This constructs the branch for c < c ∗ c<c_* c < c ∗ h > h ∗ h>h_* h > h ∗ ϵ \epsilon ϵ e e e W ( c ) > 0 W(c)>0 W ( c ) > 0 sign X ( c ) = sign ( A X ( c ) ) = sign e \operatorname{sign}X(c)=\operatorname{sign}(AX(c))=\operatorname{sign}e sign X ( c ) = sign ( A X ( c )) = sign e δ > 0 \delta>0 δ > 0 (22) then shows that G G G P 1 ( c ) P_1(c) P 1 ( c ) P 2 ( c ) P_2(c) P 2 ( c )
We next use the Hessian of Φ c \Phi_c Φ c
Lemma 5.6 (Nondegeneracy of the period-two branch). Under the assumptions and notation of Lemma 5.5, the critical points ± X ( c ) \pm X(c) ± X ( c ) Φ c \Phi_c Φ c (23) for sufficiently small δ > 0 \delta>0 δ > 0
λ ( δ ) = − δ + O ( δ 2 ) . \lambda(\delta)=-\delta+O(\delta^2). λ ( δ ) = − δ + O ( δ 2 ) . Proof. Recall that δ = L A − 2 c \delta=L_A-2c δ = L A − 2 c H ∗ = A − L A I H_*=A-L_AI H ∗ = A − L A I B = I − ( 2 / σ 3 ) diag ( ∣ e ∣ ) B=I-(2/\sigma_3)\operatorname{diag}(|e|) B = I − ( 2/ σ 3 ) diag ( ∣ e ∣ ) (33) from Lemma 5.5 gives
( 1 − 2 ∣ X i ( c ) ∣ ) − 2 = 1 + 4 ∣ X i ( c ) ∣ + O ( δ 2 ) = 1 + 2 δ ∣ e i ∣ L A σ 3 + O ( δ 2 ) . (35) (1-2|X_i(c)|)^{-2}=1+4|X_i(c)|+O(\delta^2)=1+\frac{2\delta|e_i|}{L_A\sigma_3}+O(\delta^2). \tag*{(35)} ( 1 − 2∣ X i ( c ) ∣ ) − 2 = 1 + 4∣ X i ( c ) ∣ + O ( δ 2 ) = 1 + L A σ 3 2 δ ∣ e i ∣ + O ( δ 2 ) . (35) Substituting (35) and 2 c = L A − δ 2c=L_A-\delta 2 c = L A − δ Φ c \Phi_c Φ c (23) yields
H ( δ ) : = D 2 Φ c ( X ( c ) ) = A − 2 c diag ( ( 1 − 2 ∣ X i ( c ) ∣ ) − 2 ) = H ∗ + δ B + O ( δ 2 ) . (36) \mathcal{H}(\delta):=D^2\Phi_c(X(c))=A-2c\operatorname{diag}((1-2|X_i(c)|)^{-2})=H_*+\delta B+O(\delta^2). \tag*{(36)} H ( δ ) := D 2 Φ c ( X ( c )) = A − 2 c diag (( 1 − 2∣ X i ( c ) ∣ ) − 2 ) = H ∗ + δ B + O ( δ 2 ) . (36) The matrix H ∗ H_* H ∗ e e e
H ( δ ) ξ ( δ ) = λ ( δ ) ξ ( δ ) , ξ ( 0 ) = e , λ ( 0 ) = 0. (37) \mathcal{H}(\delta)\xi(\delta)=\lambda(\delta)\xi(\delta),\qquad\xi(0)=e,\quad\lambda(0)=0. \tag*{(37)} H ( δ ) ξ ( δ ) = λ ( δ ) ξ ( δ ) , ξ ( 0 ) = e , λ ( 0 ) = 0. (37) Differentiating (37) with respect to δ \delta δ
H ′ ( 0 ) ξ ( 0 ) + H ( 0 ) ξ ′ ( 0 ) = λ ′ ( 0 ) ξ ( 0 ) + λ ( 0 ) ξ ′ ( 0 ) . \mathcal{H}'(0)\xi(0)+\mathcal{H}(0)\xi'(0)=\lambda'(0)\xi(0)+\lambda(0)\xi'(0). H ′ ( 0 ) ξ ( 0 ) + H ( 0 ) ξ ′ ( 0 ) = λ ′ ( 0 ) ξ ( 0 ) + λ ( 0 ) ξ ′ ( 0 ) . Substitute H ′ ( 0 ) = B \mathcal{H}'(0)=B H ′ ( 0 ) = B H ( 0 ) = H ∗ \mathcal{H}(0)=H_* H ( 0 ) = H ∗ (36) , and ξ ( 0 ) = e \xi(0)=e ξ ( 0 ) = e λ ( 0 ) = 0 \lambda(0)=0 λ ( 0 ) = 0 (37) :
B e + H ∗ ξ ′ ( 0 ) = λ ′ ( 0 ) e . (38) Be+H_*\xi'(0)=\lambda'(0)e. \tag*{(38)} B e + H ∗ ξ ′ ( 0 ) = λ ′ ( 0 ) e . (38) Multiplying (38) by e ⊤ e^\top e ⊤ ξ ′ ( 0 ) \xi'(0) ξ ′ ( 0 ) e ⊤ H ∗ = 0 e^\top H_*=0 e ⊤ H ∗ = 0 e ⊤ e = 1 e^\top e=1 e ⊤ e = 1 σ 3 = ∑ i ∣ e i ∣ 3 \sigma_3=\sum_i|e_i|^3 σ 3 = ∑ i ∣ e i ∣ 3
λ ′ ( 0 ) = e ⊤ B e = ∑ i e i 2 − 2 σ 3 ∑ i e i 2 ∣ e i ∣ = 1 − 2 = − 1. \lambda'(0)=e^\top Be=\sum_i e_i^2-\frac{2}{\sigma_3}\sum_i e_i^2|e_i|=1-2=-1. λ ′ ( 0 ) = e ⊤ B e = i ∑ e i 2 − σ 3 2 i ∑ e i 2 ∣ e i ∣ = 1 − 2 = − 1. Hence λ ( δ ) = − δ + O ( δ 2 ) < 0 \lambda(\delta)=-\delta+O(\delta^2)<0 λ ( δ ) = − δ + O ( δ 2 ) < 0 ± X ( c ) \pm X(c) ± X ( c ) Φ c \Phi_c Φ c
Theorem 5.7 (Stability of the period-two branch for h ≳ h ∗ h\gtrsim h_* h ≳ h ∗ Under the assumptions of Lemma 5.5, the period-two branch { P 1 ( c ) , P 2 ( c ) } \{P_1(c),P_2(c)\} { P 1 ( c ) , P 2 ( c )} h > h ∗ h>h_* h > h ∗ h ∗ h_* h ∗ G 2 G^2 G 2
μ c r i t ( c ) = 1 − 4 ( L A − 2 c ) L A + O ( ( L A − 2 c ) 2 ) , \mu_{\mathrm{crit}}(c)=1-\frac{4(L_A-2c)}{L_A}+O((L_A-2c)^2), μ crit ( c ) = 1 − L A 4 ( L A − 2 c ) + O (( L A − 2 c ) 2 ) , and all other Floquet multipliers remain strictly inside the unit disk as h ↓ h ∗ h\downarrow h_* h ↓ h ∗ ± X ( c ) \pm X(c) ± X ( c )
Proof. Use the branch and notation from Lemma 5.5, with c c c
T ( X , W ) = ( − X , W ) , D T = T = [ − I d 0 0 I d ] , T 2 = I 2 d . T(X,W)=(-X,W), \qquad DT=T=\begin{bmatrix}-I_d & 0 \\ 0 & I_d\end{bmatrix}, \qquad T^2=I_{2d}. T ( X , W ) = ( − X , W ) , D T = T = [ − I d 0 0 I d ] , T 2 = I 2 d . The map satisfies G ∘ T = T ∘ G G\circ T=T\circ G G ∘ T = T ∘ G P 2 ( c ) = T P 1 ( c ) = G ( P 1 ( c ) ) P_2(c)=TP_1(c)=G(P_1(c)) P 2 ( c ) = T P 1 ( c ) = G ( P 1 ( c )) P 1 ( c ) P_1(c) P 1 ( c )
D G ( P 2 ( c ) ) T = T D G ( P 1 ( c ) ) , D G ( P 2 ( c ) ) = T D G ( P 1 ( c ) ) T . DG(P_2(c))T=TDG(P_1(c)), \qquad DG(P_2(c))=TDG(P_1(c))T. D G ( P 2 ( c )) T = T D G ( P 1 ( c )) , D G ( P 2 ( c )) = T D G ( P 1 ( c )) T . Thus, by the chain rule, the derivative of the second iterate is
D ( G 2 ) ( P 1 ( c ) ) = D G ( P 2 ( c ) ) D G ( P 1 ( c ) ) = T D G ( P 1 ( c ) ) T D G ( P 1 ( c ) ) = K ( c ) 2 , D(G^2)(P_1(c))=DG(P_2(c))\,DG(P_1(c))
=TDG(P_1(c))TDG(P_1(c))=K(c)^2, D ( G 2 ) ( P 1 ( c )) = D G ( P 2 ( c )) D G ( P 1 ( c )) = T D G ( P 1 ( c )) T D G ( P 1 ( c )) = K ( c ) 2 , where K ( c ) : = T D G ( P 1 ( c ) ) K(c):=TDG(P_1(c)) K ( c ) := T D G ( P 1 ( c ))
The eigenvalues of K ( c ) 2 K(c)^2 K ( c ) 2 K ( c ) K(c) K ( c ) K K K
det ( K 2 − κ 2 I ) = det ( K − κ I ) det ( K + κ I ) ⟹ spec ( K 2 ) = { κ 2 : κ ∈ spec ( K ) } \det(K^2-\kappa^2I)=\det(K-\kappa I)\det(K+\kappa I)
\quad\Longrightarrow\quad
\operatorname{spec}(K^2)=\{\kappa^2:\kappa\in\operatorname{spec}(K)\} det ( K 2 − κ 2 I ) = det ( K − κ I ) det ( K + κ I ) ⟹ spec ( K 2 ) = { κ 2 : κ ∈ spec ( K )} holds over C \mathbb{C} C K ( c ) K(c) K ( c ) K ( c ) K(c) K ( c )
We compute D G ( X , W ) DG(X,W) D G ( X , W ) G 2 ( X , W ) G_2(X,W) G 2 ( X , W )
q = A X , D = diag ( G 2 ( X , W ) ) , ζ = diag ( G 2 ( X , W ) + c 1 ) = D + c I . q=AX,\qquad D=\operatorname{diag}(G_2(X,W)),\qquad\zeta=\operatorname{diag}(G_2(X,W)+c\mathbf{1})=D+cI. q = A X , D = diag ( G 2 ( X , W )) , ζ = diag ( G 2 ( X , W ) + c 1 ) = D + c I . Differentiating the square-root memory update and then G 1 = X − q ⊘ ( G 2 + c 1 ) G_1=X-q\mathbin{\oslash}(G_2+c\mathbf{1}) G 1 = X − q ⊘ ( G 2 + c 1 )
D X G 2 = ( 1 − β ) D − 1 diag ( q ) A , D W G 2 = β D − 1 diag ( W ) , D X G 1 = I − ζ − 1 A + ζ − 2 diag ( q ) D X G 2 , D W G 1 = ζ − 2 diag ( q ) D W G 2 . (39) \begin{aligned}
D_XG_2&=(1-\beta)D^{-1}\operatorname{diag}(q)A, &
D_WG_2&=\beta D^{-1}\operatorname{diag}(W),\\
D_XG_1&=I-\zeta^{-1}A+\zeta^{-2}\operatorname{diag}(q)D_XG_2, &
D_WG_1&=\zeta^{-2}\operatorname{diag}(q)D_WG_2. \tag*{(39)}
\end{aligned} D X G 2 D X G 1 = ( 1 − β ) D − 1 diag ( q ) A , = I − ζ − 1 A + ζ − 2 diag ( q ) D X G 2 , D W G 2 D W G 1 = β D − 1 diag ( W ) , = ζ − 2 diag ( q ) D W G 2 . (39) Substituting (39) into (40), we obtain
D G ( X , W ) = [ I − ζ − 1 A + ( 1 − β ) ζ − 2 D − 1 diag ( q ⊙ 2 ) A β ζ − 2 D − 1 diag ( q ⊙ W ) ( 1 − β ) D − 1 diag ( q ) A β D − 1 diag ( W ) ] . (40) DG(X,W)=
\begin{bmatrix}
I-\zeta^{-1}A+(1-\beta)\zeta^{-2}D^{-1}\operatorname{diag}(q^{\odot2})A
&
\beta\zeta^{-2}D^{-1}\operatorname{diag}(q\odot W)\\
(1-\beta)D^{-1}\operatorname{diag}(q)A
&
\beta D^{-1}\operatorname{diag}(W)
\end{bmatrix}. \tag*{(40)} D G ( X , W ) = [ I − ζ − 1 A + ( 1 − β ) ζ − 2 D − 1 diag ( q ⊙ 2 ) A ( 1 − β ) D − 1 diag ( q ) A β ζ − 2 D − 1 diag ( q ⊙ W ) β D − 1 diag ( W ) ] . (40) At P 1 ( c ) P_1(c) P 1 ( c ) X = X ( c ) X=X(c) X = X ( c ) W = W ( c ) W=W(c) W = W ( c ) ζ c = diag ( W + c 1 ) \zeta_c=\operatorname{diag}(W+c\mathbf{1}) ζ c = diag ( W + c 1 )
G 2 ( P 1 ( c ) ) = W = ∣ A X ∣ > 0 , D − 1 diag ( A X ) = diag ( sign e ) , ζ c − 1 A X = 2 X . (41) G_2(P_1(c))=W=|AX|>0,\qquad
D^{-1}\operatorname{diag}(AX)=\operatorname{diag}(\operatorname{sign}e),\qquad
\zeta_c^{-1}AX=2X. \tag*{(41)} G 2 ( P 1 ( c )) = W = ∣ A X ∣ > 0 , D − 1 diag ( A X ) = diag ( sign e ) , ζ c − 1 A X = 2 X . (41) Using (42) in (41) yields
D G ( P 1 ( c ) ) = [ I − ζ c − 1 A + 2 ( 1 − β ) ζ c − 1 diag ( X ) diag ( sign e ) A 2 β ζ c − 1 diag ( X ) ( 1 − β ) diag ( sign e ) A β I ] . (42) DG(P_1(c))=
\begin{bmatrix}
I-\zeta_c^{-1}A+2(1-\beta)\zeta_c^{-1}\operatorname{diag}(X)\operatorname{diag}(\operatorname{sign}e)A
&
2\beta\zeta_c^{-1}\operatorname{diag}(X)\\
(1-\beta)\operatorname{diag}(\operatorname{sign}e)A
&
\beta I
\end{bmatrix}. \tag*{(42)} D G ( P 1 ( c )) = [ I − ζ c − 1 A + 2 ( 1 − β ) ζ c − 1 diag ( X ) diag ( sign e ) A ( 1 − β ) diag ( sign e ) A 2 β ζ c − 1 diag ( X ) β I ] . (42) Finally, multiplying (43) by T T T
K ( c ) = [ − I + ζ c − 1 A − 2 ( 1 − β ) ζ c − 1 diag ( X ) diag ( sign e ) A − 2 β ζ c − 1 diag ( X ) ( 1 − β ) diag ( sign e ) A β I ] . (43) K(c)=
\begin{bmatrix}
-I+\zeta_c^{-1}A-2(1-\beta)\zeta_c^{-1}\operatorname{diag}(X)\operatorname{diag}(\operatorname{sign}e)A
&
-2\beta\zeta_c^{-1}\operatorname{diag}(X)\\
(1-\beta)\operatorname{diag}(\operatorname{sign}e)A
&
\beta I
\end{bmatrix}. \tag*{(43)} K ( c ) = [ − I + ζ c − 1 A − 2 ( 1 − β ) ζ c − 1 diag ( X ) diag ( sign e ) A ( 1 − β ) diag ( sign e ) A − 2 β ζ c − 1 diag ( X ) β I ] . (43) Along the branch, X ( c ) , W ( c ) → 0 X(c),W(c)\to0 X ( c ) , W ( c ) → 0 ζ c → c ∗ I \zeta_c\to c_*I ζ c → c ∗ I c ↑ c ∗ = L A / 2 c\uparrow c_*=L_A/2 c ↑ c ∗ = L A /2 K ( c ) K(c) K ( c )
K ∗ : = lim c ↑ c ∗ K ( c ) = [ − I + 2 A / L A 0 ( 1 − β ) diag ( sign e ) A β I ] . (44) K_*:=\lim_{c\uparrow c_*}K(c)=
\begin{bmatrix}
-I+2A/L_A & 0\\
(1-\beta)\operatorname{diag}(\operatorname{sign}e)A & \beta I
\end{bmatrix}. \tag*{(44)} K ∗ := c ↑ c ∗ lim K ( c ) = [ − I + 2 A / L A ( 1 − β ) diag ( sign e ) A 0 β I ] . (44) The eigenvalues of K ∗ K_* K ∗ β \beta β 2 λ i ( A ) / L A − 1 2\lambda_i(A)/L_A-1 2 λ i ( A ) / L A − 1 1 1 1
r = ( e L A ∣ e ∣ ) , ℓ = ( e 0 ) , K ∗ r = r , ℓ T K ∗ = ℓ T , ℓ T r = 1. (45) r=\begin{pmatrix}e\\L_A|e|\end{pmatrix},\qquad
\ell=\begin{pmatrix}e\\0\end{pmatrix},\qquad
K_*r=r,\qquad\ell^{\mathsf T}K_*=\ell^{\mathsf T},\qquad\ell^{\mathsf T}r=1. \tag*{(45)} r = ( e L A ∣ e ∣ ) , ℓ = ( e 0 ) , K ∗ r = r , ℓ T K ∗ = ℓ T , ℓ T r = 1. (45) The reduced equation (22) and the branch expansions (33)–(34) give
W i + c = c 1 − 2 ∣ X i ∣ , ∑ i e i 2 ∣ X i ∣ = δ 2 L A + O ( δ 2 ) , L A c = 2 + 2 δ L A + O ( δ 2 ) . (46) W_i+c=\frac{c}{1-2\lvert X_i\rvert},\qquad\sum_i e_i^2\lvert X_i\rvert=\frac{\delta}{2L_A}+O(\delta^2),\qquad\frac{L_A}{c}=2+\frac{2\delta}{L_A}+O(\delta^2). \tag*{(46)} W i + c = 1 − 2 ∣ X i ∣ c , i ∑ e i 2 ∣ X i ∣ = 2 L A δ + O ( δ 2 ) , c L A = 2 + L A 2 δ + O ( δ 2 ) . (46) Substituting (46) into the first block row of (43) , with the eigenvectors (45) , yields
ℓ ⊤ ( K ( c ) − K ∗ ) r = L A c ∑ i e i 2 ( 1 − 2 ∣ X i ∣ ) 2 − 2 = ( 2 + 2 δ L A ) ( 1 − 2 δ L A ) − 2 + O ( δ 2 ) = − 2 δ L A + O ( δ 2 ) . (47) \ell^\top(K(c)-K_*)r=\frac{L_A}{c}\sum_i e_i^2(1-2\lvert X_i\rvert)^2-2=\left(2+\frac{2\delta}{L_A}\right)\left(1-\frac{2\delta}{L_A}\right)-2+O(\delta^2)=-\frac{2\delta}{L_A}+O(\delta^2). \tag*{(47)} ℓ ⊤ ( K ( c ) − K ∗ ) r = c L A i ∑ e i 2 ( 1 − 2 ∣ X i ∣ ) 2 − 2 = ( 2 + L A 2 δ ) ( 1 − L A 2 δ ) − 2 + O ( δ 2 ) = − L A 2 δ + O ( δ 2 ) . (47) Since K ( c ) − K ∗ = O ( δ ) K(c)-K_*=O(\delta) K ( c ) − K ∗ = O ( δ ) (47) give
κ c r i t = 1 + ℓ ⊤ ( K ( c ) − K ∗ ) r + O ( δ 2 ) = 1 − 2 δ L A + O ( δ 2 ) , μ c r i t = κ c r i t 2 = 1 − 4 δ L A + O ( δ 2 ) ∈ ( 0 , 1 ) . \kappa_{\mathrm{crit}}=1+\ell^\top(K(c)-K_*)r+O(\delta^2)=1-\frac{2\delta}{L_A}+O(\delta^2),\qquad\mu_{\mathrm{crit}}=\kappa_{\mathrm{crit}}^2=1-\frac{4\delta}{L_A}+O(\delta^2)\in(0,1). κ crit = 1 + ℓ ⊤ ( K ( c ) − K ∗ ) r + O ( δ 2 ) = 1 − L A 2 δ + O ( δ 2 ) , μ crit = κ crit 2 = 1 − L A 4 δ + O ( δ 2 ) ∈ ( 0 , 1 ) . All other Floquet multipliers remain inside the unit disk by continuity, and symmetry gives the same spectrum at P 2 P_2 P 2 (34) , W ( c ) > 0 W(c)>0 W ( c ) > 0 G G G [16 Chapter 4] therefore proves local exponential asymptotic stability.
Period-two orbits for h ≫ h ∗ h\gg h_* h ≫ h ∗
For large learning rates, we look for period-two orbits whose normalized spatial coordinates are close to the vertices of the cube [ − 1 / 2 , 1 / 2 ] d [-1/2,1/2]^d [ − 1/2 , 1/2 ] d corners and write them as s / 2 s/2 s /2 s ∈ { − 1 , 1 } d s\in\{-1,1\}^d s ∈ { − 1 , 1 } d A A A h h h
Definition 5.8 (Admissible corner). A sign vector s ∈ { − 1 , 1 } d s\in\{-1,1\}^d s ∈ { − 1 , 1 } d weakly admissible for A A A s i ( A s ) i ≥ 0 s_i(As)_i\geq0 s i ( A s ) i ≥ 0 i = 1 , … , d i=1,\ldots,d i = 1 , … , d admissible (or strictly admissible ) for A A A s i ( A s ) i > 0 s_i(As)_i>0 s i ( A s ) i > 0 i = 1 , … , d i=1,\ldots,d i = 1 , … , d A s As A s s s s
Theorem 5.9 (Corner asymptotics). If a solution branch satisfies X ( c ) → s / 2 X(c)\to s/2 X ( c ) → s /2 as c ↓ 0 c\downarrow0 c ↓ 0 then s i ( A s ) i ≥ 0 s_i(As)_i\geq0 s i ( A s ) i ≥ 0 for all i i i so s s s s s s s / 2 s/2 s /2
X s ( c ) = s 2 − c ( s ⊘ ∣ A s ∣ ) + O ( c 2 ) ⟺ X ~ s ( h ) = s 2 − ϵ h ( s ⊘ ∣ A s ∣ ) + O ( h − 2 ) ( h → ∞ ) . (48) X_s(c)=\frac{s}{2}-c(s\mathbin{\oslash}\lvert As\rvert)+O(c^2)\quad\Longleftrightarrow\quad\widetilde{X}_s(h)=\frac{s}{2}-\frac{\epsilon}{h}(s\mathbin{\oslash}\lvert As\rvert)+O(h^{-2})\quad(h\to\infty). \tag*{(48)} X s ( c ) = 2 s − c ( s ⊘ ∣ A s ∣) + O ( c 2 ) ⟺ X s ( h ) = 2 s − h ϵ ( s ⊘ ∣ A s ∣) + O ( h − 2 ) ( h → ∞ ) . (48) Proof. If X ( c ) → s / 2 X(c)\to s/2 X ( c ) → s /2 sign ( X ( c ) ) = sign ( A X ( c ) ) = s \operatorname{sign}(X(c))=\operatorname{sign}(AX(c))=s sign ( X ( c )) = sign ( A X ( c )) = s c > 0 c>0 c > 0 c ↓ 0 c\downarrow0 c ↓ 0 s i ( A s ) i ≥ 0 s_i(As)_i\geq0 s i ( A s ) i ≥ 0 s s s sign ( A X ) = s \operatorname{sign}(AX)=s sign ( A X ) = s s / 2 s/2 s /2 E ( X , c ) E(X,c) E ( X , c ) (22) . This residual is smooth with m a t h c a l E ( s / 2 , 0 ) = 0 mathcal{E}(s/2,0)=0 ma t h c a l E ( s /2 , 0 ) = 0 D X E ( s / 2 , 0 ) = 2 I d D_XE(s/2,0)=2I_d D X E ( s /2 , 0 ) = 2 I d X s ( c ) X_s(c) X s ( c ) E ( X s ( c ) , c ) = 0 E(X_s(c),c)=0 E ( X s ( c ) , c ) = 0 c = 0 c=0 c = 0 X s ′ ( 0 ) = − s ⊘ ∣ A s ∣ X_s'(0)=-s\mathbin{\oslash}\lvert As\rvert X s ′ ( 0 ) = − s ⊘ ∣ A s ∣ (48) .
Theorem 5.10 (Existence of admissible corners). Every SPD matrix A ≻ 0 A\succ0 A ≻ 0 possesses at least one admissible corner. Consequently, a corner period-two orbit exists for all sufficiently large h h h
Proof. Let s ∗ ∈ arg max s ∈ { − 1 , 1 } d s ⊤ A s s^*\in\operatorname*{arg\,max}_{s\in\{-1,1\}^d}s^\top As s ∗ ∈ arg max s ∈ { − 1 , 1 } d s ⊤ A s i i i s ^ = s ∗ − 2 s i ∗ e i \widehat{s}=s^*-2s_i^*e_i s = s ∗ − 2 s i ∗ e i 0 ≤ ( s ∗ ) ⊤ A s ∗ − s ^ ⊤ A s ^ = 4 ( s i ∗ ( A s ∗ ) i − A i i ) 0\leq(s^*)^\top As^*-\widehat{s}^\top A\widehat{s}=4(s_i^*(As^*)_i-A_{ii}) 0 ≤ ( s ∗ ) ⊤ A s ∗ − s ⊤ A s = 4 ( s i ∗ ( A s ∗ ) i − A ii ) A ≻ 0 A\succ0 A ≻ 0 A i i > 0 A_{ii}>0 A ii > 0 s i ∗ ( A s ∗ ) i ≥ A i i > 0 s_i^*(As^*)_i\geq A_{ii}>0 s i ∗ ( A s ∗ ) i ≥ A ii > 0 i i i
Under the condition, which we call diagonal dominance , where A i i > ∑ j ≠ i ∣ A i j ∣ A_{ii}>\sum_{j\ne i}\lvert A_{ij}\rvert A ii > ∑ j = i ∣ A ij ∣
s i ( A s ) i ≥ A i i − ∑ j ≠ i ∣ A i j ∣ > 0 s_i(As)_i\geq A_{ii}-\sum_{j\ne i}\lvert A_{ij}\rvert>0 s i ( A s ) i ≥ A ii − j = i ∑ ∣ A ij ∣ > 0 for every sign vector s s s 2 d 2^d 2 d 2 d − 1 2^{d-1} 2 d − 1 c > 0 c>0 c > 0
Numerical conventions and methods
To study the dynamics in the trapping region, we obtain numerical parameterizations of invariant circles and tori using the parameterization method (see § 2.5 ), then follow their branches as the learning rate varies. We study their stability through normal Lyapunov exponents and compare the results with the behavior of nearby trajectories.
Parameters, map convention, and orientation
For the three-dimensional case study, we use the normalized epsilon-outside map G G G
d = 3 , β = 0.98 , ϵ = 10 − 7 . (49) d = 3,\qquad\beta= 0.98,\qquad\epsilon= 10^{-7}. \tag*{(49)} d = 3 , β = 0.98 , ϵ = 1 0 − 7 . (49) We write Z = ( X , W ) ∈ R 6 Z = (X,W) \in\mathbb{R}^{6} Z = ( X , W ) ∈ R 6 P = G 2 P = G^2 P = G 2 G G G P P P
Unless stated otherwise, the trajectory examples and parameter surveys use the normalized initial state
X 0 = ( 0.3 , 0.4 , − 0.3 ) , W 0 = ( 0.2 , 0.2 , 0.2 ) . (50) X_0 = (0.3, 0.4, -0.3),\qquad W_0 = (0.2, 0.2, 0.2). \tag*{(50)} X 0 = ( 0.3 , 0.4 , − 0.3 ) , W 0 = ( 0.2 , 0.2 , 0.2 ) . (50) For autonomous Adam (§ 10 ), we initialize the momentum coordinates at zero. We state separately when a trajectory is continued from the final state at the preceding parameter value.
Because RMSProp rescales gradients coordinatewise, we specify both the eigenvalues and the orientation of the Hessian. We construct a positive-definite matrix from a diagonal matrix Λ = diag ( λ 1 , λ 2 , λ 3 ) \Lambda= \operatorname{diag}(\lambda_1,\lambda_2,\lambda_3) Λ = diag ( λ 1 , λ 2 , λ 3 )
A ( ϑ ) = Q ( ϑ ) Λ Q ( ϑ ) T , Q ( ϑ ) = R z ( ϑ z ) R y ( ϑ y ) R x ( ϑ x ) . (51) A(\vartheta) = Q(\vartheta)\Lambda Q(\vartheta)^{\mathsf{T}},\qquad Q(\vartheta) = R_z(\vartheta_z)R_y(\vartheta_y)R_x(\vartheta_x). \tag*{(51)} A ( ϑ ) = Q ( ϑ ) Λ Q ( ϑ ) T , Q ( ϑ ) = R z ( ϑ z ) R y ( ϑ y ) R x ( ϑ x ) . (51) Here Q Q Q [38 Chapter 2] :
R x ( α ) = ( 1 0 0 0 cos α − sin α 0 sin α cos α ) , R y ( α ) = ( cos α 0 sin α 0 1 0 − sin α 0 cos α ) , R z ( α ) = ( cos α − sin α 0 sin α cos α 0 0 0 1 ) . (52) R_x(\alpha) =
\begin{pmatrix}
1 & 0 & 0 \\
0 & \cos\alpha& -\sin\alpha\\
0 & \sin\alpha& \cos\alpha
\end{pmatrix},\qquad
R_y(\alpha) =
\begin{pmatrix}
\cos\alpha& 0 & \sin\alpha\\
0 & 1 & 0 \\
-\sin\alpha& 0 & \cos\alpha
\end{pmatrix},\qquad
R_z(\alpha) =
\begin{pmatrix}
\cos\alpha& -\sin\alpha& 0 \\
\sin\alpha& \cos\alpha& 0 \\
0 & 0 & 1
\end{pmatrix}. \tag*{(52)} R x ( α ) = 1 0 0 0 cos α sin α 0 − sin α cos α , R y ( α ) = cos α 0 − sin α 0 1 0 sin α 0 cos α , R z ( α ) = cos α sin α 0 − sin α cos α 0 0 0 1 . (52) For the main case study, we use a fixed approximation to the matrix in (51) . We take ( ϑ x , ϑ y , ϑ z ) = ( 60 ∘ , − 60 ∘ , − 143.5 ∘ ) (\vartheta_x,\vartheta_y,\vartheta_z) = (60^\circ,-60^\circ,-143.5^\circ) ( ϑ x , ϑ y , ϑ z ) = ( 6 0 ∘ , − 6 0 ∘ , − 143. 5 ∘ ) Λ = diag ( 0.5 , 0.8 , 1.18 ) \Lambda= \operatorname{diag}(0.5,0.8,1.18) Λ = diag ( 0.5 , 0.8 , 1.18 ) A A A
A ≈ ( 0.7621 − 0.0964 0.0886 − 0.0964 1.1191 0.1679 0.0886 0.1679 0.5988 ) . (53) A \approx
\begin{pmatrix}
0.7621 & -0.0964 & 0.0886 \\
-0.0964 & 1.1191 & 0.1679 \\
0.0886 & 0.1679 & 0.5988
\end{pmatrix}. \tag*{(53)} A ≈ 0.7621 − 0.0964 0.0886 − 0.0964 1.1191 0.1679 0.0886 0.1679 0.5988 . (53) This gives h ∗ ≈ 1.6949 × 10 − 7 h_* \approx1.6949 \times10^{-7} h ∗ ≈ 1.6949 × 1 0 − 7 (2) . Finally, in the orientation survey of § 8.1 , we vary the angles in (51) while keeping the eigenvalues fixed.
Computation and continuation of invariant circles and tori
We numerically compute Fourier series expansions for attracting invariant circles and invariant tori using the parameterization method.
In the numerical examples below, we repeatedly find invariant sets in disjoint pairs (chaotic attractors, circles, periodic points), with iterates alternating between the two, as illustrated in Figure 4 and Figure 5 . We call a pair of disjoint circles exchanged by G G G period-two invariant circle and parameterize its two components together. Here period two refers to the exchange of components, independently of the motion within each component.
Figure 4. The appearance of the primary invariant circle for the matrix (53) and parameters (49) in the X 1 X_1 X 1 X 3 X_3 X 3 h = 5.4946 h = 5.4946 h = 5.4946 5.5793 5.5793 5.5793 6.0431 6.0431 6.0431 10 − 7 10^{-7} 1 0 − 7 C N S C_{\mathrm{NS}} C NS
Figure 5. Four distinct period-two circles at h = 1.0514 × 10 − 6 h = 1.0514 \times10^{-6} h = 1.0514 × 1 0 − 6 (53) and parameters (49) . Each panel shows both circle components in the X 1 X_1 X 1 X 2 X_2 X 2 G G G C 0 C_0 C 0 C 1 C_1 C 1 C 2 C_2 C 2 C N S C_{\mathrm{NS}} C NS
The two components K 1 , K 2 : T → R 6 K_1,K_2:\mathbb{T}\to\mathbb{R}^6 K 1 , K 2 : T → R 6 T = R / ( 2 π Z ) \mathbb{T}=\mathbb{R}/(2\pi\mathbb{Z}) T = R / ( 2 π Z )
G ( K 1 ( θ ) ) = K 2 ( θ + ω ) , G ( K 2 ( θ ) ) = K 1 ( θ + ω ) . (54) G(K_1(\theta))=K_2(\theta+\omega),\qquad G(K_2(\theta))=K_1(\theta+\omega). \tag*{(54)} G ( K 1 ( θ )) = K 2 ( θ + ω ) , G ( K 2 ( θ )) = K 1 ( θ + ω ) . (54) Thus each component is invariant under G 2 G^2 G 2
G 2 ( K j ( θ ) ) = K j ( θ + 2 ω ) , ρ P = ω / π ( m o d 1 ) . (55) G^2(K_j(\theta))=K_j(\theta+2\omega),\qquad\rho_P=\omega/\pi\pmod{1}. \tag*{(55)} G 2 ( K j ( θ )) = K j ( θ + 2 ω ) , ρ P = ω / π ( mod 1 ) . (55) Here ω \omega ω (54) . We define the rotation number ρ P \rho_P ρ P P P P P = G 2 P=G^2 P = G 2 ( − 1 / 2 , 1 / 2 ] (-1/2,1/2] ( − 1/2 , 1/2 ]
We solve the above invariance equations, together with phase conditions, for the Fourier coefficients and rotation number ρ P \rho_P ρ P
Since the parameterization method uses an iterative scheme to solve a functional equation describing the invariant circle, an essential step is obtaining a good enough initial approximation of the rotation number and Fourier coefficients. For this we exploit the Birkhoff averaging method of Das, Saiki, Sander, and Yorke [11 , 12 , 13 ] (see also [43 ] ). In this sense, our work extends the approach of Blessing and Mireles James [2 ] to dissipative maps in higher dimensions. We then improve these approximations with Newton’s method and follow the branches by fixed-parameter or pseudo-arclength continuation [27 ] .
Another advantage of solving a functional equation is that, by treating the learning rate as a bifurcation parameter, we follow attracting invariant circles through the loss of stability via numerical continuation. Section C.2 gives the continuation procedure and accuracy checks.
Lyapunov exponents and orbit classification
We estimate the leading finite-time Lyapunov exponent (FTLE) λ ^ 1 \hat{\lambda}_1 λ ^ 1 [1 ] , with an unweighted average over N N N λ ^ ⊥ , i \hat{\lambda}_{\perp,i} λ ^ ⊥ , i [24 ] . Appendix A gives their construction and the weighted averages that we use for the circle and torus stability comparisons in § 7.2 and § 7.3 .
We report full-state Lyapunov exponent estimates per iterate of G G G G β 1 G_{\beta_1} G β 1 P = G 2 P=G^2 P = G 2 G G G Table 4 in Appendix C gives the full-state estimates and the numbers of discarded and measured iterations for the individual trajectory examples.
Table 4. Leading full-state Lyapunov exponent estimates for the individual examples.
Recurrence tests
For N N N Z 0 , … , Z N − 1 Z_0,\ldots,Z_{N-1} Z 0 , … , Z N − 1 normalized recurrence error using the final m m m
e p = median N − m + p ≤ k < N ∥ Z k − Z k − p ∥ 2 max { 1 , max N − m ≤ k < N ∥ Z k ∥ 2 } , p = 1 , 2. (56) e_p=\frac{\operatorname{median}_{N-m+p\le k<N}\lVert Z_k-Z_{k-p}\rVert_2}{\max\left\{1,\max_{N-m\le k<N}\lVert Z_k\rVert_2\right\}},\qquad p=1,2. \tag*{(56)} e p = max { 1 , max N − m ≤ k < N ∥ Z k ∥ 2 } median N − m + p ≤ k < N ∥ Z k − Z k − p ∥ 2 , p = 1 , 2. (56) For example, small e 2 e_2 e 2 e 1 e_1 e 1 m = 128 m=128 m = 128 m = 1,024 m=1{,}024 m = 1 , 024
Criteria used in the parameter surveys
We group invariant circles, which are one-dimensional tori, and higher-dimensional tori in a single torus category for the parameter surveys. We first screen for fixed-point and period-two behavior using recurrence, then classify the remaining trajectories using the largest finite-time Lyapunov exponent. For the parameter surveys of § 8.1 , 8.2, and § 10 , we use the following four classification criteria, which do not cover every finite-window outcome:
Fixed point: e 1 ≤ 10 − 6 e_1 \le10^{-6} e 1 ≤ 1 0 − 6
Period two: e 1 > 10 − 6 e_1 > 10^{-6} e 1 > 1 0 − 6 e 2 ≤ 10 − 6 e_2 \le10^{-6} e 2 ≤ 1 0 − 6
Torus: ∣ λ ^ 1 ∣ ≤ 2 × 10 − 4 |\hat{\lambda}_1| \le2 \times10^{-4} ∣ λ ^ 1 ∣ ≤ 2 × 1 0 − 4 e 1 , e 2 > 10 − 6 e_1,e_2 > 10^{-6} e 1 , e 2 > 1 0 − 6
Chaotic: λ ^ 1 > 2 × 10 − 4 \hat{\lambda}_1 > 2 \times10^{-4} λ ^ 1 > 2 × 1 0 − 4 e 1 , e 2 > 10 − 6 e_1,e_2 > 10^{-6} e 1 , e 2 > 1 0 − 6
For each survey, we discard 50,000 iterates of G G G G β 1 G_{\beta_1} G β 1 m m m
Remark 6.1. We use fairly coarse classification criteria for our parameter surveys. The recurrence checks do not systematically identify longer-period orbits, and the Lyapunov exponent cutoff limits our ability to distinguish quasiperiodic motion from weakly chaotic trajectories. We flag samples in the torus class with estimated Lyapunov exponents above 10 − 5 10^{-5} 1 0 − 5
Computation of period-two orbits
The analytical results suggest two ways to find period-two orbits numerically. Near h ∗ h_* h ∗ Theorem 5.3 to search for them by optimization. For h ≫ h ∗ h \gg h_* h ≫ h ∗ Theorem 5.9 to obtain initial approximations, then continue the branches toward smaller learning rates. For the corner continuation below, we use the parameters in § 6.1 .
Finding period-two orbits near the threshold
For h > h ∗ h > h_* h > h ∗ Φ c \Phi_c Φ c (23) on ( − 1 / 2 , 1 / 2 ) d (-1/2,1/2)^d ( − 1/2 , 1/2 ) d A A A λ max ( A ) \lambda_{\max}(A) λ m a x ( A ) Theorem 5.7 . Starting from these vectors, we minimize − Φ c -\Phi_c − Φ c [4 , 47 ] , with bounds strictly inside the open hypercube. We refine the resulting critical points with Newton’s method applied to (24) , then recover the orbit points from (21) . After identifying X X X − X -X − X
Continuation from admissible corners
The corners are s / 2 ∈ 1 2 { − 1 , 1 } d s/2 \in\frac{1}{2}\{-1,1\}^d s /2 ∈ 2 1 { − 1 , 1 } d s ∈ { − 1 , 1 } d s \in\{-1,1\}^d s ∈ { − 1 , 1 } d s s s − s -s − s s 1 = ( + + + ) s_1=(+++) s 1 = ( + + + ) s 2 = ( + + − ) s_2=(++-) s 2 = ( + + − ) s 3 = ( + − + ) s_3=(+-+) s 3 = ( + − + ) s 4 = ( + − − ) s_4=(+--) s 4 = ( + − − ) A A A § 6.1 . For each representative corner, we initialize at large h h h (48) and track branches toward smaller h h h [27 ] . Along each path, we compute the Floquet multipliers to determine orbital stability.
We plot the continued branches against the normalized learning rate
ρ = h / h ∗ . (57) \rho= h/h_*. \tag*{(57)} ρ = h / h ∗ . (57) For the matrix A A A § 6.1 , branch s 4 s_4 s 4 s 1 , s 2 , s 3 s_1,s_2,s_3 s 1 , s 2 , s 3 ρ ≈ 2.18839 , 13.35115 , 4.41042 \rho\approx2.18839, 13.35115, 4.41042 ρ ≈ 2.18839 , 13.35115 , 4.41042
Numerical bifurcation analysis
We now return to the question raised in the introduction: what kinds of dynamics occur in the trapping region? We study the normalized RMSProp map with the matrix (54) and parameters (50), using P = G 2 P=G^2 P = G 2 Figure 3 gives an overview as we increase the learning rate. We then examine three features in more detail: coexistence of invariant circles, a secondary Neimark–Sacker bifurcation, and a circle-doubling sequence. We also find intervals of learning rates for which the computed trajectories have positive finite-time Lyapunov exponents. Appendix § 16 collects additional sampled transitions and period-two circles. Figure 21 and Figure 22 show further projections of the coexisting sets, and Appendix § 18 gives additional trajectories.
!Bifurcation diagram with panels (a) Bifurcation diagram, (b) Leading full-state exponent, (c) h ∈ [ 0.9 , 1.4 ] h \in[0.9,1.4] h ∈ [ 0.9 , 1.4 ] h ∈ [ 1.8 , 2.0 ] h \in[1.8,2.0] h ∈ [ 1.8 , 2.0 ] X 1 X_1 X 1 G G G h ( 10 − 6 ) h(10^{-6}) h ( 1 0 − 6 ) bbox=0.18,0.10,0.87,0.56
Figure 3. Bifurcation diagram for the matrix (53) and parameters (49) . For the upper panel, we begin with X 0 = W 0 = 10 − 81 X_0 = W_0 = 10^{-81} X 0 = W 0 = 1 0 − 81 10 − 8 10^{-8} 1 0 − 8 G G G X 1 X_1 X 1 G G G (50) , so they may approach different attractors from those in the upper panel when several coexist. The bottom panels enlarge h ∈ [ 0.9 , 1.4 ] × 10 − 6 h \in[0.9,1.4] \times10^{-6} h ∈ [ 0.9 , 1.4 ] × 1 0 − 6 h ∈ [ 1.8 , 2.0 ] × 10 − 6 h \in[1.8,2.0] \times10^{-6} h ∈ [ 1.8 , 2.0 ] × 1 0 − 6 h ∗ h_\ast h ∗
Figure 21. Four distinct period-two circles at h = 1.0514 × 10 − 6 h = 1.0514 \times10^{-6} h = 1.0514 × 1 0 − 6 X 1 X_1 X 1 X 2 X_2 X 2 X 1 X_1 X 1 X 3 X_3 X 3 X 2 X_2 X 2 X 3 X_3 X 3 W 1 W_1 W 1 W 2 W_2 W 2 W 1 W_1 W 1 W 3 W_3 W 3 W 2 W_2 W 2 W 3 W_3 W 3 C 0 C_0 C 0 C 1 C_1 C 1 C 2 C_2 C 2 C N S C_{\mathrm{NS}} C NS
Figure 22. Four period-two circles and two invariant-torus pairs coexist at h = 1.0514 × 10 − 6 h = 1.0514 \times10^{-6} h = 1.0514 × 1 0 − 6 Figure 6 shows their cross-sections.
Figure 6. Attracting and saddle two-tori for the matrix (53) and parameters (49) at h = 1.0514 × 10 − 6 h = 1.0514 \times10^{-6} h = 1.0514 × 1 0 − 6 P P P X 1 X_1 X 1 X 2 X_2 X 2 X 3 X_3 X 3 G G G X 1 = 0.18905212 X_1 = 0.18905212 X 1 = 0.18905212 X 2 > − 0.5 X_2 > -0.5 X 2 > − 0.5 X 1 X_1 X 1 θ 1 \theta_1 θ 1
Figure 2. We begin near the corners s / 2 ∈ 1 2 { − 1 , 1 } d s/2 \in\frac{1}{2}\{-1,1\}^{d} s /2 ∈ 2 1 { − 1 , 1 } d h h h (53) and parameters (49) . The upper and lower rows show X X X W W W ρ = h / h ∗ \rho= h/h_* ρ = h / h ∗ G 2 G^2 G 2 X 2 X_2 X 2 ρ = 2.19 \rho= 2.19 ρ = 2.19 X 1 X_1 X 1 X 2 X_2 X 2
The primary invariant circle
As we increase the learning rate, we find a conjugate pair of Floquet multipliers of the period-two orbit crossing the unit circle near h c = 5.5222246 × 10 − 7 h_c = 5.5222246 \times10^{-7} h c = 5.5222246 × 1 0 − 7 h = 5.4946135 × 10 − 7 h = 5.4946135 \times10^{-7} h = 5.4946135 × 1 0 − 7 h = 5.5792602 × 10 − 7 h = 5.5792602 \times10^{-7} h = 5.5792602 × 1 0 − 7 C N S C_{\mathrm{NS}} C NS Figure 4 .
We call a circle component q q q q q q h h h C N S C_{\mathrm{NS}} C NS h = 0.99245872 × 10 − 6 h = 0.99245872 \times10^{-6} h = 0.99245872 × 1 0 − 6 h = 1.0514 × 10 − 6 h = 1.0514 \times10^{-6} h = 1.0514 × 1 0 − 6 C N S C_{\mathrm{NS}} C NS
Multistability
At h = 1.0514 × 10 − 6 h = 1.0514 \times10^{-6} h = 1.0514 × 1 0 − 6 C 0 C_0 C 0 C 1 C_1 C 1 C 2 C_2 C 2 C N S C_{\mathrm{NS}} C NS Table 1 gives their rotation numbers and leading normal Lyapunov exponents. For each of the three attracting circles, we find trajectories from nearby initial conditions that remain close for five million iterations of P P P
Table 1. Four period-two circles for A A A h = 1.0514 × 10 − 6 h = 1.0514 \times10^{-6} h = 1.0514 × 1 0 − 6 P = G 2 P = G^2 P = G 2 10 − 7 10^{-7} 1 0 − 7 P P P
A secondary Neimark–Sacker bifurcation
Visual evidence suggests a secondary Neimark–Sacker bifurcation of C 0 C_0 C 0 h = 1.05147235 × 10 − 6 h = 1.05147235 \times10^{-6} h = 1.05147235 × 1 0 − 6 h = 1.05147230 × 10 − 6 h = 1.05147230 \times10^{-6} h = 1.05147230 × 1 0 − 6 h = 1.05147240 × 10 − 6 h = 1.05147240 \times10^{-6} h = 1.05147240 × 1 0 − 6 C 0 C_0 C 0 h = 1.0514 × 10 − 6 h = 1.0514 \times10^{-6} h = 1.0514 × 1 0 − 6 C 0 C_0 C 0 C 0 C_0 C 0
Chaotic dynamics
At h = 1.06 × 10 − 6 h = 1.06 \times10^{-6} h = 1.06 × 1 0 − 6 Figure 8 shows the trajectory from the initial state (50) after discarding 250,000 iterations of P P P
Figure 8. Finite-window chaotic dynamics for the matrix (53) and parameters (49) at h = 1.06 × 10 − 6 h = 1.06 \times10^{-6} h = 1.06 × 1 0 − 6 X 1 X_1 X 1 X 2 X_2 X 2 X 3 X_3 X 3 W 1 W_1 W 1 W 2 W_2 W 2 W 3 W_3 W 3 P P P (50) .
A circle-doubling sequence
Near h = 1.95 × 10 − 6 h = 1.95 \times10^{-6} h = 1.95 × 1 0 − 6 Figure 9 . Here we again call a circle component q q q q q q P P P h h h Figure 3 . The observed sequence suggests a torus-doubling route to chaos, although we have not identified the individual bifurcations.
Figure 9. Circle doubling for the matrix (54) and parameters (50), shown in the X 1 X_1 X 1 X 2 X_2 X 2 P = G 2 P = G^2 P = G 2 ( q , h ) = ( 3 , 1.950 × 10 − 6 ) , ( 6 , 1.953 × 10 − 6 ) , ( 12 , 1.960 × 10 − 6 ) (q,h) = (3, 1.950 \times10^{-6}), (6, 1.953 \times10^{-6}), (12, 1.960 \times10^{-6}) ( q , h ) = ( 3 , 1.950 × 1 0 − 6 ) , ( 6 , 1.953 × 1 0 − 6 ) , ( 12 , 1.960 × 1 0 − 6 ) ( 24 , 1.9615 × 10 − 6 ) (24, 1.9615 \times10^{-6}) ( 24 , 1.9615 × 1 0 − 6 ) P P P
Figure 7. The saddle torus and C 0 C_0 C 0 (53) and parameters (49) at h = 1.0514 × 10 − 6 h = 1.0514 \times10^{-6} h = 1.0514 × 1 0 − 6 X 1 = 0.18905212 X_1 = 0.18905212 X 1 = 0.18905212 W 1 W_1 W 1 W 3 W_3 W 3
Survey of parameters
We now vary the Hessian orientation, memory decay, and dimension to examine their effects on the dynamics. We then test a simple learning-rate reduction schedule.
Dependence on the Hessian orientation
We use a 6 ∘ 6^\circ 6 ∘ ( 0.5 , 0.8 , 1.18 ) (0.5, 0.8, 1.18) ( 0.5 , 0.8 , 1.18 ) β = 0.98 \beta= 0.98 β = 0.98 ϵ = 10 − 7 \epsilon= 10^{-7} ϵ = 1 0 − 7 h = 3 × 10 − 6 h = 3 \times10^{-6} h = 3 × 1 0 − 6 § 6.3.2 . Section C.4 gives the angular ranges. The Euler-angle representation of S O ( 3 ) SO(3) S O ( 3 ) [34 ] .
Using the finite-window criteria in (58), we classify 103,248 trajectories as period two, 6,043 as torus motion, and 2,309 as chaotic. The non-periodic cases in Figure 10 form narrow bands. Thus, even with the eigenvalues fixed, changing the Hessian orientation changes the observed motion.
Figure 10. We find different types of motion as we change the Hessian orientation while keeping its spectrum fixed. We use a uniform ZYX Euler-angle grid of 111,600 triples at h = 3 × 10 − 6 h = 3 \times10^{-6} h = 3 × 1 0 − 6 β \beta β ϵ \epsilon ϵ ϑ y ∈ { − 90 ∘ , − 60 ∘ , − 30 ∘ , 0 ∘ , 30 ∘ , 60 ∘ , 90 ∘ } \vartheta_y \in\{-90^\circ, -60^\circ, -30^\circ, 0^\circ, 30^\circ, 60^\circ, 90^\circ\} ϑ y ∈ { − 9 0 ∘ , − 6 0 ∘ , − 3 0 ∘ , 0 ∘ , 3 0 ∘ , 6 0 ∘ , 9 0 ∘ } § 6.3.2 for each trajectory. Blue, orange, and red denote the finite-window period-two, torus, and chaotic classifications, respectively.
Dependence on memory decay beta
We vary h h h β \beta β A A A 196 × 161 196 \times161 196 × 161 Figure 11 and the protocol of § 6.3.2 . Using these criteria, we classify 806 trajectories as approaching the origin, 28,089 as period two, 1,833 as torus motion, and 825 as chaotic. We examine the remaining three samples further because they fit none of the criteria. In longer computations, we find periodic motion with periods 108, 214, and 1,814 iterates of G G G β \beta β
Figure 11. The ( h , β ) (h,\beta) ( h , β ) A A A (53) and ϵ \epsilon ϵ (49) . The grid is specified in Section C.4. The dashed line is h = h ∗ h=h^{*} h = h ∗ G G G
Dimension
To examine how the dynamics change with dimension, we consider examples in dimensions one, two, four, and five.
Examples in dimensions one and two
In the one-dimensional computations, we find convergence to the origin below h ∗ h_* h ∗ h ∗ h_* h ∗
For d = 2 d = 2 d = 2 A 2 = Q ( 0.84975 ) diag ( 0.5 , 0.8 ) Q ( 0.84975 ) ⊤ A_2 = Q(0.84975)\operatorname{diag}(0.5, 0.8)Q(0.84975)^\top A 2 = Q ( 0.84975 ) diag ( 0.5 , 0.8 ) Q ( 0.84975 ) ⊤ β = 0.5 \beta= 0.5 β = 0.5
A 2 ≈ [ 0.669252296871 − 0.148759366311 − 0.148759366311 0.630747703129 ] , Q ( θ ) = [ cos θ − sin θ sin θ cos θ ] . A_2 \approx\begin{bmatrix} 0.669252296871 & -0.148759366311 \\ -0.148759366311 & 0.630747703129 \end{bmatrix}, \qquad Q(\theta) = \begin{bmatrix} \cos\theta& -\sin\theta\\ \sin\theta& \cos\theta\end{bmatrix}. A 2 ≈ [ 0.669252296871 − 0.148759366311 − 0.148759366311 0.630747703129 ] , Q ( θ ) = [ cos θ sin θ − sin θ cos θ ] . Here Q ( θ ) Q(\theta) Q ( θ ) θ \theta θ
Equations (2) and (57) , with λ max = 0.8 \lambda_{\max} = 0.8 λ m a x = 0.8 ϵ = 10 − 7 \epsilon= 10^{-7} ϵ = 1 0 − 7
h ∗ = 2.5 × 10 − 7 . h_* = 2.5 \times10^{-7}. h ∗ = 2.5 × 1 0 − 7 . The sweep contains 2,259 uniformly spaced values of ρ ∈ [ 0.2 , 30 ] \rho\in[0.2, 30] ρ ∈ [ 0.2 , 30 ] h ∈ [ 0.05 , 7.5 ] × 10 − 6 h \in[0.05, 7.5] \times10^{-6} h ∈ [ 0.05 , 7.5 ] × 1 0 − 6 Δ ρ ≈ 0.01320 \Delta\rho\approx0.01320 Δ ρ ≈ 0.01320
X 0 ( h ) = ϵ h A 2 − 1 b 2 , W 0 ( h ) = ϵ h a 2 , X_0(h) = \frac{\epsilon}{h}A_2^{-1}b_2, \qquad W_0(h) = \frac{\epsilon}{h}a_2, X 0 ( h ) = h ϵ A 2 − 1 b 2 , W 0 ( h ) = h ϵ a 2 , where b 2 = ( 0.3 , 0.4 ) ⊤ b_2 = (0.3, 0.4)^\top b 2 = ( 0.3 , 0.4 ) ⊤ a 2 = ( 0.2 , 0.2 ) ⊤ a_2 = (0.2, 0.2)^\top a 2 = ( 0.2 , 0.2 ) ⊤
In the two-dimensional example, we find higher-period cycles and coexistence of different types of motion. At ρ = 15 \rho= 15 ρ = 15 Figure 12 shows a stable period-two orbit together with a trajectory whose geometry and positive finite-time Lyapunov exponent suggest chaotic motion.
Figure 12. Epsilon-outside RMSProp in dimension two with A 2 A_2 A 2 β = 0.5 \beta= 0.5 β = 0.5 ϵ = 10 − 7 \epsilon= 10^{-7} ϵ = 1 0 − 7 ρ ∈ [ 0.2 , 30 ] \rho\in[0.2, 30] ρ ∈ [ 0.2 , 30 ] X 1 X_1 X 1 G G G ρ = 15 \rho= 15 ρ = 15 h = 3.75 × 10 − 6 h = 3.75 \times10^{-6} h = 3.75 × 1 0 − 6 X 0 = A − 1 ( 0.3 , 0.4 ) ⊤ / 37.5 X_0 = A^{-1}(0.3, 0.4)^\top/37.5 X 0 = A − 1 ( 0.3 , 0.4 ) ⊤ /37.5 W 0 = ( 0.2 , 0.2 ) ⊤ / 37.5 W_0 = (0.2, 0.2)^\top/37.5 W 0 = ( 0.2 , 0.2 ) ⊤ /37.5 X 0 = ( 0.4 , − 0.4 ) ⊤ X_0 = (0.4, -0.4)^\top X 0 = ( 0.4 , − 0.4 ) ⊤ W 0 = ( 0.3 , 0.3 ) ⊤ W_0 = (0.3, 0.3)^\top W 0 = ( 0.3 , 0.3 ) ⊤ G G G P 1 P_1 P 1 P 2 P_2 P 2
We continue the period-two branch emerging from the origin over ρ ∈ [ 11.2 , 20 ] \rho\in[11.2, 20] ρ ∈ [ 11.2 , 20 ]
Coupled extensions in dimensions four and five
We test whether the three-dimensional behavior persists when we add one or two coordinates and couple them to the original variables through
A d ( η ) = [ A η B d η B d ⊤ D d ] , η ≥ 0. A_d(\eta) =
\begin{bmatrix}
A & \eta B_d \\
\eta B_d^\top& D_d
\end{bmatrix},
\qquad\eta\ge0. A d ( η ) = [ A η B d ⊤ η B d D d ] , η ≥ 0. At η = 0 \eta= 0 η = 0 A ⊕ D d A \oplus D_d A ⊕ D d X 4 = ⋯ = X d = W 4 = ⋯ = W d = 0 X_4 = \cdots= X_d = W_4 = \cdots= W_d = 0 X 4 = ⋯ = X d = W 4 = ⋯ = W d = 0 d = 4 d = 4 d = 4
D 4 = [ 0.15 ] , B 4 = 1 0.40 2 + 0.70 2 + 0.55 2 [ 0.40 − 0.70 0.55 ] . D_4 = [0.15], \qquad B_4 = \frac{1}{\sqrt{0.40^2 + 0.70^2 + 0.55^2}}
\begin{bmatrix}
0.40 \\
-0.70 \\
0.55
\end{bmatrix}. D 4 = [ 0.15 ] , B 4 = 0.4 0 2 + 0.7 0 2 + 0.5 5 2 1 0.40 − 0.70 0.55 . For d = 5 d = 5 d = 5 D 5 = diag ( 0.15 , 0.12 ) D_5 = \operatorname{diag}(0.15, 0.12) D 5 = diag ( 0.15 , 0.12 )
B 5 = [ 0.50 / ν 1 − 0.20 / ν 2 − 0.35 / ν 1 0.75 / ν 2 0.60 / ν 1 0.45 / ν 2 ] , ν 1 = 0.50 2 + 0.35 2 + 0.60 2 , ν 2 = 0.20 2 + 0.75 2 + 0.45 2 . B_5 =
\begin{bmatrix}
0.50/\nu_1 & -0.20/\nu_2 \\
-0.35/\nu_1 & 0.75/\nu_2 \\
0.60/\nu_1 & 0.45/\nu_2
\end{bmatrix},
\qquad
\nu_1 = \sqrt{0.50^2 + 0.35^2 + 0.60^2}, \qquad
\nu_2 = \sqrt{0.20^2 + 0.75^2 + 0.45^2}. B 5 = 0.50/ ν 1 − 0.35/ ν 1 0.60/ ν 1 − 0.20/ ν 2 0.75/ ν 2 0.45/ ν 2 , ν 1 = 0.5 0 2 + 0.3 5 2 + 0.6 0 2 , ν 2 = 0.2 0 2 + 0.7 5 2 + 0.4 5 2 . We use η = 0.10 \eta= 0.10 η = 0.10 A d ( 0.10 ) A_d(0.10) A d ( 0.10 )
To obtain the bifurcation diagrams in Figure 13 , we initialize a trajectory at each learning rate using the first three spatial and memory coordinates from (51) and perturbations in the additional coordinates:
Figure 13. Bifurcation diagrams for d = 4 d = 4 d = 4 d = 5 d = 5 d = 5 η = 0.10 \eta= 0.10 η = 0.10 β = 0.98 \beta= 0.98 β = 0.98 ϵ = 10 − 7 \epsilon= 10^{-7} ϵ = 1 0 − 7 2 ϵ / λ max ( A d ) 2\epsilon/\lambda_{\max}(A_d) 2 ϵ / λ m a x ( A d ) G G G
X ⊥ ( 0 ) = 10 − 5 ( d = 4 ) , X ⊥ ( 0 ) = 10 − 5 ( 1 , − 1 / 2 ) ( d = 5 ) , W ⊥ ( 0 ) = 10 − 5 1. X_{\perp}^{(0)} = 10^{-5} \quad(d = 4), \qquad
X_{\perp}^{(0)} = 10^{-5}(1,-1/2) \quad(d = 5), \qquad
W_{\perp}^{(0)} = 10^{-5}\mathbf{1}. X ⊥ ( 0 ) = 1 0 − 5 ( d = 4 ) , X ⊥ ( 0 ) = 1 0 − 5 ( 1 , − 1/2 ) ( d = 5 ) , W ⊥ ( 0 ) = 1 0 − 5 1 . With this coupling, the trajectories at h = 1.06 × 10 − 6 h = 1.06 \times10^{-6} h = 1.06 × 1 0 − 6
Adaptive learning rates
In practice, we want to reduce the loss while using as few iterates as possible. For fixed h ≥ h ∗ h \ge h_* h ≥ h ∗ h h h Theorem 4.13 . The observed post-transient radii in Figure 1 agree with this scaling in the original coordinates.
We therefore try beginning with a relatively large learning rate and decreasing it when the average loss stops decreasing. Figure 14 illustrates the initial descent and subsequent oscillating loss at fixed h = 0.1 h = 0.1 h = 0.1
(a) Initial descent
(b) Final 20,000 states
(c) Initial loss
Figure 14. With the learning rate fixed at h = 0.1 h = 0.1 h = 0.1 x 0 = ( 400 , − 300 , 200 ) x_0 = (400, -300, 200) x 0 = ( 400 , − 300 , 200 ) v 0 = 0 v_0 = 0 v 0 = 0 x x x k = 7700 k = 7700 k = 7700
To decide when to reduce h h h
L ‾ k = 1 300 ∑ j = 0 299 L ( x k − j ) , ∣ L ‾ k − L ‾ k − 300 ∣ ≤ 0.10 max { L ‾ k , L ‾ k − 300 } . \overline{L}_k = \frac{1}{300}\sum_{j=0}^{299} L(x_{k-j}), \qquad\left|\overline{L}_k-\overline{L}_{k-300}\right| \le0.10\max\{\overline{L}_k,\overline{L}_{k-300}\}. L k = 300 1 j = 0 ∑ 299 L ( x k − j ) , L k − L k − 300 ≤ 0.10 max { L k , L k − 300 } . After each update, we shift both averaging windows of 300 iterates forward by one iterate. When the inequality holds for 300 consecutive updates, we divide h h h ( x k , v k ) (x_k,v_k) ( x k , v k ) h = 0.1 h = 0.1 h = 0.1 h = 10 − 5 h = 10^{-5} h = 1 0 − 5
For the parameters and initial state of Figure 14 , the reductions occur at k = 5474 , 6764 , 8048 , 9326 k = 5474, 6764, 8048, 9326 k = 5474 , 6764 , 8048 , 9326 h = 10 − 5 h = 10^{-5} h = 1 0 − 5 Figure 15 , the average loss decreases after each reduction, although the loss continues to oscillate.
(a) Loss-triggered schedule
(b) Loss
Figure 15. Adaptive learning-rate reduction. (a) Learning rate versus iterate, where dashed vertical lines mark reductions in learning rate h h h k ≥ 299 k \ge299 k ≥ 299
Autonomous Adam
Adam was introduced by Kingma and Ba [30 ] . We write their step size as h h h m 0 = v 0 = 0 m_0=v_0=0 m 0 = v 0 = 0 x 0 x_0 x 0 k ≥ 1 k\ge1 k ≥ 1
g k = ∇ f k ( x k − 1 ) , m k = β 1 m k − 1 + ( 1 − β 1 ) g k , v k = β 2 v k − 1 + ( 1 − β 2 ) g k ⊙ 2 , m ^ k = m k 1 − β 1 k , v ^ k = v k 1 − β 2 k , x k = x k − 1 − h m ^ k ⊙ ( v ^ k + ϵ 1 ) . \begin{aligned}
g_k&=\nabla f_k(x_{k-1}),\qquad m_k=\beta_1m_{k-1}+(1-\beta_1)g_k,\qquad v_k=\beta_2v_{k-1}+(1-\beta_2)g_k^{\odot2},\\
\widehat{m}_k&=\frac{m_k}{1-\beta_1^k},\qquad\widehat{v}_k=\frac{v_k}{1-\beta_2^k},\qquad x_k=x_{k-1}-h\widehat{m}_k\odot\left(\sqrt{\widehat{v}_k}+\epsilon\mathbf{1}\right).
\end{aligned} g k m k = ∇ f k ( x k − 1 ) , m k = β 1 m k − 1 + ( 1 − β 1 ) g k , v k = β 2 v k − 1 + ( 1 − β 2 ) g k ⊙ 2 , = 1 − β 1 k m k , v k = 1 − β 2 k v k , x k = x k − 1 − h m k ⊙ ( v k + ϵ 1 ) . The denominator convention is epsilon outside, as in Kingma–Ba Algorithm 1 and the cited PyTorch and Keras implementations [30 , 40 , 28 ] . Here f k f_k f k k k k f k = L f_k=L f k = L k k k L L L (3) , so that g k = A x k − 1 g_k=Ax_{k-1} g k = A x k − 1 ( 1 − β 1 k ) − 1 (1-\beta_1^k)^{-1} ( 1 − β 1 k ) − 1 ( 1 − β 2 k ) − 1 (1-\beta_2^k)^{-1} ( 1 − β 2 k ) − 1 k k k
For 0 ≤ β 1 , β 2 < 1 0\le\beta_1,\beta_2<1 0 ≤ β 1 , β 2 < 1 F V F_V F V (4) evaluated at β = β 2 \beta=\beta_2 β = β 2
m k + 1 = β 1 m k + ( 1 − β 1 ) A x k , v k + 1 = F V ( x k , v k ) , x k + 1 = x k − h m k + 1 ⊙ ( v k + 1 + ϵ 1 ) . (58) m_{k+1}=\beta_1m_k+(1-\beta_1)Ax_k,\qquad v_{k+1}=F_V(x_k,v_k),\qquad x_{k+1}=x_k-hm_{k+1}\odot\left(\sqrt{v_{k+1}}+\epsilon\mathbf{1}\right). \tag*{(58)} m k + 1 = β 1 m k + ( 1 − β 1 ) A x k , v k + 1 = F V ( x k , v k ) , x k + 1 = x k − h m k + 1 ⊙ ( v k + 1 + ϵ 1 ) . (58) Using X , W X,W X , W (11) and scaling the momentum as M = m / h M=m/h M = m / h G β 1 G_{\beta_1} G β 1 G 2 G_2 G 2 (12) with β = β 2 \beta=\beta_2 β = β 2
M k + 1 = β 1 M k + ( 1 − β 1 ) A X k , W k + 1 = G 2 ( X k , W k ) , X k + 1 = X k − M k + 1 ⊙ ( W k + 1 + c 1 ) . (59) \begin{aligned}
M_{k+1}&=\beta_1M_k+(1-\beta_1)AX_k,\qquad W_{k+1}=G_2(X_k,W_k),\\
X_{k+1}&=X_k-M_{k+1}\odot(W_{k+1}+c\mathbf{1}). \tag*{(59)}
\end{aligned} M k + 1 X k + 1 = β 1 M k + ( 1 − β 1 ) A X k , W k + 1 = G 2 ( X k , W k ) , = X k − M k + 1 ⊙ ( W k + 1 + c 1 ) . (59) Continuation of invariant circles into Adam
At β 1 = 0 \beta_1=0 β 1 = 0 M M M (59) gives (12) with β = β 2 \beta=\beta_2 β = β 2 β 1 \beta_1 β 1 K j = ( X j , W j ) K_j=(X_j,W_j) K j = ( X j , W j )
The momentum relation in (59) at β 1 = 0 \beta_1=0 β 1 = 0 (54) , the predecessor of K j ( θ ) K_j(\theta) K j ( θ ) K 3 − j ( θ − ω ) K_{3-j}(\theta-\omega) K 3 − j ( θ − ω )
K ~ 1 ( θ ) = ( X 1 ( θ ) , A X 2 ( θ − ω ) , W 1 ( θ ) ) , K ~ 2 ( θ ) = ( X 2 ( θ ) , A X 1 ( θ − ω ) , W 2 ( θ ) ) . \widetilde{K}_1(\theta)=\left(X_1(\theta),AX_2(\theta-\omega),W_1(\theta)\right),\qquad\widetilde{K}_2(\theta)=\left(X_2(\theta),AX_1(\theta-\omega),W_2(\theta)\right). K 1 ( θ ) = ( X 1 ( θ ) , A X 2 ( θ − ω ) , W 1 ( θ ) ) , K 2 ( θ ) = ( X 2 ( θ ) , A X 1 ( θ − ω ) , W 2 ( θ ) ) . We write G 0 G_0 G 0 G β 1 G_{\beta_1} G β 1 β 1 = 0 \beta_1=0 β 1 = 0
G 0 ( K ~ j ( θ ) ) = ( X 3 − j ( θ + ω ) , A X j ( θ ) , W 3 − j ( θ + ω ) ) = K ~ 3 − j ( θ + ω ) . G_0\left(\widetilde{K}_j(\theta)\right)=\left(X_{3-j}(\theta+\omega),AX_j(\theta),W_{3-j}(\theta+\omega)\right)=\widetilde{K}_{3-j}(\theta+\omega). G 0 ( K j ( θ ) ) = ( X 3 − j ( θ + ω ) , A X j ( θ ) , W 3 − j ( θ + ω ) ) = K 3 − j ( θ + ω ) . Using these parameterizations as initial approximations at h = 7.80889203 × 10 − 7 h=7.80889203\times10^{-7} h = 7.80889203 × 1 0 − 7 β 1 = 0.01 , 0.02 \beta_1=0.01,0.02 β 1 = 0.01 , 0.02 Figure 16 compares these continuations with the RMSProp case β 1 = 0 \beta_1=0 β 1 = 0
Figure 16. An RMSProp period-two circle at h = 7.80889203 × 10 − 7 h = 7.80889203 \times10^{-7} h = 7.80889203 × 1 0 − 7 β 1 = 0.01 , 0.02 \beta_1 = 0.01, 0.02 β 1 = 0.01 , 0.02 β 2 = β \beta_2 = \beta β 2 = β
We then follow the period-two circle as β 1 \beta_1 β 1 Figure 17 shows it shrinking toward the period-two orbit near β 1 ≈ 0.21550 \beta_1\approx0.21550 β 1 ≈ 0.21550
Figure 17. Autonomous-Adam circle continuation at h = 7.8089 × 10 − 7 h = 7.8089 \times10^{-7} h = 7.8089 × 1 0 − 7 X 1 X_1 X 1 X 3 X_3 X 3 M 1 M_1 M 1 M 3 M_3 M 3 β 1 \beta_1 β 1 β 1 = 0.21550 \beta_1 = 0.21550 β 1 = 0.21550
At the larger learning rate h = 1.92 × 10 − 6 h=1.92\times10^{-6} h = 1.92 × 1 0 − 6 β 1 = 0 , 0.01 , 0.02 \beta_1=0,0.01,0.02 β 1 = 0 , 0.01 , 0.02 Figure 18 shows their geometry.
Figure 18. Autonomous-Adam trajectories at h = 1.92 × 10 − 6 h = 1.92 \times10^{-6} h = 1.92 × 1 0 − 6 β 1 = 0 , 0.01 , 0.02 \beta_1 = 0, 0.01, 0.02 β 1 = 0 , 0.01 , 0.02 § 19 .
We vary h h h β 1 \beta_1 β 1 195 × 181 195\times181 195 × 181 § 6.3.2 and the parameter ranges in Section C.4.
We classify 6,868 samples as approaching the origin, 2,030 as period two, 7,768 as invariant circles or tori, and 18,629 as chaotic (Figure 19 ). These are finite-window classifications, subject to the limitations in § 6.3.2 . The dashed curve gives the stability boundary of the origin for epsilon-outside Adam without bias correction [3 Corollary III.5] ,
Figure 19. Autonomous Adam in the ( h , β 1 ) (h,\beta_1) ( h , β 1 ) Figure 16 . The dashed curve is the stability threshold of the origin. Panels A–F show X 1 X_1 X 1 X 2 X_2 X 2 ( h / 10 − 6 , β 1 ) = ( 1.05 , 0 ) , ( 1.92 , 0 ) , ( 1.92 , 0.02 ) , ( 3.00 , 0.30 ) , ( 2.40 , 0.75 ) , ( 3.00 , 0.40 ) (h/10^{-6},\beta_1) = (1.05,0),(1.92,0),(1.92,0.02),(3.00,0.30),(2.40,0.75),(3.00,0.40) ( h /1 0 − 6 , β 1 ) = ( 1.05 , 0 ) , ( 1.92 , 0 ) , ( 1.92 , 0.02 ) , ( 3.00 , 0.30 ) , ( 2.40 , 0.75 ) , ( 3.00 , 0.40 )
h = 2 ϵ ( 1 + β 1 ) ( 1 − β 1 ) λ max ( A ) . h=\frac{2\epsilon(1+\beta_1)}{(1-\beta_1)\lambda_{\max}(A)}. h = ( 1 − β 1 ) λ m a x ( A ) 2 ϵ ( 1 + β 1 ) . Appendix G provides additional trajectories, circle projections, and the trajectory settings for panels E and F of Figure 19 .
Conclusion
We begin with the question of what happens to RMSProp when the learning rate is too large for convergence to the minimizer. For positive-definite quadratic losses, we show that the origin is Lyapunov unstable whenever h > h ∗ = 2 ϵ / λ max ( A ) h > h_* = 2\epsilon/\lambda_{\max}(A) h > h ∗ = 2 ϵ / λ m a x ( A ) h h h h 2 h^2 h 2
The next question is what actually happens inside this region. We establish that period-two orbits exist precisely above the critical learning rate. Under additional spectral assumptions, we obtain an attracting branch emerging from the origin near h ∗ h_* h ∗ h h h
Our numerical experiments show that period-two motion is only a small part of the picture. As we increase the learning rate, we find invariant circles, invariant tori, coexistence of attracting sets, and circle-doubling sequences. We also find trajectories with positive finite-time Lyapunov exponents, consistent with chaotic motion. Even for a quadratic loss with a unique minimizer, there is therefore quite a bit of dynamics to understand. Our parameter surveys show that the observed behavior depends on Hessian orientation, memory decay, and dimension. We also continue selected invariant circles into autonomous Adam, beginning from the RMSProp case β 1 = 0 \beta_1 = 0 β 1 = 0
A more systematic study of the parameter space, with broader hyperparameter searches in arbitrary dimensions, could examine how often multistability and chaotic motion occur beyond the selected examples. The work of Mischaikow and collaborators on DSGRN [10 ] suggests one possible starting point, since their framework describes recurrent sets and partitions parameter space according to qualitative dynamics.
Another direction is to study the numerically observed invariant circles and tori, their bifurcations, and the chaotic regimes using rigorous methods. Computer-assisted methods such as those developed by Tucker [45 ] and Wilczak [46 ] provide examples of rigorous proofs of chaotic attractors and suggest approaches to selected examples in RMSProp.
Finally, the results raise questions about the optimization problem that motivates this study. Our learning-rate reduction example shows how decreasing h h h
Acknowledgments
We used OpenAI Codex and Anthropic Claude to assist with mathematical derivations, developing and reviewing numerical code, preparing figures, and editing the manuscript. We also used Google Gemini to review manuscript wording and formatting. The authors assume responsibility for all content.
Fourier parameterization and Newton's method
For a period-two circle, we solve the coupled equations (54) for G G G P = G 2 P = G^2 P = G 2 θ ℓ = 2 π ℓ / M \theta_\ell= 2\pi\ell/M θ ℓ = 2 π ℓ / M ℓ = 0 , … , M − 1 \ell= 0,\ldots,M-1 ℓ = 0 , … , M − 1
K j ( M ) ( θ ) = Re ∑ k ∈ I M K ^ j , k e i k θ , I M = { − ⌊ M / 2 ⌋ , … , ⌈ M / 2 ⌉ − 1 } . K_j^{(M)}(\theta) = \operatorname{Re}\sum_{k\in I_M}\widehat{K}_{j,k}e^{ik\theta}, \qquad I_M = \left\{-\left\lfloor M/2\right\rfloor,\ldots,\left\lceil M/2\right\rceil-1\right\}. K j ( M ) ( θ ) = Re k ∈ I M ∑ K j , k e ik θ , I M = { − ⌊ M /2 ⌋ , … , ⌈ M /2 ⌉ − 1 } . The coefficients are the discrete Fourier transform divided by M M M
Φ ( K 1 ) = ⟨ K 1 ( 0 ) − v r e f , η r e f ⟩ = 0. \Phi(K_1) = \left\langle K_1(0)-v_{\mathrm{ref}},\eta_{\mathrm{ref}}\right\rangle= 0. Φ ( K 1 ) = ⟨ K 1 ( 0 ) − v ref , η ref ⟩ = 0. Here v r e f v_{\mathrm{ref}} v ref η r e f \eta_{\mathrm{ref}} η ref ⟨ K 1 ′ ( 0 ) , η r e f ⟩ ≠ 0 \left\langle K_1'(0),\eta_{\mathrm{ref}}\right\rangle\ne0 ⟨ K 1 ′ ( 0 ) , η ref ⟩ = 0
We initialize Newton’s method from a trajectory or a previously computed parameterization. When starting from a trajectory, we estimate the rotation number and Fourier coefficients using weighted Birkhoff averages, following the references in Section 6.2.
With these initial approximations, we truncate the Fourier series and apply Newton’s method to (55) and the phase condition, solving for ω \omega ω K 1 K_1 K 1 K 2 K_2 K 2 M M M
For a two-torus, we similarly use a two-angle Fourier parameterization K : T 2 → R 6 K:\mathbb{T}^2\to\mathbb{R}^6 K : T 2 → R 6
G 2 ( K ( θ 1 , θ 2 ) ) − K ( θ 1 + Ω 1 , θ 2 + Ω 2 ) = 0 , ( θ 1 , θ 2 ) ∈ T 2 , (60) G^2\left(K(\theta_1,\theta_2)\right)-K(\theta_1+\Omega_1,\theta_2+\Omega_2)=0,\qquad(\theta_1,\theta_2)\in\mathbb{T}^2, \tag*{(60)} G 2 ( K ( θ 1 , θ 2 ) ) − K ( θ 1 + Ω 1 , θ 2 + Ω 2 ) = 0 , ( θ 1 , θ 2 ) ∈ T 2 , (60) with unknown rotation numbers ρ P , j = Ω j / ( 2 π ) ( m o d 1 ) \rho_{P,j}=\Omega_j/(2\pi)\pmod{1} ρ P , j = Ω j / ( 2 π ) ( mod 1 ) j = 1 , 2 j=1,2 j = 1 , 2 P = G 2 P=G^2 P = G 2 Ω j \Omega_j Ω j P P P ω \omega ω G G G
Φ j ( K ) = 1 ( 2 π ) 2 ∫ T 2 ⟨ K ( θ ) − K r e f ( θ ) , ∂ θ j K r e f ( θ ) ⟩ d θ = 0 , j = 1 , 2. (61) \Phi_j(K)=\frac{1}{(2\pi)^2}\int_{\mathbb{T}^2}\left\langle K(\theta)-K_{\mathrm{ref}}(\theta),\partial_{\theta_j}K_{\mathrm{ref}}(\theta)\right\rangle\,d\theta=0,\qquad j=1,2. \tag*{(61)} Φ j ( K ) = ( 2 π ) 2 1 ∫ T 2 ⟨ K ( θ ) − K ref ( θ ) , ∂ θ j K ref ( θ ) ⟩ d θ = 0 , j = 1 , 2. (61) where K r e f K_{\mathrm{ref}} K ref θ = ( θ 1 , θ 2 ) \theta=(\theta_1,\theta_2) θ = ( θ 1 , θ 2 ) ∂ θ j K r e f \partial_{\theta_j}K_{\mathrm{ref}} ∂ θ j K ref j j j 2 × 2 2 \times2 2 × 2 K K K Ω 1 , Ω 2 \Omega_1,\Omega_2 Ω 1 , Ω 2 K K K G G G K ( T 2 ) K(\mathbb{T}^2) K ( T 2 )
Normal Lyapunov exponents
We measure growth perpendicular to a circle or torus by removing the tangent directions at each iteration of P = G 2 P=G^2 P = G 2 [24 ] . Let K K K
z j = K ( θ ( j ) ) , z j + 1 = P ( z j ) . z_j=K(\theta^{(j)}), \qquad z_{j+1}=P(z_j). z j = K ( θ ( j ) ) , z j + 1 = P ( z j ) . Here θ ( 0 ) \theta^{(0)} θ ( 0 ) θ ( j ) \theta^{(j)} θ ( j ) 2 ω 2\omega 2 ω ( Ω 1 , Ω 2 ) (\Omega_1,\Omega_2) ( Ω 1 , Ω 2 ) 2 π 2\pi 2 π z j z_j z j K ′ ( θ ( j ) ) K'(\theta^{(j)}) K ′ ( θ ( j ) ) ∂ θ 1 K ( θ ( j ) ) , ∂ θ 2 K ( θ ( j ) ) \partial_{\theta_1}K(\theta^{(j)}),\partial_{\theta_2}K(\theta^{(j)}) ∂ θ 1 K ( θ ( j ) ) , ∂ θ 2 K ( θ ( j ) ) Q j Q_j Q j R 6 \mathbb{R}^6 R 6 Q j Q_j Q j D P ( z j ) DP(z_j) D P ( z j ) z j + 1 z_{j+1} z j + 1
B j = Q j + 1 T D P ( z j ) Q j . (62) B_j=Q_{j+1}^{\mathsf{T}}DP(z_j)Q_j. \tag*{(62)} B j = Q j + 1 T D P ( z j ) Q j . (62) Definition A.1 (Normal Lyapunov exponents). When the limits exist, we define the normal Lyapunov exponents along the orbit by
λ ⊥ , i ( z 0 ) = lim N → ∞ 1 N log σ i ( B N − 1 ⋯ B 0 ) . \lambda_{\perp,i}(z_0)=\lim_{N\to\infty}\frac{1}{N}\log\sigma_i(B_{N-1}\cdots B_0). λ ⊥ , i ( z 0 ) = N → ∞ lim N 1 log σ i ( B N − 1 ⋯ B 0 ) . Here σ i \sigma_i σ i i i i P P P
We evaluate K K K [1 ] , we use the QR iteration
B j U j = U j + 1 R j , U j T U j = I , ( R j ) i i > 0 , B_jU_j=U_{j+1}R_j,\qquad U_j^{\mathsf{T}}U_j=I,\qquad(R_j)_{ii}>0, B j U j = U j + 1 R j , U j T U j = I , ( R j ) ii > 0 , where R j R_j R j
After discarding initial iterations to reduce dependence on the starting perturbation directions, we estimate the normal Lyapunov exponents over the next N N N P P P
λ ^ ⊥ , i u n w e i g h t e d = 1 N ∑ j = 0 N − 1 log ( R j ) i i , \widehat{\lambda}_{\perp,i}^{\mathrm{unweighted}}=\frac{1}{N}\sum_{j=0}^{N-1}\log(R_j)_{ii}, λ ⊥ , i unweighted = N 1 j = 0 ∑ N − 1 log ( R j ) ii , λ ^ ⊥ , i w e i g h t e d = ∑ j = 0 N − 1 w j log ( R j ) i i . \widehat{\lambda}_{\perp,i}^{\mathrm{weighted}}=\sum_{j=0}^{N-1}w_j\log(R_j)_{ii}. λ ⊥ , i weighted = j = 0 ∑ N − 1 w j log ( R j ) ii . For the weighted estimates, we use the smooth weighting function of Das et al. [12 ] , evaluated at the midpoints
s j = j + 1 / 2 N , w j = exp [ − 1 / ( s j ( 1 − s j ) ) ] ∑ k = 0 N − 1 exp [ − 1 / ( s k ( 1 − s k ) ) ] , j = 0 , … , N − 1. s_j=\frac{j+1/2}{N},\qquad
w_j=\frac{\exp[-1/(s_j(1-s_j))]}{\sum_{k=0}^{N-1}\exp[-1/(s_k(1-s_k))]},\qquad j=0,\ldots,N-1. s j = N j + 1/2 , w j = ∑ k = 0 N − 1 exp [ − 1/ ( s k ( 1 − s k ))] exp [ − 1/ ( s j ( 1 − s j ))] , j = 0 , … , N − 1. Section C.6 gives the computational settings and compares estimates obtained with different iteration counts and initial points on the circles.
Proof of the general trapping theorem
Theorem 4.9 (Trapping by exterior energy decrease; restated). Let F : R d × R ≥ 0 d → R d × R ≥ 0 d F : \mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d} \to\mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d} F : R d × R ≥ 0 d → R d × R ≥ 0 d V : R d × R ≥ 0 d → R V : \mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d} \to\mathbb{R} V : R d × R ≥ 0 d → R V V V
V ( z ) ⟶ + ∞ as ∥ z ∥ ⟶ ∞ , z ∈ R d × R ≥ 0 d . (63) V(z) \longrightarrow+\infty\quad\text{as } \lVert z\rVert\longrightarrow\infty,\qquad z \in\mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d}. \tag*{(63)} V ( z ) ⟶ + ∞ as ∥ z ∥ ⟶ ∞ , z ∈ R d × R ≥ 0 d . (63) Let R ⊂ R d × R ≥ 0 d R \subset\mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d} R ⊂ R d × R ≥ 0 d
V ( F ( z ) ) < V ( z ) for every z ∈ ( R d × R ≥ 0 d ) ∖ R . (64) V(F(z)) < V(z) \qquad\text{for every } z \in\left(\mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d}\right) \setminus R. \tag*{(64)} V ( F ( z )) < V ( z ) for every z ∈ ( R d × R ≥ 0 d ) ∖ R . (64) Then there exists C 0 > 0 C_0 > 0 C 0 > 0 C ≥ C 0 C \geq C_0 C ≥ C 0 Ω C = { z ∈ R d × R ≥ 0 d : V ( z ) ≤ C } \Omega_C = \{z \in\mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d} : V(z) \leq C\} Ω C = { z ∈ R d × R ≥ 0 d : V ( z ) ≤ C }
(i) Compactness. Ω C \Omega_C Ω C
(ii) Containment of the region. R ⊂ int R d × R ≥ 0 d ( Ω C ) R \subset\operatorname{int}_{\mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d}}(\Omega_C) R ⊂ int R d × R ≥ 0 d ( Ω C )
(iii) Strict trapping. F ( Ω C ) ⊂ int R d × R ≥ 0 d ( Ω C ) F(\Omega_C) \subset\operatorname{int}_{\mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d}}(\Omega_C) F ( Ω C ) ⊂ int R d × R ≥ 0 d ( Ω C ) Ω C \Omega_C Ω C
(iv) Finite-time entry. For every z 0 ∈ R d × R ≥ 0 d z_0 \in\mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d} z 0 ∈ R d × R ≥ 0 d N = N ( z 0 , C ) ≥ 0 N = N(z_0,C) \geq0 N = N ( z 0 , C ) ≥ 0 F n ( z 0 ) ∈ Ω C F^n(z_0) \in\Omega_C F n ( z 0 ) ∈ Ω C n ≥ N n \geq N n ≥ N
The interiors in parts (ii) and (iii) are relative to R d × R ≥ 0 d \mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d} R d × R ≥ 0 d
Proof. We refer to V V V R ‾ \overline{R} R R R R R d × R ≥ 0 d \mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d} R d × R ≥ 0 d R R R R ‾ \overline{R} R R d × R ≥ 0 d \mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d} R d × R ≥ 0 d V ∘ F V \circ F V ∘ F F F F V V V R R R
M 1 = max z ∈ R ‾ V ( z ) , M 2 = max z ∈ R ‾ V ( F ( z ) ) . M_1 = \max_{z \in\overline{R}} V(z),\qquad M_2 = \max_{z \in\overline{R}} V(F(z)). M 1 = z ∈ R max V ( z ) , M 2 = z ∈ R max V ( F ( z )) . Choose any C 0 > max { 0 , M 1 , M 2 } C_0 > \max\{0,M_1,M_2\} C 0 > max { 0 , M 1 , M 2 } C 0 > 0 C_0 > 0 C 0 > 0 z ∈ R ‾ z \in\overline{R} z ∈ R
V ( z ) ≤ M 1 < C 0 , V ( F ( z ) ) ≤ M 2 < C 0 . V(z) \leq M_1 < C_0,\qquad V(F(z)) \leq M_2 < C_0. V ( z ) ≤ M 1 < C 0 , V ( F ( z )) ≤ M 2 < C 0 . If R R R C 0 > 0 C_0 > 0 C 0 > 0 C ≥ C 0 C \geq C_0 C ≥ C 0
(i) Compactness. By continuity of V V V Ω C = V − 1 ( ( − ∞ , C ] ) \Omega_C = V^{-1}((-\infty,C]) Ω C = V − 1 (( − ∞ , C ]) R d × R ≥ 0 d \mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d} R d × R ≥ 0 d R 2 d \mathbb{R}^{2d} R 2 d Ω C \Omega_C Ω C R 2 d \mathbb{R}^{2d} R 2 d M > 0 M > 0 M > 0 V ( z ) > C V(z) > C V ( z ) > C z ∈ R d × R ≥ 0 d z \in\mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d} z ∈ R d × R ≥ 0 d ∥ z ∥ > M \lVert z\rVert> M ∥ z ∥ > M z ∈ Ω C z \in\Omega_C z ∈ Ω C ∥ z ∥ ≤ M \lVert z\rVert\leq M ∥ z ∥ ≤ M Ω C \Omega_C Ω C
(ii) Containment of the bounded region. The choice of C 0 C_0 C 0 V ( z ) < C 0 ≤ C V(z) < C_0 \leq C V ( z ) < C 0 ≤ C z ∈ R z \in R z ∈ R R ⊂ { V < C } R \subset\{V < C\} R ⊂ { V < C } V V V { z ∈ R d × R ≥ 0 d : V ( z ) < C } \{z \in\mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d} : V(z) < C\} { z ∈ R d × R ≥ 0 d : V ( z ) < C } R d × R ≥ 0 d \mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d} R d × R ≥ 0 d Ω C \Omega_C Ω C Ω C \Omega_C Ω C
R ⊂ { V < C } ⊂ int R d × R ≥ 0 d ( Ω C ) . R \subset\{V < C\} \subset\operatorname{int}_{\mathbb{R}^{d} \times\mathbb{R}_{\geq0}^{d}}(\Omega_C). R ⊂ { V < C } ⊂ int R d × R ≥ 0 d ( Ω C ) . (iii) Strict trapping. To prove that F F F Ω C \Omega_C Ω C z ∈ Ω C z \in\Omega_C z ∈ Ω C
Case 1: z ∈ R z \in R z ∈ R The choice of C 0 C_0 C 0
V ( F ( z ) ) < C 0 ≤ C . V(F(z)) < C_0 \leq C. V ( F ( z )) < C 0 ≤ C . No comparison between V ( F ( z ) ) V(F(z)) V ( F ( z )) V ( z ) V(z) V ( z ) R R R C C C Case 2: z ∈ Ω C ∖ R z \in\Omega_C \setminus R z ∈ Ω C ∖ R Since z z z R R R z ∈ Ω C z \in\Omega_C z ∈ Ω C V ( z ) ≤ C V(z) \le C V ( z ) ≤ C
Thus
V ( F ( z ) ) < V ( z ) ≤ C . V(F(z)) < V(z) \le C. V ( F ( z )) < V ( z ) ≤ C . In both cases, V ( F ( z ) ) < C V(F(z)) < C V ( F ( z )) < C F ( Ω C ) ⊂ { V < C } ⊂ int R ≥ 0 d × R ≥ 0 d ( Ω C ) F(\Omega_C) \subset\{V < C\} \subset\operatorname{int}_{\mathbb{R}_{\geq0}^d \times\mathbb{R}_{\geq0}^d}(\Omega_C) F ( Ω C ) ⊂ { V < C } ⊂ int R ≥ 0 d × R ≥ 0 d ( Ω C )
(iv) Finite-time entry. Suppose an orbit z n = F n ( z 0 ) z_n = F^n(z_0) z n = F n ( z 0 ) Ω C \Omega_C Ω C K = { z ∈ R ≥ 0 d × R ≥ 0 d : C ≤ V ( z ) ≤ V ( z 0 ) } K = \{z \in\mathbb{R}_{\geq0}^d \times\mathbb{R}_{\geq0}^d : C \le V(z) \le V(z_0)\} K = { z ∈ R ≥ 0 d × R ≥ 0 d : C ≤ V ( z ) ≤ V ( z 0 )} R R R V − V ∘ F V - V \circ F V − V ∘ F K K K δ \delta δ V ( z n ) ≤ V ( z 0 ) − n δ V(z_n) \le V(z_0) - n\delta V ( z n ) ≤ V ( z 0 ) − n δ n n n V ( z n ) < C V(z_n) < C V ( z n ) < C Ω C \Omega_C Ω C
Additional numerical methods
Matrix used in the three-dimensional examples
For the fixed three-dimensional examples, we use the matrix below, with ϵ = 10 − 7 \epsilon= 10^{-7} ϵ = 1 0 − 7 β = 0.98 \beta= 0.98 β = 0.98 β 2 = 0.98 \beta_2 = 0.98 β 2 = 0.98
A = 1 144115188075855872 ( 109836008153513888 − 13893213032827795 12762003043063982 − 13893213032827795 161280689414189920 24193156135327376 12762003043063982 24193156135327376 86288968860418720 ) . A = \frac{1}{144115188075855872}
\begin{pmatrix}
109836008153513888 & -13893213032827795 & 12762003043063982 \\
-13893213032827795 & 161280689414189920 & 24193156135327376 \\
12762003043063982 & 24193156135327376 & 86288968860418720
\end{pmatrix}. A = 144115188075855872 1 109836008153513888 − 13893213032827795 12762003043063982 − 13893213032827795 161280689414189920 24193156135327376 12762003043063982 24193156135327376 86288968860418720 . Continuation and numerical checks
Following a branch. We vary h h h β 1 \beta_1 β 1
Circle checks. We evaluate the invariance equations (54) on angular grids with more points than the grid used by Newton’s method, and offset the evaluation points to check accuracy between the original nodes. We define the invariance error as the ℓ ∞ \ell^\infty ℓ ∞ M M M max ( 2048 , 4 M ) \max(2048, 4M) max ( 2048 , 4 M ) L L L θ ℓ = 2 π ( ℓ + s ) / L \theta_\ell= 2\pi(\ell+ s)/L θ ℓ = 2 π ( ℓ + s ) / L ℓ = 0 , … , L − 1 \ell= 0, \ldots, L - 1 ℓ = 0 , … , L − 1 s = 2 / 10 s = \sqrt{2}/10 s = 2 /10 0.371 0.371 0.371 ( 5 − 1 ) / 2 (\sqrt{5} - 1)/2 ( 5 − 1 ) /2 s s s s s s 2 π / L 2\pi/L 2 π / L max ( 4096 , 4 M ) \max(4096, 4M) max ( 4096 , 4 M )
For the primary period-two circle in Figure 4 , the invariance error is below 5 × 10 − 14 5 \times10^{-14} 5 × 1 0 − 14 h = 1.0514 × 10 − 6 h = 1.0514 \times10^{-6} h = 1.0514 × 1 0 − 6 2 × 10 − 11 2 \times10^{-11} 2 × 1 0 − 11 10 − 10 10^{-10} 1 0 − 10
As a check that we have used enough Fourier modes for each circle component, we require the sum of squared coefficient magnitudes for ∣ k ∣ > 3 4 ⌊ M / 2 ⌋ |k| > \frac{3}{4}\lfloor M/2 \rfloor ∣ k ∣ > 4 3 ⌊ M /2 ⌋ 10 − 8 10^{-8} 1 0 − 8 ⌊ M / 2 ⌋ \lfloor M/2 \rfloor ⌊ M /2 ⌋ ∥ K j ′ ( θ ) ∥ 2 \lVert K_j'(\theta) \rVert_2 ∥ K j ′ ( θ ) ∥ 2 θ \theta θ 10 − 8 10^{-8} 1 0 − 8
Torus checks. For a two-torus, we evaluate the invariance equations (60) at angles between the computational nodes. We also check the magnitudes of the Fourier coefficients near the largest frequencies in each angle, the linear independence of the two tangent vectors, and the positivity of the memory coordinates.
For the two-tori at h = 1.0514 × 10 − 6 h = 1.0514 \times10^{-6} h = 1.0514 × 1 0 − 6 10 − 12 10^{-12} 1 0 − 12 10 − 14 10^{-14} 1 0 − 14 M 1 × M 2 M_1 \times M_2 M 1 × M 2 ( 2 M 1 + 1 ) × ( 2 M 2 + 1 ) (2M_1 + 1) \times(2M_2 + 1) ( 2 M 1 + 1 ) × ( 2 M 2 + 1 ) 3 M 1 × 3 M 2 3M_1 \times3M_2 3 M 1 × 3 M 2 ( 0.371 , 0.237 ) (0.371, 0.237) ( 0.371 , 0.237 ) ( 0.1414 , 0.618 ) (0.1414, 0.618) ( 0.1414 , 0.618 ) ( θ 1 , θ 2 ) (\theta_1,\theta_2) ( θ 1 , θ 2 ) Table 2 gives the computational grids and invariance errors on these denser, offset grids.
Table 2: Fourier grids and errors measured on denser, offset grids at h = 1.0514 × 10 − 6 h = 1.0514 \times10^{-6} h = 1.0514 × 1 0 − 6 (54) , and torus invariance errors to (60) .
Table 2.
Distance to an invariant circle
We measure Euclidean distance in the full normalized phase space. For a state Z = ( X , W ) Z = (X, W) Z = ( X , W ) C = K 1 ( T ) ∪ K 2 ( T ) C = K_1(\mathbb{T}) \cup K_2(\mathbb{T}) C = K 1 ( T ) ∪ K 2 ( T )
d ( Z , C ) = min j ∈ { 1 , 2 } θ ∈ [ 0 , 2 π ) ∥ Z − K j ( θ ) ∥ 2 . (65) d(Z,C) = \min_{\substack{j \in\{1,2\} \\ \theta\in[0,2\pi)}} \lVert Z - K_j(\theta) \rVert_2. \tag*{(65)} d ( Z , C ) = j ∈ { 1 , 2 } θ ∈ [ 0 , 2 π ) min ∥ Z − K j ( θ ) ∥ 2 . (65) We first search a dense angular grid on both circle components, then refine several candidate angles using the Fourier parameterization and retain the smallest distance. To check sensitivity to the angular grid, we repeat the calculation on a second grid.
Survey measurements and longer checks
We use the classification rule in (58) for all three surveys. Table 3 gives the grids, discarded and measured lengths, recurrence window m m m § 6.1 . We compute the leading full-state Lyapunov exponent as described in § 6.3 , reporting it per iterate of the map listed in the final column of Table 3 .
Table 3: Survey grids and measurement settings. Discarded and measured lengths count iterates of the map in the last column. Recurrence uses the final m m m
Table 3.
The first two surveys use the fixed matrix (53) and ϵ = 10 − 7 \epsilon= 10^{-7} ϵ = 1 0 − 7 β 2 = 0.98 \beta_2 = 0.98 β 2 = 0.98 z z z y y y x x x ( 0.5 , 0.8 , 1.18 ) (0.5, 0.8, 1.18) ( 0.5 , 0.8 , 1.18 ) β = 0.98 \beta= 0.98 β = 0.98 ϵ = 10 − 7 \epsilon= 10^{-7} ϵ = 1 0 − 7 h = 3 × 10 − 6 h = 3 \times10^{-6} h = 3 × 1 0 − 6
We selected the tolerances by comparing the classifications produced by candidate rules with those obtained using more post-transient iterations. For samples outside the four categories, we compute longer trajectories and search for higher-period orbits by comparing states separated by more than two iterates.
Bifurcation diagram. For the upper panel of Figure 3 , before each run we perturb every spatial coordinate by an amount in [ − 5 × 10 − 9 , 5 × 10 − 9 ] [-5 \times10^{-9}, 5 \times10^{-9}] [ − 5 × 1 0 − 9 , 5 × 1 0 − 9 ] W i W_i W i max { W i , 10 − 8 } \max\{W_i, 10^{-8}\} max { W i , 1 0 − 8 } (50) .
Individual-trajectory Lyapunov exponent measurements
Table 4 gives the leading full-state estimates per iterate of G G G G β 1 G_{\beta_1} G β 1
The three-dimensional and coupled RMSProp examples use h = 1.06 × 10 − 6 h = 1.06 \times10^{-6} h = 1.06 × 1 0 − 6 η = 0.10 \eta= 0.10 η = 0.10 § 9.0.2 , with all entries of its initial perturbation vector equal to 1 / 2 d 1/\sqrt{2d} 1/ 2 d ρ = 15 \rho= 15 ρ = 15 β = 0.5 \beta= 0.5 β = 0.5 ϵ = 10 − 7 \epsilon= 10^{-7} ϵ = 1 0 − 7 h = 1.92 × 10 − 6 h = 1.92 \times10^{-6} h = 1.92 × 1 0 − 6
Normal Lyapunov exponent measurements
At h = 1.0514 × 10 − 6 h = 1.0514 \times10^{-6} h = 1.0514 × 1 0 − 6 C 0 C_0 C 0 § 13 .
We start the circle calculations at K ( 0 ) K(0) K ( 0 ) K ( 0.713 ) K(0.713) K ( 0.713 ) K ( 0 , 0 ) K(0,0) K ( 0 , 0 ) Table 5 , checking selected torus values against direct Fourier evaluation.
Table 5. Settings for the normal Lyapunov-exponent comparisons at h = 1.0514 × 10 − 6 h = 1.0514 \times10^{-6} h = 1.0514 × 1 0 − 6 P P P
We count positive and negative weighted normal Lyapunov exponent estimates using a threshold of 10 − 7 10^{-7} 1 0 − 7 P P P C 0 C_0 C 0
When we increase the number of measured iterations or change the initial point on a circle, we find that the weighted estimates retain their signs, while some unweighted estimates change sign (Table 6 ). With the larger number of measured iterations, the attracting and saddle tori have leading weighted normal Lyapunov exponents − 2.289 × 10 − 6 -2.289 \times10^{-6} − 2.289 × 1 0 − 6 + 1.178 × 10 − 6 +1.178 \times10^{-6} + 1.178 × 1 0 − 6 P P P
Table 6. Ranges of normal Lyapunov exponent estimates at h = 1.0514 × 10 − 6 h = 1.0514 \times10^{-6} h = 1.0514 × 1 0 − 6 10 − 6 10^{-6} 1 0 − 6 P P P C 0 C_0 C 0
Higher-period motion
We use longer recurrence checks to find periodic motion missed by the fixed-point and period-two tests. For a proposed period p p p z 0 , … , z p − 1 z_0,\ldots,z_{p-1} z 0 , … , z p − 1
R p ( z 0 , … , z p − 1 ) = ( G ( z j ) − z j + 1 ) j = 0 p − 1 , z p = z 0 , (66) \mathcal{R}_p(z_0,\ldots,z_{p-1})=\left(G(z_j)-z_{j+1}\right)_{j=0}^{p-1},\qquad z_p=z_0, \tag*{(66)} R p ( z 0 , … , z p − 1 ) = ( G ( z j ) − z j + 1 ) j = 0 p − 1 , z p = z 0 , (66) R p = 0. \mathcal{R}_p=0. R p = 0. For Adam, we replace G G G G β 1 G_{\beta_1} G β 1 ℓ ∞ \ell^\infty ℓ ∞
∥ R p ∥ ℓ ∞ = max 0 ≤ j < p max i ∣ ( G ( z j ) − z j + 1 ) i ∣ . \|\mathcal{R}_p\|_{\ell^\infty}=\max_{0\leq j<p}\max_i\left|(G(z_j)-z_{j+1})_i\right|. ∥ R p ∥ ℓ ∞ = 0 ≤ j < p max i max ∣ ( G ( z j ) − z j + 1 ) i ∣ . We also check proper divisors of p p p Table 7 lists the six refined orbits, of which the first three are the magenta samples in Figure 11 . We find the other three through longer recurrence searches on samples flagged by the survey criteria.
Table 7. Numerically refined higher-period orbits. Periods count iterates of G G G G β 1 G_{\beta_1} G β 1 ∥ R p ∥ ℓ ∞ \|\mathcal{R}_p\|_{\ell^\infty} ∥ R p ∥ ℓ ∞ (66) .
Additional bifurcation observations
Table 8 records the sequences and invariant objects underlying the examples in § 7 . Here A A A (53) , and, unless stated otherwise, we report Lyapunov exponents per iteration of P = G 2 P=G^2 P = G 2 P P P
Table 8. Observed invariant objects and transition diagnostics for A A A 10 − 6 10^{-6} 1 0 − 6
Table 8 (continued).
Additional doubling observations in the CNS family
In Figure 20 , we compare four computed circles in the C N S C_{\mathrm{NS}} C NS q q q
Figure 20. Additional doubling observations in the C N S C_{\mathrm{NS}} C NS h = 1.001845 × 10 − 6 h = 1.001845 \times10^{-6} h = 1.001845 × 1 0 − 6 q = 4 , 8 , 16 , 32 q = 4, 8, 16, 32 q = 4 , 8 , 16 , 32 X 1 X_1 X 1 X 2 X_2 X 2
We compute distinct 16-fold and 32-fold circles at h = 1.001845 × 10 − 6 h = 1.001845 \times10^{-6} h = 1.001845 × 1 0 − 6 h = 1.0018475 × 10 − 6 h = 1.0018475 \times10^{-6} h = 1.0018475 × 1 0 − 6 Table 9 . At the latter learning rate, we find numerical evidence of a non-orientable leading normal line bundle. A unit vector chosen continuously along this bundle returns with the opposite sign after one circuit of each component. Together, these observations suggest another circle doubling.
Table 9. Normal Lyapunov exponent estimates and final distances for the 32-fold circle. Estimates are per iteration of P P P P P P
We also check whether nearby trajectories approach the 32-fold circle on either side of the sign change. At each learning rate, we follow 16 states displaced by 10 − 7 10^{-7} 1 0 − 7 P P P (65) and the angular refinement in Section C.3. The distances agree closely on the two evaluation grids. We compare the final distances at the two learning rates in Table 9 .
Additional numerical figures
Iterates alternate between the two components of each circle or torus pair, and P = G 2 P = G^2 P = G 2
Figure 23. At large h h h X 1 X_1 X 1 X 2 X_2 X 2 h ∈ { 2.4 × 10 − 6 , 6.654 × 10 − 6 , 10 − 5 , 10 − 4 , 10 − 2 , 1 } h \in\{2.4 \times10^{-6}, 6.654 \times10^{-6}, 10^{-5}, 10^{-4}, 10^{-2}, 1\} h ∈ { 2.4 × 1 0 − 6 , 6.654 × 1 0 − 6 , 1 0 − 5 , 1 0 − 4 , 1 0 − 2 , 1 } G G G
Epsilon-inside comparison
For the epsilon-inside map G I G_I G I h ∗ , I = 2 ϵ / λ max ( A ) h_{*,I}=2\sqrt{\epsilon}/\lambda_{\max}(A) h ∗ , I = 2 ϵ / λ m a x ( A ) ρ I = h / h ∗ , I \rho_I=h/h_{*,I} ρ I = h / h ∗ , I A A A ρ I = 2 \rho_I=2 ρ I = 2 ρ I = 4 \rho_I=4 ρ I = 4
Figure 24. Bifurcation diagram for the normalized epsilon-inside map G I G_I G I (49) . Panel (a) uses a linear axis for h h h h ∗ , I = 2 ϵ / λ max ( A ) h_{*,I}=2\sqrt{\epsilon}/\lambda_{\max}(A) h ∗ , I = 2 ϵ / λ m a x ( A ) ρ I = h / h ∗ , I \rho_I=h/h_{*,I} ρ I = h / h ∗ , I 10 − 8 10^{-8} 1 0 − 8 G I G_I G I X 1 X_1 X 1
Additional RMSProp trajectories
In the six examples below, we plot bounded trajectories with positive leading finite-time Lyapunov exponents after discarding the initial transient.
Figure 25. RMSProp chaotic-dynamics gallery in the X 1 X_1 X 1 X 2 X_2 X 2 A A A β = 0.98 \beta= 0.98 β = 0.98 ϵ = 10 − 7 \epsilon= 10^{-7} ϵ = 1 0 − 7 X = ( 0.3 , 0.4 , − 0.3 ) X = (0.3, 0.4, -0.3) X = ( 0.3 , 0.4 , − 0.3 ) W = ( 0.2 , 0.2 , 0.2 ) W = (0.2, 0.2, 0.2) W = ( 0.2 , 0.2 , 0.2 ) G G G G G G G G G
Additional autonomous-Adam plots
Additional circle continuations
We continue autonomous-Adam period-two circles in β 1 \beta_1 β 1 h = 6.30261051 × 10 − 7 h = 6.30261051 \times10^{-7} h = 6.30261051 × 1 0 − 7 7.80889203 × 10 − 7 7.80889203 \times10^{-7} 7.80889203 × 1 0 − 7 Table 10 lists the last pair computed on each branch. For one component, we define the spatial amplitude as the largest coordinate range along the circle, max i ( max θ X i ( θ ) − min θ X i ( θ ) ) \max_i(\max_\theta X_i(\theta) - \min_\theta X_i(\theta)) max i ( max θ X i ( θ ) − min θ X i ( θ )) G β 1 2 G_{\beta_1}^2 G β 1 2
Table 10. Last computed period-two circles on two autonomous-Adam continuation branches. The invariance errors are measured for the second iterate G β 1 2 G_{\beta_1}^2 G β 1 2
At the smaller learning rate, the circles approach the period-two orbit near β 1 ≈ 0.085913 \beta_1 \approx0.085913 β 1 ≈ 0.085913 β 1 = 0.08591318628 \beta_1 = 0.08591318628 β 1 = 0.08591318628 h = 6.30261051 × 10 − 7 h = 6.30261051 \times10^{-7} h = 6.30261051 × 1 0 − 7 β 1 = 0.21550405864 \beta_1 = 0.21550405864 β 1 = 0.21550405864 h = 7.80889203 × 10 − 7 h = 7.80889203 \times10^{-7} h = 7.80889203 × 1 0 − 7 Figure 26 gives further spatial and momentum projections.
Figure 26. Additional X X X M M M h = 6.3026 × 10 − 7 h = 6.3026 \times10^{-7} h = 6.3026 × 1 0 − 7 7.8089 × 10 − 7 7.8089 \times10^{-7} 7.8089 × 1 0 − 7 β 2 = 0.98 \beta_2 = 0.98 β 2 = 0.98 ϵ = 10 − 7 \epsilon= 10^{-7} ϵ = 1 0 − 7 β 1 \beta_1 β 1 β 1 = 0.085913 \beta_1 = 0.085913 β 1 = 0.085913 0.21550 0.21550 0.21550
Trajectory examples
The fixed-learning-rate trajectory figures below, Figure 18 , and panels E and F of Figure 19 use A A A β 2 = 0.98 \beta_{2}=0.98 β 2 = 0.98 ϵ = 10 − 7 \epsilon=10^{-7} ϵ = 1 0 − 7 ( X , M , W ) (X,M,W) ( X , M , W ) § 6.1 . Table 11 gives the trajectory lengths.
Table 11: Trajectory lengths for the fixed-learning-rate Adam trajectory figures. All lengths count iterates of G β 1 G_{\beta_{1}} G β 1
Table 11.
For both trajectory lengths in Table 11 , we plot two consecutive states out of every four from the final 50,000 states, retaining both even and odd iterates. The finite-time Lyapunov exponent estimates use different numbers of discarded and measured iterations, given in Table 4 for the individual examples and in the caption of Figure 33 for the β 1 \beta_{1} β 1
Figure 33. Leading finite-time Lyapunov exponent versus β 1 \beta_1 β 1 h = 1.05 × 10 − 6 h = 1.05 \times10^{-6} h = 1.05 × 1 0 − 6 β 2 = 0.98 \beta_2 = 0.98 β 2 = 0.98 ϵ = 10 − 7 \epsilon= 10^{-7} ϵ = 1 0 − 7 β 1 = 0 , 0.01 , … , 0.80 \beta_1 = 0, 0.01, \ldots, 0.80 β 1 = 0 , 0.01 , … , 0.80 G β 1 G_{\beta_1} G β 1 § 6.1 and the first coordinate vector as the initial perturbation. The indicated stability boundary of the origin is β 1 ≈ 0.72203 \beta_1 \approx0.72203 β 1 ≈ 0.72203
Figure 27. Spatial projections of autonomous-Adam trajectories at h = 1.92 × 10 − 6 h=1.92\times10^{-6} h = 1.92 × 1 0 − 6 β 1 ∈ { 0 , 0.01 , 0.02 , 0.05 , 0.10 , 0.15 , 0.20 , 0.25 , 0.30 } \beta_{1}\in\{0,0.01,0.02,0.05,0.10,0.15,0.20,0.25,0.30\} β 1 ∈ { 0 , 0.01 , 0.02 , 0.05 , 0.10 , 0.15 , 0.20 , 0.25 , 0.30 }
Figure 28. Spatial projections of autonomous-Adam trajectories at h = 1.92 × 10 − 6 h = 1.92 \times10^{-6} h = 1.92 × 1 0 − 6 β 1 ∈ { 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 } \beta_1 \in\{0.4, 0.5, 0.6, 0.7, 0.8, 0.9\} β 1 ∈ { 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 } β 1 = 0.8 \beta_1 = 0.8 β 1 = 0.8 β 1 = 0.9 \beta_1 = 0.9 β 1 = 0.9
Figure 29. Spatial projections of autonomous-Adam trajectories at h = 3 × 10 − 6 h = 3 \times10^{-6} h = 3 × 1 0 − 6 β 1 ∈ { 0 , 0.01 , 0.02 , 0.05 , 0.10 , 0.15 , 0.20 , 0.25 , 0.30 } \beta_1 \in\{0, 0.01, 0.02, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30\} β 1 ∈ { 0 , 0.01 , 0.02 , 0.05 , 0.10 , 0.15 , 0.20 , 0.25 , 0.30 }
Figure 30. Spatial projections of autonomous-Adam trajectories at h = 3 × 10 − 6 h = 3 \times10^{-6} h = 3 × 1 0 − 6 β 1 ∈ { 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 } \beta_1 \in\{0.4, 0.5, 0.6, 0.7, 0.8, 0.9\} β 1 ∈ { 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 } β 1 = 0.9 \beta_1 = 0.9 β 1 = 0.9
Figure 31. Spatial projections of autonomous-Adam trajectories at h = 1.05 × 10 − 6 h = 1.05 \times10^{-6} h = 1.05 × 1 0 − 6 β 1 ∈ { 0 , 0.02 , 0.10 , 0.30 , 0.50 , 0.65 , 0.70 , 0.72 , 0.73 } \beta_1 \in\{0, 0.02, 0.10, 0.30, 0.50, 0.65, 0.70, 0.72, 0.73\} β 1 ∈ { 0 , 0.02 , 0.10 , 0.30 , 0.50 , 0.65 , 0.70 , 0.72 , 0.73 } β 1 = 0.73 \beta_1 = 0.73 β 1 = 0.73 β 1 = 0.72 \beta_1 = 0.72 β 1 = 0.72 β 1 = 0.73 \beta_1 = 0.73 β 1 = 0.73
Figure 32. Spatial projections of autonomous-Adam trajectories near β 1 = 0.275 \beta_1 = 0.275 β 1 = 0.275 h = 1.05 × 10 − 6 h = 1.05 \times10^{-6} h = 1.05 × 1 0 − 6 β 1 ∈ { 0.270 , 0.272 , 0.274 , 0.276 , 0.278 , 0.280 } \beta_1 \in\{0.270, 0.272, 0.274, 0.276, 0.278, 0.280\} β 1 ∈ { 0.270 , 0.272 , 0.274 , 0.276 , 0.278 , 0.280 }
Figure 34. Spatial projections at h = 1.92 × 10 − 6 h = 1.92 \times10^{-6} h = 1.92 × 1 0 − 6 β 1 = 0 \beta_1 = 0 β 1 = 0 0.025 0.025 0.025 0.05 0.05 0.05
