Year: 2021

Abstract:
We theoretically analyze the Feedback Alignment (FA) algorithm, an efficient alternative to backpropagation for training neural networks. We provide convergence guarantees with rates for deep linear networks for both continuous and discrete dynamics. Additionally, we study incremental learning phenomena for shallow linear networks. Interestingly, certain specific initializations imply that negligible components are learned before the principal ones, thus potentially negatively affecting the effectiveness of such a learning algorithm; a phenomenon we classify as implicit anti-regularization. We also provide initialization schemes where the components of the problem are approximately learned by decreasing order of importance, thus providing a form of implicit regularization.

heading: 1 INTRODUCTION
text: In order to train a Machine Learning (ML) architecture in the context of supervised learning, it is a consolidated practice to resort to the Gradient Descent (GD) algorithm and its variants (e.g. stochastic GD) and to calculate the gradient of the loss function via backpropagation (Rumelhart et al., 1986). However, one of its limitations is the so-called weight transport problem (Grossberg, 1987): the update of each neuron during the learning phase relies on the downstream weights (which are therefore "transported" along the backward pass). This implies that 1) the synaptic weights are the same, up to transposition, for both the forward (inference) and backward (learning) pass and 2) each neuron has a precise knowledge of all of its forward synapses. It is known that such a model does not accurately reflect the behaviour of a human brain, thus it is biologically implausible (Crick, 1989). Feedback Alignment (FA), originally proposed by Lillicrap (Lillicrap et al., 2016), has been an attempt to mitigate this discrepancy between the neuroscience side and the engineering side of AI and to tackle the weight transport problem in a systematic way. The simple, yet powerful idea is to replace the transposed weight matrices in the backpropagation algorithm with arbitrary, fixed matrices: via this strategy, each neuron during training receives a random projection of the error vector.
In order to illustrate the algorithm and for the purpose of our study, let us consider a fully-connected L-layer neural network (NN): given the input, x = h0 ∈ Rd, the predicted label ŷ ∈ Ro is computed as a` = W`h`−1 , h` = σ(a`) , ` ∈ {1, . . . , L} := [L], (1) σ being a (potentially) nonlinear function applied entry-wise to the vector a`, and finally ŷ = f(aL) for some output nonlinear function f . We consider a regression task with loss function L being the standard Mean Square Error (MSE). The backpropagation strategy for GD prescribes the following updates for each layer in the network: given a step size η > 0 and the error vector e = δaL = ∂L∂ŷ = ŷ − y where y is the true label, we have
δWL = −ηeh>L−1, δW` = −ηδa`h>`−1 and δa` = ( W>`+1δa`+1 ) σ′(a`), ` ∈ [L−1] (GD)
where denotes the Hadamard product. The FA recipe prescribes the following:
Feedback Alignment: in (GD) use δa` = (M`δa`+1) σ′(a`) , ` ∈ [L− 1] . (2)
where {M`} are a collection of arbitrarily chosen matrices, which do not evolve during training.
Related and previous work. Inspired by Lillicrap et al. (2016) and shortly after its publication, Nøkland (2016) proposed two new algorithms called Direct and Inverse FA. In particular, Direct FA reduces to the FA algorithm in the case of a 2-layer NN or multi-layer linear NN; the difference between these two algorithms in the general multi-layer nonlinear NN resides in injecting the error vector directly into the update of each hidden layer of the network: δa` = (M`e) σ′(a`), ` ∈ [1, L − 1]. Such a modification of the FA algorithm can be interpreted as a noisy version of layer-wise training (Gilmer et al., 2017).
Since its first formulation, there has been an extensive numerical analysis on testing whether FA and Direct FA algorithms can be successfully applied to ML problems, in particular in the presence of complex datasets and deep architectures (Bartunov et al., 2018; Moskovitz et al., 2018; Launay et al., 2019). We mention the recent survey paper Launay et al. (2020) where the focus is on Direct FA and it is shown that the algorithm performs well in modern deep architectures (Transformers and Graph NN, among others) applied to problems like neural view synthesis or natural language processing. Furthermore, more on the applied (biologically-inspired) side, there have been recent pushes in furthering the idea of FA: spike-train level direct feedback alignment (Lee et al., 2020), memory-efficient direct feedback alignment (Chu et al., 2020) and direct random target projection (Frenkel et al., 2021) to mention a few examples.
On the other hand, rigorous theoretical results are scarce. Preliminary asymptotic results can be found in the original papers by Lillicrap et al. (2016) and Nøkland (2016). More recently, Refinetti et al. (2020) presented the derivation of a system of differential equations that governs the Direct FA dynamics for a 2-layer non-linear NN in the regime where the input dimension goes to infinity. There is no explicit formula for the solution of the system, but via a careful local analysis it is possible to identify two separate phases of learning (alignment and memorization) and a "degeneracy breaking" effect: at the end of the FA training, the selected solution maximizes the overlap between the second weight matrix W2 and the FA matrix M . Such an "alignment" dynamics is also evident in our analysis (Section 2). Convergence guarantees for deep networks are still missing from the literature and one of the goals of this paper is indeed to provide rigorous proofs of convergence, together with convergence rates.
Our analysis is in a similar vein as some recent works on implicit regularization in linear neural networks Saxe et al. (2018); Gidel et al. (2019); Arora et al. (2018). Even though our setting is very similar, our analysis is different: because of the feedback alignment matrices, the autonomous differential equation obtained significantly differ and lead to different dynamics. Moreover our focus is slightly different from Saxe et al. (2018); Gidel et al. (2019); Arora et al. (2018) which are fully dedicated to implicit regularization of GD while we also aim at proving the convergence of FA.
Our contributions. We aim to provide a theoretical understanding of how and when FA works. We focus on two aspects: 1) proving convergence of the algorithm and 2) analyzing implicit regularization phenomena. We will extensively study the case of deep linear NNs. Despite being sometimes regarded as too simplistic, the study of linear networks is a crucial starting point for a systematic analysis of the FA algorithm: it will allow us to rigorously analyze a tractable model and it will give us insights and intuition on the general nonlinear setting (Section 7).
Our main results are the following.
• We prove convergence of the FA continuous dynamics, with rates, for L-layer NNs with L ≥ 2 (Sections 2 and 3). Such results hold for any input dimension and number of neurons: the network is not necessarily overparametrized.
• We prove that, for certain initialization schemes, the continuous FA dynamics sequentially learn the solutions of a reduced-rank regression problem, but in a backward fashion: smaller components first, dominant components later (Section 4). Additionally, we provide guidelines for an initialization scheme that avoids such phenomena.
• Finally (Section 5), we analyze the corresponding discrete FA dynamics and we prove convergence to the true labels with linear convergence rates.
We assumed some conditions on the data and on the model in order to render the exposition and the reading more clear. However, some of the "strong" assumptions that will appear in the paper can be easily relaxed (see Section 2).

heading: 2 WARM-UP: SHALLOW NETWORKS
text: In this section, we consider a 2-layer linear NN with vector output ŷ = W2W1x ∈ Ro and input vector x ∈ Rd such that (x,y) ∼ D for some distribution D. Define Σxx := E[xx>] ∈ Rd×d, that we assume to be positive definite, and Σxy := E[yx>] ∈ Ro×d, then the following result holds: Proposition 1. For any distribution D, the FA dynamics that dictates the learning process is
Ẇ1 = M (Σxy −W2W1Σxx) , Ẇ2 = (Σxy −W2W1Σxx)W>1 , (3)
where the dot refers to the time derivative.
Note that this result particularly holds for empirical distributions corresponding to finite datasets.
In this work, we will assume that the two matrices Σxx, Σxy can be simultaneously decomposed as Σxx = V ΛxxV > and Σxy = UΛxyV >, where the latter corresponds to the singular value decomposition (SVD) of Σxy. This assumption holds in many situations, as discussed in details by Gidel et al. (2019). Moreover, a perturbation analysis similar as the one done in Gidel et al. (2019, Theorem 1) could be performed in the general case in order to handle the non-commutative case. Next, we introduce the change of variablesW1 = RW̃1V > andW2 = UW̃2R>, for some arbitrary left-orthogonal matrixR, and we choose the FA matrixM to be decomposable in a similar fashion: M = RDU>, for some matrixD with positive diagonal entriesDii > 0 and zero entries otherwise. In particular, this implies that M is full rank, therefore guaranteeing convergence to a global minimum of the dynamics. The equations for the FA flow will then become
˙̃W1 = D ( Λxy − W̃2W̃1Λxx ) ˙̃W2 = ( Λxy − W̃2W̃1Λxx ) W̃>1 . (4)
Moreover, by the change of variableW ′1 := W̃1Λ −1/2 xx ,W ′2 := Λ −1/2 xx W̃2,D′ := DΛ −1/2 xx where Λ 1/2 xx is the coordinate-wise square root of Λxx,1 we get after some simple calculations that
Ẇ ′1 = D ′ ( Λ1/2xx ΛxyΛ −1/2 xx −W ′2W ′1 ) , Ẇ ′2 = ( Λ1/2xx ΛxyΛ −1/2 xx −W ′2W ′1 ) W ′1 > . (5)
Therefore, by considering Λ′xy := Λ 1/2 xx ΛxyΛ −1/2 xx , we may equivalently assume Σxx = Id (isotropic features), without any loss of generality.
It is now crucial to observe that ifW ′1 ∈ Rh×d,W ′2 ∈ Ro×h are diagonal matrices, then the dynamics decouples into k independent equations with k = min{d, h, o}. In particular, let θi1 = (W ′1)ii and θi2 = (W ′ 2)ii and assume that at time t = 0 the matricesW ′ 1(0) andW ′ 2(0) are diagonal. Then, for each i = 1, . . . , k we obtain the following scalar system:
θ̇i1 = d i ( λi − θi2θi1 ) θ̇i2 = ( λi − θi2θi1 ) θi1, (6)
where λi = (Λxy)ii and di = Dii. Note that for i > k one could define diagonal coefficients for the matricesW ′1 orW ′ 2, but their derivative will be 0, thus non-trivial dynamics only occur for i ∈ [k].

heading: 2.1 THE SCALAR DYNAMICS
text: We now turn our attention to analyzing a nonlinear system of the form (6) for the scalar functions θ1, θ2 ∈ C1(R) (we suppressed the indices). Theorem 2. For any FA constant d ∈ R>0, and ∀ θ0 ∈ R initial value with initialization scheme
θ1(0) = θ0 and θ2(0) = θ0 2
2d , (7)
∃ C > 0 such that the product θ2(t)θ1(t) of the solution of the system (6) converges exponentially to the signal λ in a monotonic fashion: |θ2(t)θ1(t)− λ| < Ce− 3 2 (dλ) 2 3 t as t→ +∞. Moreover the
first layer gets "aligned" with the feedback matrix such that θ1(t)→ 3 √ 2rd.
1Note that we used a slight abuse of notation: depending on the matrix it is multiplied to Λ−1/2xx corresponds to the diagonal matrix either of size d×d or o×o with the diagonal coefficients [Λ−1/2xx ]ii , 1 ≤ i ≤ rank(Λxx).
The proof can be found in Appendix A.2: along with the proof, it appears that the FA weights drive the learning of the forward weights for the first layers to the point that their limit is "aligned" to the value of the FA constants (in a nonlinear fashion). Similarly, this holds for the deep network case (Section 3 and Appendix B.2). One can also observe that components with a large FA constant d are learned faster (see (35)–(36) in Appendix A.2), boosting the alignment.
From Theorem 2 it is clear that the choice of the FA constant will force the optimization to select one particular pair of solutions (θ1(t), θ2(t)) among the infinitely many possible solutions (the equation θ1θ2 = λ is a 1-dimensional manifold on the R2-plane). However, all of these solutions represent global minima (those are the only stationary points of FA), and therefore via the FA dynamics, we can find the global minima of the loss function.
We can easily notice that zero-initialization does guarantee convergence for FA, unlike in the case of GD dynamics (Figure 1b). Furthermore, we want to highlight the ideological novelty of the proposed initialization scheme (7): unlike considering all weight matrices as independent players at the beginning of the optimization process and initializing them at random, we are imposing that once the weight matrixW1 of the first layers is initialized, the initial value of the weight matrixW2 of the second layer is deterministic, and it depends directly from the first layer.
More generally, it is easy (though lengthy) to prove that the product of the solutions (θ1(t), θ2(t)) of the system (6) always converges to the signal λ, regardless of the initialization scheme considered. Some intuitions can be built from the phase diagram in Fig 1a where one can see that the trajectories of the dynamics lie on parabolas with equation 2dθ2 = θ21 +2dK while the solution set is represented by the hyperbola θ1θ2 = λ. The difficulties in the analysis arise when the initialization correspond to a parabola that intersect the solution set more than once. Moreover, the non-mononicity of the quantity |θ1(t)θ2(t)− λ| can be seen on this diagram as the result of some trajectories getting close to the bottom left part of the hyperbola but eventually converging to the top right one.
Theorem 3. For any FA constant d ∈ R>0, and for any initial values θ1(0) and θ2(0), the system (6) admits a unique solution (θ1(t), θ2(t)) such that θ1(t), θ2(t) are bounded for all times t ≥ 0 and θ2(t)θ1(t)→ λ as t→ +∞. The convergence rate is exponential for any pair of initial values (θ1(0), θ2(0)) except for a set of (Lebesgue) measure zero in R2.
For the proof, see Appendix A.3. The main consequence of Theorem 3 is that even by initializing all layers randomly and independently (as it is common practice), the FA dynamics still converges to the true signal. However, we want to stress that the initialization scheme (7) is more suitable for ML applications: for arbitrary initial conditions, the product of the solutions θ2(t)θ1(t) may not be monotonic (compare with Theorem 2). Thus it may trigger inefficient early stopping practices (see Figure 1c). Furthermore, for a randomly and independently sampled pair of initial conditions, we may observe a highly undesirable "reversed" incremental learning phenomenon (Section 4), which is avoided if (7) is put into place.
Relaxing some assumptions. Considering the the assumption that the weight matrices W1,W2 share the left/right singular vectors with Σxy is not strictly needed: with more work, the more general case can be analyzed via a perturbation analysis, along the same lines as in Gidel et al. (2019) (Assumption 1). Furthermore, even though the initialization of the matrices W̃1, W̃2 are diagonal
Saxe et al. (2018) argues that the diagonal solutions provide good approximations to the full solutions when W̃1(0), W̃2(0) are full matrices, initialized with small random weights.

heading: 3 DEEP NETWORKS
text: In this section we extend the results of Section 2 to deep linear NNs: given a distribution of input-output pairs (x,y) ∼ D, the predicted output vector is now ŷ = WL . . .W2W1x with weight matrices W` ∈ Rh`×h`−1 (h0 = d, hL = o) and we note Σxx = E[xx>] ∈ Rd×d and Σxy = E[yx>] ∈ Ro×d. Proposition 4. For any distribution D, the FA equations of motion are
Ẇ` = M`(Σxy −W[1:L]Σxx)W>[1:`−1] ` ∈ [L], (8)
for some FA weight matrices {M`}`∈[L], withML = I , and withW[1:`] = W` . . .W1.
Notice that the legitimate FA update would be Ẇ` = M̃[`:L−1](Σxy−W[1:L]Σxx)W>[1:`−1], however since we are requiring the FA matrices to be full rank, we can redefine the FA matrices as M` = M̃[`:L−1]. We now perform a similar change of variables as in Section 2: W` = R`W̃`R>`−1 with R0 = V ,RL = U and {R`}`=1,...,L−1 arbitrary left-orthogonal matrices. We additionally assume that the FA matrices M` have the following prescribed form M` = R`D`U> (with D` a matrix with zero entries except for entries on the diagonal, which are positive). Then, the system becomes
˙̃W` = D`(Λxy − W̃[1:L]Λxx)W̃>[1:`−1] ` ∈ [L]. (9)
Considering the change of variable W ′1 := W̃1Λ −1/2 xx , W ′L := Λ −1/2 xx W̃L, D′` := D`Λ −1/2 xx , Λ′xy := Λ 1/2 xx ΛxyΛ −1/2 xx and assuming as in §2 thatW ′j(0) are diagonal, the system decouples:
θ̇i` = d i `(λ i − θiL . . . θi1)θi`−1 . . . θi1, ` ∈ [L], (10)
and for each entry i = 1, . . . , k, with k = min{d, h1, . . . , hL−1, o}, the analysis reduces to studying a system of scalar functions {θi1, . . . , θiL}. Theorem 5. For any set of FA constants {d`}`∈[L−1], d` ∈ R>0 and for any θ0 ∈ R, there exists an initialization scheme of the form
θ1(0) = θ0, θ`(0) = C`θ 2`−1 0 , ` = 2, . . . , L, (11)
where C` ∈ R>0 depends on {d1, . . . , d`}, such that the system of differential equations (10) has a unique solution (θ1(t), . . . , θL(t)), which converges exponentially to the signal λ as t→ +∞ in a monotone fashion.
As in the 2-layer case, the proposed initialization scheme is characterized by the fact that only the first layer is initialized at random, while all the other layers depend deterministically on it. Note that, along the same lines as Theorem 3, it is also possible to prove convergence to the true signal for any initialization set {θ`(0)}`∈[L]. However, some desirable properties of the dynamics (monotonicity, exponential convergence) may be lost. For the proofs and the explicit construction of the initialization scheme (11) see Appendix B.2.

heading: 4 IMPLICIT REGULARIZATION PHENOMENA
text: In this section, we analyze the phenomenon of hierarchical learning (or implicit regularization) for FA, as described in Gidel et al. (2019) and Gissin et al. (2019) in the case of GD: along the optimization path, features of the model are learned in a sequential way, so that one feature is fully learned before the dynamics starts to learn the next feature. For FA, we may observe a similar behavior as GD, however, in certain settings the dynamics seems to favor the learning of secondary component prior than the learning of principal components, i.e. signals λi’s with small magnitude are learned earlier than those of big magnitude. This phenomenon is potentially highly problematic since smaller singular values λi’s are normally tied to noisy measurements and do not represent nor
give meaningful information on the structure of the data. Privileging negligible features goes against the standard principles of reasoned data analysis (Principal Component Analysis, for example). We will see that the occurrence of this behavior depends on particular initializations of the algorithm and it can be therefore successfully avoided.
Theorem 6 (Informal). Let λ1 > . . . > λk be the principal components of Λxy. There exists initialization schemes such that the component λi is incrementally learned at time Ti. Moreover, there exists two different initialization schemes such that
Anti-Regularization(§4.1): T1 > . . . > Tk and Regularization(§4.2): T1 < . . . < Tk (12)
To illustrate such a phenomenon and rigorously describe it, consider the 2-layer case in the SVDsetting as in Section 2: ∀ i = 1, . . . , k, the system (6) can be equivalently rewritten as
θ̇i1 = − 1
2
( (θi1) 3 + 2dKiθi1 − 2dλi ) , θi1(0) = θ i 0, (13)
and θi2(t) = 1 2dθ i 1(t) 2 +Ki, with Ki = θi2(0)− 12di (θ i 1(0)) 2 ∈ R and d, λi > 0 (for simplicity we assumed di = d, ∀ i). Depending on the value of Ki (i.e. the initial condition of the i-th system), the cubic polynomial on the right-hand-side of (13) may have between 1 and 3 distinct real roots. The value of the discriminant ∆ = −4d2 ( 8d(Ki)3 + 27(λi)2 ) will dictate which regime we fall into.
We will be considering the case where ∆ > 0 and the case Ki = 0 (i.e. the initialization scheme (7)), which is a special instance of the case ∆ < 0. The case ∆ = 0 can be equally analyzed, but it will not be considered in this paper, as it has zero probability of occurring.2 The sequential learning phenomenon becomes evident when we perform a double scaling limit of the time variable t 7→ δt and the initial condition θi0 7→ f(δ, θi0), δ → +∞, as it will be clear below. Note that the rescaling is uniform across the components i. Proof of Theorem 7 and lengthier discussion can be found in Appendix C.

heading: 4.1 CASE ∆ > 0: BACKWARD LEARNING
text: If Ki < − 32 (λ i)3/2 < 0, the cubic polynomial in (13) has 3 real, distinct root ri1 < r i 2 < r i 3. The intuition is that by initializing θi0 at a value that is very close to r i 2 and applying the time scaling t 7→ δt, δ 1 (i.e. we are "fast-forwarding" the solution θ1(δt)), we can clearly identify the time threshold T i after which the signal is fully learned. More formally,
Theorem 7. ∀ i = 1, . . . , k, if θi1(0) = ri2 + e−δ (δ > 0), then the solution θi1(δt) converges to a step function:
θi1(δt)→ ri2I{t<T i} + αiI{t=T i} + ri3I{t>T i} as δ → +∞, (14)
where I{t∈A} is the indicator function on the set A,
T i = 2
(ri3 − ri2)(ri2 − ri1) , αi ∈ (ri2, ri3). (15)
Furthermore, T i is an increasing function of λi: for λ1 > λ2 > . . . > λk and fixed FA constant d, we have T 1 > T 2 > . . . > T k.
Similar conclusions hold when θi1(0) = r i 2 − e−δ . Therefore, as the value of λi increases, the value of T i increases as well: negligible features of the data are learned sooner than the dominant ones.
On the other hand, Theorem 7 shows that such a backward incremental learning does not happen with high probability: indeed, we are requiring all the parameters Ki (which depends on the initialization) to be sufficiently negative and the probability of sampling initial conditions that will yield antiregularization depends on the value of the signals λi (i.e., it is data-dependent): the bigger the value of λi, the smaller the probability. Nevertheless, we will see in the next subsection how such a phenomenon can be safely bypassed and we discuss some properties that suggest that the initialization scheme (7) can guarantee the correct sequential learning of features as GD exhibits.
2The claim is correct if we sample the initial conditions θ1(0), θ2(0) independently and generically: for example, from a uniform distribution U([−R,R])× U([−R,R]), for some R ∈ R>0.

heading: 4.2 CASE K = 0: EVIDENCE OF INCREMENTAL LEARNING
text: Imposing Ki = 0 for all i’s means to recur to the initialization scheme (7): θi1(0) = θ i 0, θ i 2(0) =
1 2d (θ i 0) 2. In this case, Theorem 2 already provides exponential convergence rates that depend on the magnitude of the signal (|θi2(t)θi1(t)− λi| < Ce− 3 2 (dλ i) 2 3 t as t→ +∞), suggesting that the bigger the magnitude of the signal λi, the faster the convergence.
If θi0 < 0, there exists a unique time T i 0 (vanishing time) such that θ i 1(T i 0)θ i 2(T i 0) = 0 and θi1(t)θ i 2(t) ≶ 0 for t ≶ T i 0 respectively, where
T i0 = π
(ri)23 √ 3 +
1
3(ri)2 ln
( (ri − θi0)2
(θi0) 2 + riθi0 + (r i)2
) − 2
(ri)2 √ 3 arctan
( 2θi0 + r i
ri √ 3
) . (16)
We can interpret such a vanishing time as the first milestone of the learning process. In the regime as θi0 → −∞ ∀ i = 1, . . . , k, we have
T i0 = 4π
3(ri)2 √ 3 +O
( 1
(θi0) 2
) ; (17)
where ri = 3 √
2dλi, therefore for bigger λi’s (i.e. dominant singular values) the vanishing time happens sooner: for λ1 > λ2 > . . . > λk and fixed FA constant d > 0, T 10 < T 2 0 < . . . < T k 0 .

heading: 5 DISCRETE DYNAMICS

heading: 5.1 FORWARD EULER DISCRETIZATION
text: In this section, we are interested in analyzing the discrete FA dynamics and its convergence properties. For the discrete time setting, we focus again on a 2-layer linear NN. The standard forward Euler method yields the following discretized system:
W (t+1) 1 = W (t) 1 + ηM ( Σxy −W (t)2 W (t) 1 ) (18)
W (t+1) 2 = W (t+1) 2 + η ( Σxy −W (t)2 W (t) 1 )( W (t) 1 )> (19)
and performing again the SVD change of variables, the system becomes
W̃ (t+1) 1 = W̃ (t) 1 + ηD ( Λxy − W̃ (t)2 W̃ (t) 1 ) (20)
W̃ (t+1) 2 = W̃ (t) 2 + η ( Λxy − W̃ (t)2 W̃ (t) 1 )( W̃ (t) 1 )> (21)
If we set the initial conditionsW (0)1 ,W (0) 2 as diagonal, the matrices remain diagonal for all times and the dynamic decouples. For each i = 1, . . . , k, we have( θ (t+1) 1 )i = ( θ (t) 1 )i + ηdi ( λi − ( θ (t) 2 )i ( θ (t) 1
)i) (22)(
θ (t+1) 2 )i = ( θ (t) 2 )i + η ( λi − ( θ (t) 2 )i ( θ (t) 1 )i)( θ (t) 1 )i (23)
The following result holds (the index i is suppressed for simplicity).
Theorem 8. Consider the system (22)–(23). With zero initialization θ(0)1 = θ (0) 2 = 0, ∀ λ, d > 0 and with constant step size η small enough, the product θ(t)2 θ (t) 1 converges to the true signal λ with linear
rates: ∣∣∣θ(t)2 θ(t)1 − λ∣∣∣ ≤ Cqt, for some C > 0.
For the sake of clarity, we deferred the precise prescription for the magnitude of η and the linear rate q to Appendix D, where the full proof of the Theorem can be found.

heading: 5.2 MODIFIED FORWARD DISCRETIZATION
text: In addition to the standard discretization strategy, we are proposing here a new discretization technique that is particularly instrumental for the FA setting, as it will allow a thorough inspection of its convergence guarantees and convergence rates for NNs of any depth. Inspired by the fact that the behaviour of the second layer strictly depend on the choices that we impose on the first layer in the continuous dynamics (Theorem 2 and initialization (7)), we slightly modify the update rule for second weight matrixW2 in the following way:
W (t+1) 1 = W (t) 1 + ηM ( Σxy −W (t)2 W (t) 1 ) (24)
W (t+1) 2 = W (t+1) 2 + η ( Σxy −W (t)2 W (t) 1 )( W (t+ 12 ) 1 )> (25)
withW ( t+ 12 )
1 = 1 2
( W
(t+1) 1 +W (t) 1
) . The above discretization technique can be considered as a
hybrid between the Euler forward method and the trapezoidal rule (Iserles, 1996), although we are still evaluating the error Σxy −W (t)2 W (t) 1 at the present time t. The underlying idea is that we want to already use the information encoded in the new (t+ 1)-iterate ofW1 to update the second layer by averagingW (t)1 andW (t+1) 1 .
Performing again the SVD change of variables, and setting W (0)1 , W (0) 2 as diagonal, we obtain a scalar system of the form (∀ i = 1, . . . , k)( θi1 )(t+1) = ( θi1 )(t) + ηdi ( λi − ( θi2 )(t) ( θi1 )(t))
(26)( θi2 )(t+1) = ( θi2 )(t) + η ( λi − ( θi2 )(t) ( θi1 )(t)) ( θi1 )(t+ 12 ) . (27)
A convergence result follows (suppressing the i-dependence).
Theorem 9. Consider the system (26)–(27). With zero initialization θ(0)1 = θ (0) 2 = 0, ∀ λ, d > 0 and with constant step size η < 2 3(2dλ) 2 3 , the product θ(t)2 θ (t) 1 linearly converges to the true signal λ: |θ(t)2 θ (t) 1 − λ| ≤ 3λ(1− 3η 2 (2dλ) 2 3 )t.
The discretization scheme easily generalizes to the case of L layers for L ≥ 3 (we again suppress the i-index for readability):
θ (t+1) ` = θ (t) ` + ηd`θ (t+ 12 ) [1:`−1] ( λ− θ(t)L . . . θ (t) 1 ) , ` = 1, . . . , L, (28)
with dL = 1. Also in this case, the dynamics closely mimic the continuous counterpart, reducing the system to a single recurrence equation for {θ(t)1 }. By selecting a step size η small enough, the sequence of products {θ(t)L . . . θ (t) 2 θ (t) 1 } will converge to the signal λ with exponential rate.
Theorem 10. Consider the system (28) with λ, d` > 0 and zero initialization θ (0) ` = 0, ∀ ` = 1, . . . , L. If the constant step size satisfies η < (d1γ(Kλγ−1)1/γ)−1, where γ = 2L − 1 and K is a suitable positive constant depending on the FA constants d1, . . . , dL−1, then the product ∏L `=1 θ (t) ` linearly converges to the true signal λ: |θ(t)L . . . θ (t) 1 − λ| ≤ Cλ(1 − ηd1γ ( Kλγ−1 ) 1 γ )t, for some constant Cλ ∈ R>0 only dependent on λ.
The proofs of Theorems 9 and 10 can be found in Appendix D.

heading: 6 EXPERIMENTS
text: In this last section, we discuss one ML application of our theoretical results. The set-up is borrowed from Gidel et al. (2019) (Section 4.2), where a linear autoencoder is considered.
For an autoencoder, the (true) output is equal to its input y = x ∈ Rd and we want to compare the reconstruction properties ofW (t)2 W (t) 1 (2 layer NN) andW (t) 3 W (t) 2 W (t) 1 (3 layer NN) computed via
FA versus GD. In this experiment, we have d = o = 20 and k = 5. We generate n = 1000 synthetic data {xi} in the following way: each data point is given as xi = Azi + i, where A ∈ Rd×k is a fixed matrix with entries sampled as Akl ∼ U([0, 1]), zi ∼ N ( 0,Λ := diag { 4, 2, 1, 12 , 1 4 }) and the noise i ∼ 10−3N (0, Id). We initialize the weight matrices to a close-to-zero initialization W (0)1 ,W (0) 2 ,W (0) 3 ∼ 10−5U([0, 1]): this will guarantees the FA dynamics to avoid the anti-regularization pattern and to converge to the signal (Theorems 8, 9 and 10). We generate the FA matrices M (2-layers) and M1,M2 (3-layers) with entries Mij ∼ U([0, 1]). Keeping the same initial conditions W (0)` and same fixed matrix A, we repeat the FA training 15 times, sampling a different set of FA matrices each time: in Figure 2, we report the average (with error bars) of the trace norm ‖Θ(t)‖∗ with Θ(t) := W (t)2 W (t) 1 or Θ (t) := W (t) 3 W (t) 2 W (t) 1 respectively,
as well as thee average of the reconstruction errors ∥∥Θ(t) −AΛA>∥∥
2 , as a function of t the number
of iterations.
We can see that the FA learning dynamics is distinctly faster than GD, which is learning at a much slower pace, as we initialized the dynamics close to one of GD’s local maxima. In particular, the 3-layer NN trained with GD hasn’t yet started to learn at the end of the selected training period of the experiment. Notably, the FA training for the 3-layer NN hints at the presence of an implicit regularization behavior, similar to what has been observed and proved in Gidel et al. (2019) for a 2-layer NN. Additional technical details can be found in Appendix E.

heading: 7 FINAL DISCUSSIONS AND INSIGHTS
text: In the present paper we rigorously analyze a biologically plausible optimization algorithm, we describe its potential as well as its limitations in possible ML applications. We provide convergence guarantees in both the continuous and discrete-time settings; additionally, we thoroughly discuss implicit regularization phenomena that may arise with certain initialization schemes and that may deeply affect the optimization effectiveness. The occurrence of drastically opposite regularization by only modifying an algorithm’s initialization showcases the critical importance of the latter for ML applications.
Although the analysis is conducted only for linear models, the results give already numerous insights on the behaviour of nonlinear architectures. The analysis of Theorem 5 shows that, in some regimes, the deeper layers behave as "a power" of the first layer: this result suggests that each layer is converging at a different speed. It thus defines a hierarchy (in terms of speed of convergence) between the layers of the network. Furthermore, one can notice that in the linear case FA "removes" some stationary points: unlike GD that potentially has many stationary points that are saddle points, all the stationary points of FA are global minima. Therefore, we believe that FA will also "remove" some undesirable stationary points in the non-linear case. Several challenges remain, but the above pieces of intuition will be critical for the analysis in the nonlinear setting.
Beyond this, it is enticing to conduct a comparison between FA and GD and illustrate the advantages of one algorithm over the other (and in which setting). All of these directions will be the central objects of study in future works.

heading: A.1 PROOF OF PROPOSITION 1
text: Proof. We have that L := 12E[‖y − ŷ‖ 2] and ŷ = W2W1x, thus
∇W1L = E[ 12∇W1‖y − ŷ‖ 2] (29)
= E[W>2 (y −W2W1x)x>] (30) = W>2 (Σxy −W2W1Σxx) , (31)
where the last line comes from the linearity of the expectation. Similarly we have,
∇W2L = E[ 12∇W2‖y − ŷ‖ 2] (32)
= E[(y −W2W1x)x>W>1 ] (33) = (Σxy −W2W1Σxx)W>1 . (34)
Now we conclude by noticing that by definition of the FA dynamics (2), the matrix W>2 is replaced by a random matrixM .

heading: A.2 PROOF OF THEOREM 2
text: Consider the following system of ODEs for the functions θ1, θ2 ∈ C1(R):
θ̇1 = d (λ− θ2θ1) (35) θ̇2 = (λ− θ2θ1) θ1 (36)
for some d, λ ∈ R>0. From equation (35), we substitute θ̇1d = (λ − θ2θ1) in the second equation (36) to obtain
θ̇2 = 1
2d θ̇21, i.e. θ2(t) =
1
2d θ21(t) +K, (37)
for some integration constant K ∈ R. We then substitute such expression for θ2(t) back into the equation (35) in order to obtain a closed form expression for θ̇1:
θ̇1 = d
( λ− 1
2d θ31 −Kθ1
) (38)
We now choose the initial conditions for θ1(t) and θ2(t) in the following way:
θ1(0) = θ0 and θ2(0) = θ0
2
2d , (39)
so that the equation for θ1 becomes
θ̇1 = 1
2
( 2dλ− θ31 ) (40)
We are now ready to prove Theorem 2, which we recall here for convenience: Theorem 11. For any FA constant d ∈ R>0, and ∀ θ0 ∈ R initial value with initialization scheme (39), ∃ C > 0 such that the product θ2(t)θ1(t) of the solution of the system (35)–(36) converges exponentially to the signal λ and in a monotonic fashion:
|θ2(t)θ1(t)− λ| < Ce− 3 2 (dλ) 2 3 t as t→ +∞. (41)
Proof. By Picard–Lindelöf Theorem, equation (40) admits a unique local solution θ1 ∈ C1([0, T ]) (T > 0), for any initial condition θ0 ∈ R. The local solution can be then easily extended to a global solution on [0,+∞). The cubic polynomial on the right hand side of (40) has a unique real root r := 3 √
2dλ: 2dλ− θ31 = (r − θ1−) ( θ21 + θ1r + r 2 ) where the quadratic θ21 + θ1r+ r 2 is irreducible (indeed, its discriminant is ∆ = −3r2 < 0).
If θ0 = r, the solution of the ODE (40) is the constant function θ1(t) = r, ∀ t ∈ R≥0. In this case, it is trivial to see that θ2(t) = r 2
2d and θ2(t)θ1(t) = λ, ∀ t ∈ R≥0. If θ0 6= r, we simply need to analyze a separable ODE whose solution reads
− 2 3r2 ln |θ1(t)− r|+ 1 3r2 ln ( θ21(t) + rθ1(t) + r 2 ) + 2 r2 √ 3 arctan
( 2θ1(t) + r
r √ 3
) + C0 = t
(42)
with constant of integration
C0 = 2
3r2 ln |θ0 − r| −
1 3r2 ln ∣∣θ20 + rθ0 + r2∣∣− 2 r2 √ 3 arctan
( 2θ0 + r
r √ 3
) . (43)
Equation (42) is transcendental and it cannot be solved for θ1 explicitly. However, we can still analyze the behaviour of the solution.
Lemma 12. θ1(t) is monotonic: in particular, θ̇1 < 0, if θ0 > r, and θ̇1 > 0, if θ0 < r.
The above claims easily follows by inspecting the sign of the right hand side of (40). In particular, Lemma 12 implies that θ1(t) is always bounded, for any initial condition θ0 ∈ R. As t → +∞, the left hand side of (42) may only diverge in the first logarithmic term (θ1(t) → r), since the term involving the arctangent is always bounded and the other logarithmic term has argument that is always bounded (Proposition 12) and positive (the quadratic is irreducible).
Therefore, as t→ +∞
θ1(t)→ r and θ2(t)θ1(t)→ 1
2d r3 = λ (44)
for any initial condition θ0 ∈ R. See Figure 3. Furthermore, we can derive from (42) the rate of convergence:
|θ1(t)− r| =
exp − 3r2 2 t+ 3r2C0 2 + 1 2 ln |θ1(t)2 + rθ1(t) + r2|+ √ 3 sign(r) arctan ( 2θ1(t) + r |r| √ 3 ) ︸ ︷︷ ︸
=:g(t)
 (45)
where g(t) = 3r 2C0 2 + 1 2 ln(3r 2) + π√ 3 + o(1), as t→∞; therefore,
|θ1 − r| = e− 3r2 2 t+O(1). (46)
Finally, from θ2(t) = 12dθ1(t) 2, we get the desired result:∣∣∣∣θ2(t)θ1(t)− 12dr3 ∣∣∣∣ = |θ2(t)θ1(t)− λ| ≤ Ce− 3r22 t (47) for some positive constant C.
Remark 13. All the above calculations are valid also in the general case r, λ ∈ R \ {0}. In the case where r < 0 (i.e. sign(d) 6= sign(λ)), the solution θ1(t) satisfies the implicit equation
− 2 3r2 ln |θ1(t)− r|+ 1 3r2 ln ∣∣θ21(t) + rθ1(t) + r2∣∣+ 2 sign(r) r2 √ 3 arctan
( 2θ1(t) + r
|r| √ 3
) + C0 = t
(48)
and
C0 = 2
3r2 ln |θ0 − r| −
1 3r2 ln ∣∣θ20 + rθ0 + r2∣∣− 2 sign(r) r2 √ 3 arctan
( 2θ0 + r
|r| √ 3
) . (49)
To conclude the proof, we need to show that the product θ2(t)θ1(t) is monotonic increasing (if θ0 < r) or decreasing (if θ0 > r) to the limit value λ. See Figure 3 for a comparative plot of the solutions θ1(t), θ2(t) and their product.
The monotonicity properties of θ1 are described in Lemma 12). We recall here the first derivative of the functions θ1 and θ2 = 12dθ1(t) 2
θ̇1 = 1
2
( r3 − θ31 ) , θ̇2 = 1
d θ1θ̇1, (50)
and its the second derivative
θ̈1 = 3 4 θ21 ( θ31 − r3 ) , θ̈2 = 1 d ( (θ̇1) 2 + θ1θ̈1 ) (51)
If θ0 > r, then θ1(t) is strictly convex and strictly decreasing; by construction, θ2(t) is also strictly convex and decreasing. Therefore, the product θ2(t)θ1(t) will be strictly convex and strictly decreasing ∀ t ∈ R≥0.
If θ0 < r, then θ1(t) is strictly increasing and by construction the product θ2(t)θ1(t) = 12dθ 3 1 is increasing ∀ t ∈ R≥0. The description of convexity/concavity in this case is more involved and we do not pursue this direction, as it is not necessary for our results.
Case λ = 0. If we assume the signal λ = 0, we can follow the same argument as above and perform similar calculations.
Equation (40) has a simple explicit solution (assuming θ0 6= 0) θ1(t) = ± √ 1
t+ θ−20 , θ2(t) =
1
2d(t+ θ−20 ) , (52)
where ± depends on the sign of θ0. It is easy to see that
θ2(t)θ1(t)→ 0 = λ as t→ +∞, (53)
however the convergence rate is of the order of O(t− 32 ), therefore not exponential.

heading: A.3 PROOF OF THEOREM 3
text: In the setting of Theorem 3, we are considering the initial conditions of θ1 and θ2 as independent, therefore the integration constant
K = θ2(0)− 1
2d θ1(0)
2 ∈ R
may assume any value. The cubic polynomial in the ODE for θ1 has the general form
θ̇1 = − 1
2
( θ31 + 2dKθ1 − 2dλ ) (54)
and it could have 1 real root (and 2 complex conjugate roots), 2 real roots (a simple and a double root) or 3 real distinct roots. To discern the cases, we need to evaluate the discriminant:
∆ = −4d2 ( 8dK3 + 27λ2 ) . (55)
Regardless of the sign of the discriminant, there exists a unique global solution θ1 ∈ C1([0,+∞)) for any initial condition θ1(0) ∈ R and any constant K ∈ R (Picard-Lindelöf Theorem).
Case ∆ > 0. The polynomial has 3 real, distinct root r1 < r2 < r3, that may be written explicitly with the help of trigonometric functions, if needed. Note also that ∆ > 0 implies K < 0. There are three constant solutions of the equation (54): θ1(t) = rj , ∀ t ≥ 0 (provided that θ0 = rj for some j = 1, 2, 3). Otherwise (θ0 6= rj ∀ j = 1, 2, 3), the implicit solution of (54) is the following:
ln |θ1(t)− r3| (r3 − r2)(r3 − r1) − ln |θ1(t)− r2| (r3 − r2)(r2 − r1) + ln |θ1(t)− r1| (r3 − r1)(r2 − r1) + C0 = − t 2 (56)
with constant of integration
C0 = − ln |θ0 − r3|
(r3 − r2)(r3 − r1) + ln |θ0 − r2| (r3 − r2)(r2 − r1) − ln |θ0 − r1| (r3 − r1)(r2 − r1) . (57)
Lemma 14. θ1(t) is monotonic:
θ̇1 < 0 if θ0 ∈ (r1, r2) ∪ (r3,+∞), (58) θ̇1 > 0 if θ0 ∈ (−∞, r1) ∪ (r2, r3). (59)
In particular, from Lemma 14 we can conclude that the constant solutions θ1(t) = r1 and θ1(t) = r3 are attractive and θ1(t) = r2 is repulsive; the general solution (regardless of the initial value θ0 ∈ R) is bounded for all times.
Following similar arguments as in Section A.2, we can conclude that the convergence is exponentially fast:
|θ1(t)− r3| = e− (r3−r2)(r3−r1) 2 t+O(1) as t→ +∞, (60)
for θ0 ∈ (r2, r3) ∪ (r3,+∞). A similar asymptotic expansion is valid for the quantity |θ1(t)− r1| when θ0 ∈ (−∞, r1) ∪ (r1, r2). Using the relation θ2(t) = 12dθ1(t)
2 +K, we have∣∣∣∣θ2(t)θ1(t)− 12dr33 −Kr3 ∣∣∣∣ ≤ Ce− (r3−r2)(r3−r1)2 t (61)
for some C > 0. Notice that, being r3 a root of the cubic polynomial x3 + 2dx− 2dλ, it follows that 1 2dr 3 3 −Kr3 = λ. Similar arguments hold for r1.
Case ∆ < 0. The polynomial has 1 real root r ∈ R and two complex conjugate roots. Let us assume λ > 0 (therefore, r 6= 0); the case λ = 0 is analyzed shortly below. Using Cardano’s formula, the single root has the explicit expression
r = 3 √ dλ+ √ d2λ2 + 8d3K3
27 +
3 √ dλ− √ d2λ2 + 8d3K3
27 . (62)
The setting is the similar to the one analyzed Section A.2. From a monotonicity argument, it follows that the constant solution θ1(t) = r (with θ0 = r) is attractive, therefore for any given initial condition θ0, θ1(t) is bounded for all times t and it converges to r as t→ +∞.
The implicit solution reads
r
2 (r3 + dλ) ln |θ1(t)− r| −
r
4 (r3 + dλ) ln
( θ1(t) 2 + rθ1(t) + 2dλ
r ) − 3r 3
8(r3 + dλ)(dλ− r3) arctan
( r (2θ1(t) + r)
4 (dλ− r3)
) + C0 = − t
2 (63)
with C0 the constant of integration.
From θ2(t) = 12dθ1(t) 2 + K, convergence follows: θ2(t)θ1(t) → 12dr 3 + Kr = λ as t → +∞. Furthermore, convergence is exponentially fast :
|θ2(t)θ1(t)− λ| ≤ Ce− r3+dλ r t as t→ +∞, (64)
for some C > 0.
Case ∆ = 0. The vanishing of the discriminant implies that the constant of integration K has a specific value
K = −3 2 3
√ λ2
d .
Remark 15. Notice that the set of initial conditions (θ1(0), θ2(0)) such that ∆ = 0 is a set of Lebesgue measure zero in R2.
In this case, the polynomial has 1 simple root rs and one double root rd:
rs = − 3λ
K = 2
3 √ dλ, rd = 3λ 2K = − 3 √ dλ (65)
The solution of the differential equation (54) (for θ0 6= rs, rd) reads
ln |θ1(t)− rs| − ln |θ1(t)− rd| (rs − rd)2 + 1 (rs − rd)(θ1(t)− rd) + C0 = − t 2 (66)
with C0 the constant of integration.
Lemma 16. The solution θ1(t) is monotonic: if rs > rd,
θ̇1 < 0 if θ0 ∈ (rs,+∞) (67) θ̇1 > 0 if θ0 ∈ (−∞, rd) ∪ (rd, rs) (68)
and similarly if rs < rd,
θ̇1 < 0 if θ0 ∈ (rs, rd) ∪ (rd,+∞) (69) θ̇1 > 0 if θ0 ∈ (−∞, rs). (70)
Therefore, the solution θ1(t) = rs is always attractive, while the solution θ1(t) = rd is "one-sided" attractive. The solution θ1(t) is bounded for all times and the convergence result follows:
θ2(t)θ1(t)→ 1
2d r3s +Krs = λ as t→∞; (71)
similarly, for rd. Exponential convergence holds when θ1(t)→ rs, while when θ1(t)→ rd the rate is sub-exponential.
In the special case when ∆ = 0 and λ = 0, the polynomial has a triple root r = 0. The solution is the same as the one derived in Section A.2 in the case of λ = 0 and the convergence to the signal λ = 0 is polynomial:
|θ2(t)θ1(t)| ≤ Ct− 3 2 as t→ +∞, (72)
for some C > 0.
B CONVERGENCE OF FA FOR L-LAYER NNS

heading: B.1 PROOF OF PROPOSITION 4
text: We have that L := 12E[‖y − ŷ‖ 2] and ŷ = WL · · ·W2W1x, thus noting W[a:b] := WbWb−1 · · ·Wa we have, ∇W`L = E[ 12∇W`‖y − ŷ‖ 2] (73)
= E[W>[`+1:L](y −WL · · ·W1x)x >W>[1:i−1]] (74) = W>[`+1:L](Σxy −WL · · ·W1Σxx)W > [1:`−1] , (75)
where the last line comes from the linearity of the expectation. Now we conclude by noticing that by definition of the FA dynamics (2), the matrices on the left hand-side of (y −W2W1x) are replaced by random matricesMl. Thus we get
Ẇ` = M`(Σxy −W[1:L]Σxx)W>[1:`−1] ` ∈ [L] , (76)
where we reparametrizedM` = M[`+1:L] (see discussion in the main text).

heading: B.2 PROOF OF THEOREM 5
text: We will first state Theorem 5 in a more extensive form. Given the system
θ̇L = (λ− θL . . . θ1)θL−1 . . . θ1 (77) θ̇L−1 = dL−1(λ− θL . . . θ1)θL−2 . . . θ1 (78) ... (79)
θ̇2 = d2(λ− θL . . . θ1)θ1 (80) θ̇1 = d1(λ− θL . . . θ1) (81)
where we suppressed the index i for simplicity, the following result holds. Theorem 17. For any set of FA constants {d`}`=1,...,L−1, d` ∈ R>0, there exists an initialization scheme such that the system of differential equations (77)–(81) can be reduced into an equation for θ1 of the form
θ̇1 = d1 (λ− Kθγ1 ) (82)
for some constant K ∈ R>0 that depends on {d`}, with γ = 2L − 1. The other functions depend directly on θ1 in a power-like fashion:
θ`(t) = C`
2`−1 θ1(t)
2`−1 , ` = 2, . . . , L, (83)
for some C` ∈ R>0. Furthermore, the following convergence result holds: θL(t)θL−1(t) . . . θ1(t)→ λ as t→ +∞. (84)
See Figure 4 for a numerical simulation of a 3-layer scalar network.
Proof. The proof of reducing the system of differential equations into a single equation for θ1 will be constructive, following a simple iterative procedure. We start by considering the last equation (81) and by substituting (λ− θL . . . θ1) = θ̇1d1 in the second to last equation
θ̇2 = d2 d1 θ1θ̇1 = d2 2d1
( θ21 )· , implying θ2(t) =
d2 2d1 θ1(t) 2 +K2, (85)
for some constant of integration K2 ∈ R, which can be set to zero if we impose θ2(0) = d2d1 θ1(0) 2.
The same substitution strategy can be consecutively applied to the third to last equation, fourth to last equation and so on. In general, it is easy to prove by induction that ∀ ` = 2, . . . , L we have θ̇` = C`θ 2`−1−1 1 θ̇1, for some C` = C`(d1, . . . , d`) ∈ R>0, which implies
θ`(t) = C`
2`−1 θ1(t)
2`−1 +K` (86)
where as before we can set the constant of integration K` = 0 by imposing θ`(0) = C`2`−1 θ1(0) 2`−1 .
Once we have derived the relations between the functions {θ`(t)}`=2,...,L and θ1(t), we can plug all of them into (81) to obtain a closed-form differential equation for θ1:
θ̇1 = d1 (λ− Kθγ1 ) , (87) where K = ∏L `=1 C` 2`−1 ∈ R>0 and γ = ∑L−1 `=0 2
` = 2L − 1. This is a separable equation which can be solved (with tears) by using partial fraction decomposition.
The proof of convergence of the product ∏L `=1 θ`(t) to the true signal λ follows a similar argument
as in Appendix A.2. Let r := γ √
λ K be the unique real root of the polynomial λ − Kθ γ 1 , then the
solution θ1(t) = r is a constant solution, provided that θ1(0) = r. By construction, it follows that {θ`(t)}`=2,...,L are also constant functions and their product is such that
θL(t) . . . θ1(t) = L∏ `=1 C` 2`−1 rγ = λ ∀ t ≥ 0. (88)
Lemma 18. The general solution θ1(t) is monotonic: θ̇1 > 0 if θ1(0) < r (89)
θ̇1 < 0 if θ1(0) > r (90)
Therefore, θ1(t) is bounded for all times and θ1(t)→ r as t→ +∞. By construction, it follows that θL(t) . . . θ1(t)→ ∏ ` C` 2`−1 rγ = λ, as t→ +∞. (91)
Exponential convergence follows from a careful analysis of the implicit solution of the separable equation (87), in a similar way as it has been done in Appendix A.2.
The proof of convergence of the FA algorithm for a generic initialization (θ1(0), . . . , θL(0)), where each layer is initialized independently, follows the same guidelines as the 2-layer case (Appendix A.3).
Sketch of the proof. In the same setting as the proof above, by keeping all the constants of integration K` 6= 0, we reduce the system of ODEs (81) to a single ODE for θ1 of the form
θ̇1 = Pγ(θ1; d1, . . . , dL, λ) (92) where P is a polynomial of degree γ := 2L − 1 in the variable θ1 with coefficients depending on the FA constants d1, . . . , dL and the signal λ. It is still a separable equation that can be studied via a careful analysis of its discriminant.
C IMPLICIT REGULARIZATION PHENOMENA

heading: C.1 PROOF OF THEOREM 7 (CASE ∆ > 0)
text: For readability, we are suppressing the index i from all the formulæ below. Recall the Initial Value Problem for θ1 (and θ2):
θ̇1 = − 1
2
( θ31 + 2dKθ1 − 2dλ ) (93)
θ1(0) = θ0 (94)
and θ2(t) = 12dθ 2 1(t) +K, where K = θ2(0)− 12dθ 2 0 . By assuming that the discriminant is positive ∆ = −4d2 ( 8dK3 + 27λ2 ) > 0, the cubic polynomial on the right hand side of (93) has 3 real distinct root r1 < r2 < r3 and the solution of the ODE is (see also Section A.3)
ln |θ1(t)− r3| r32r31 − ln |θ1(t)− r2| r32r21 + ln |θ1(t)− r1| r31r21 + C0 = − t 2 (95) C0 = − ln |θ0 − r3| r32r31 + ln |θ0 − r2| r32r21 − ln |θ0 − r1| r31r21 , (96)
where we set rij := ri − rj .
If we set as initial value θ0 = r2 + e−δ (δ > 0), then
C0 = − δ
r32r21 −
ln ( r32 − e−δ ) r32r31 − ln ( r21 + e −δ) r31r21 ; (97)
we additionally rescale the time variable as t 7→ δt. Rewriting equation (95), we have r3 − θ1(δt) = exp { −r32r31
2 δ
[ t− 2
r32r21
] + ln ( r32 − e−δ ) + r32 r21 ln ( r21 + e −δ) + r31 r21 ln (θ1(δt)− r2)− r32 r21 ln (θ1(δt)− r1) } , (98)
θ1(δt)− r2 = exp { −r32r21
2 δ
[ 2 r32r21 − t ] + r21 r31 ln ( r3 − θ1(δt) r32 − e−δ ) + r32 r31 ln ( θ1(δt)− r1 r21 + e−δ )} .
(99)
From (98)–(99), we have the following asymptotic expansions: as δ → +∞ θ1(δt)→ r3 for t > T (100)
θ1(δT )→ α := r3 − exp {(
1 + r31 r21
) ln r32 +
r32 r21 ln ( r21 r31 )} (101)
θ1(δt)→ r2 for t < T (102)
with T = 2r32r21 .
In the case θ0 = r2 − e−δ , we can follow a similar argument: as δ → +∞, θ1(δt)→ r1 for t > T (103)
θ1(δT )→ α̃ := r1 + exp {(
1 + r31 r32
) ln r21 +
r21 r32 ln ( r32 r31 )} (104)
θ1(δt)→ r2 for t < T. (105)
The last bit of result we need is to analyze if the value of the critical value
T = 2
r32r21 (106)
is monotonic as λ increases. Recall the cubic polynomial on the right hand side of equation (93) (without loss of generality, we are setting 2d = 1): P (x) = x3 +Kx−λ = (x−r1)(x−r2)(x−r3), where λ > 0 and K < − 32λ 2/3 < 0 depends on the initial conditions of the FA system.
(a) Numerical simulation of θ1(t) with r1 = −2, r2 = 1, r3 = 2 and initial conditions θ0 = r2 + e−δ with δ ∈ {10, 20, 30}. (b) Plot of θ1(δt) (∆ > 0), with initial condition θ0 = r2 + e −δ , δ = 20), FA constant d = 0.5, initial condition K = −4.
Figure 5: (Anti)-implicit regularization.
The quantity r32r21 can be expressed as r32r21 = −P ′(r2) = −3r22 −K, where r2 < 0 for λ > 0. It is easy to see that if λ < λ̃, then |r2| < |r̃2| and consequently r32r21 = −r22 −K > −r̃2
2 −K = r̃32r̃21, i.e.
Tλ = 2
r32r21 <
2
r̃32r̃21 = Tλ̃ (107)
Therefore the threshold time T is an increasing function of λ. See Figure 5b: in the δ-scaling limit we see indeed that for smaller values of λ, the step appears sooner.

heading: C.2 EVIDENCE OF INCREMENTAL LEARNING IN THE CASE K = 0
text: Consider again the system (93) in the case K = 0:
θ̇1 = − 1
2 (θ31 − 2dλ) (108)
θ1(0) = θ0 (109)
and θ2(t) = 12dθ 2 1(t). In this case, the solution reads (see also Section A.2)
− 2 3r2 ln |θ1(t)− r|+ 1 3r2 ln ( θ21(t) + rθ1(t) + r 2 ) + 2 r2 √ 3 arctan
( 2θ1(t) + r
r √ 3
) + C0 = t
(110)
C0 = 2
3r2 ln |θ0 − r| −
1 3r2 ln ( θ20 + rθ0 + r 2 ) − 2 r2 √ 3 arctan
( 2θ0 + r
r √ 3
) (111)
with r := 3 √ 2dλ.
As the limit value of the function θ1(t) is the positive root r, θ1(T0) = 0 (for some T0 > 0) if and only if θ0 < 0. The value of T0 can be easily derived from the implicit expression (110):
T0 = π
r23 √ 3 +
1
3r2 ln
( (r − θ0)2
θ20 + rθ0 + r 2 ) − 2 r2 √ 3 arctan ( 2θ0 + r r √ 3 ) (112)
and in the regime as θ0 → −∞
T0 = 4π
3r2 √ 3 − 1 θ20
+O ( 1
θ30
) , θ0 → −∞. (113)
It is now clear that for bigger λ’s (i.e. dominant singular values) the crossing of the x-axis happens sooner: for λ > λ̃ and fixed FA constant d > 0,
T0,λ < T0,λ̃. (114)

heading: D CONVERGENCE OF THE DISCRETE ALGORITHMS
text: We will report here all the proofs of the convergence theorems for the discrete FA dynamics. We will first discuss convergence guarantees and convergence rates for the modified (mid-point) FA dynamics for both 2- and L-layer NNs (L ≥ 3), as it is more straightforward. The convergence proof for the forward Euler discretization method can be found in Section D.2.

heading: D.1 CONVERGENCE OF FA UNDER THE "MIDPOINT" DISCRETIZATION SCHEME
text: We recall the 2-layer system that we want to analyze. For simplicity we are suppressing all the i indices and we are denoting (xt, yt) the unknown functions in the system:
xt+1 = xt + ηd (λ− xtyt) (115) yt+1 = yt + η
2 (λ− xtyt) (xt+1 + xt) (116)
We begin with proving a series of easy lemmas that will be instrumental for the proof of the theorem. Lemma 19. With initial condition (x0, y0) = (0, 0), the second component of the solution of the system (115)–(116) reads
yt = x2t 2d
∀ t ≥ 0. (117)
Proof. From equation (115), xt+1−xtηd = λ − xtyt; substituting in equation (116), we obtain the result:
yt+1 = yt + (xt+1 + xt)(xt+1 − xt)
2d = yt + x2t+1 − x2t 2d
= y0 + t∑ j=0 x2j+1 − x2j 2d = y0 − x20 2d︸ ︷︷ ︸
=0
+ x2t+1 2d
(118)
via a telescoping sum.
Substituting the expression of yt back in equation (115), we have a recurrence relation for xt of the form:
xt+1 = xt + ηdλ
( 1− x 3 t
r3
) =: f(xt). (119)
where r := 3 √
2dλ. Lemma 20. With parameter
η < 2
3r2 =
2
3(2dλ) 2 3
, (120)
the function f(x) = x+ ηdλ ( 1− x 3
r3
) maps the interval [0, r] into itself: ∀ x ∈ [0, r], f(x) ∈ [0, r].
Proof. By inspecting the function f and its derivative we can conclude that f(x) ≥ 0, ∀ x ∈ [0, r] and f is increasing on the interval [ 0, √
2 3η
] . Notice also that f(r) = r. Imposing (120) guarantees
that r < √
2 3η , therefore
ηdλ = f(0) ≤ f(x) ≤ f(r), ∀ x ∈ [0, r].
Lemma 21. In the same hypotheses as Lemmas 19 and 20, the sequence {xt}t≥0 is bounded, positive and increasing:
0 ≤ xt ≤ r and xt ≤ xt+1, ∀ t ≥ 0. (121)
Proof. For t = 0, trivially 0 ≤ x0 = 0 ≤ r. By induction, assume that 0 ≤ xt ≤ r, then by Lemma 20
0 ≤ xt+1 = xt + ηd ( λ− x 3 t
2d
) = f(xt) ≤ f(r) = r. (122)
Furthermore,
xt+1 = xt + η
2
( 2dλ− x3t ) ≥ xt + η
2
( 2dλ− ( 3 √ 2dλ)3 ) = xt. (123)
We can now prove Theorem 9 that we rewrite here for convenience: Theorem 22. With initial conditions (x0, y0) = (0, 0), ∀ λ, d > 0 and with step size η < 2
3(2dλ) 2 3
,
the following convergence result hold: |xtyt − λ| ≤ 3λ (
1− η 2 (2dλ) 2 3
)t . (124)
Remark 23. Note that if η = 2 (2dλ) 2 3 (which is outside the range prescribed in Theorem 9) and x0 = 0, we achieve convergence in one step: xt = 3 √ 2dλ ∀ t > 0.
Proof. By Lemma 21, ∃ ` ∈ R≥0 such that xt → ` as t→ +∞. Furthermore, by solving
` = `+ ηdλ ( 1− ` 3
r3
) (125)
it follows that ` = r = 3 √ 2dλ. Finally, we have
|xt − r| = ∣∣∣xt−1 − r + η
2
( r3 − x3t−1 )∣∣∣ = |xt−1 − r| ∣∣∣1− η 2 ( x2t−1 + rxt−1 + r 2 )∣∣∣
≤ |xt−1 − r| ( 1− 3η 2 r2 )
≤ r ( 1− 3η 2 r2 )t
(126)
Note that 0 < 1 − 3η2 r 2 < 1, since η < 23r2 . In conclusion, thanks to the arguments above and Lemma 19, we have: as t→∞
|xtyt − λ| = ∣∣∣∣x3t2d − λ ∣∣∣∣ = 12d ∣∣x2t + rxt + r2∣∣ |xt − r| ≤ 3r 2
2d |xt − r| ≤
3r3
2d
( 1− 3η 2 r2 )t
= 3λ ( 1− 3η 2 r2 )t . (127)
For the general L-layer case, we first recall Theorem 10. Theorem 24. Consider the following system:
θ (t+1) L = θ (t) L + ηθ
( t+
1 2 ) L−1 . . . θ ( t+ 1 2 ) 1 ( λ− θ(t)L . . . θ (t) 1 ) (128)
θ (t+1) L−1 = θ (t) L−1 + ηdL−1θ
( t+
1 2 ) L−2 . . . θ ( t+ 1 2 ) 1 ( λ− θ(t)L . . . θ (t) 1 ) (129)
... (130)
θ (t+1) 2 = θ (t) 2 + ηd2θ
( t+
1 2 ) 1 ( λ− θ(t)L . . . θ (t) 1 ) (131)
θ (t+1) 1 = θ (t) 1 + ηd1 ( λ− θ(t)L . . . θ (t) 1 ) (132)
with zero initialization θ(0)` = 0, ∀ ` = 1, . . . , L. Assume that λ, d` > 0, ∀ ` = 1, . . . , L− 1, and the constant step size satisfies
η < 1
d1γ (Kλγ−1) 1 γ
where γ = 2L − 1 and K is a suitable positive constant depending on the FA constants d1, . . . , dL−1. Then, the product ∏L `=1 θ
(t) ` linearly converges to the true signal λ:∣∣∣θ(t)L . . . θ(t)1 − λ∣∣∣ ≤ Cλ (1− ηd1γ (Kλγ−1) 1γ )t (133)
for some constant Cλ ∈ R>0 only dependent on λ.
Proof. We consider the last equation (132) and substitute
λ− θ(t)L . . . θ (t) 1 =
θ (t+1) 1 − θ (t) 1
ηd1 (134)
in the second to last one to obtain
θ (t+1) 2 = θ (t) 2 + d2 2d1
[( θ (t+1) 1 )2 − ( θ (t) 1 )2] ⇒ θ(t)2 =
d2 2d1
( θ (t) 1 )2 (135)
by telescoping sum. The same substitution strategy can be consecutively applied to the third to last equation, fourth to last equation and so on. In general, it is easy to prove by induction that ∀ ` = 2, . . . , L we have
θ (t+1) ` = θ (t) ` + C` [( θ (t+1) 1 )2`−1 − ( θ (t) 1 )2`−1] ⇒ θ(t)` = C` ( θ (t) 1 )2`−1 (136)
for some C` = C`(d1, . . . , d`) ∈ R>0.
Once derived the relations between the sequences {θ(t)` }`=2,...,L and {θ (t) 1 }, we can go back to the recurrence equation (132) for θ(t)1 to obtain
θ (t+1) 1 = θ (t) 1 + ηd1
[ λ− K ( θ (t) 1 )γ] (137)
where K = ∏L `=1 C` ∈ R>0 and γ = ∑L−1 `=0 2 ` = 2L − 1.
Lemma 25. With initial condition θ(t)1 = 0 and for η < 1
d1γ(Kλγ−1) 1 γ
, the sequence {θ(t)1 } is positive, increasing and bounded by r := γ √
λ K .
Proof. The proof follows the same argument as in the 2-layer case. The odd polynomial f(x) = x+ ηd1 (λ− Kxγ) has at least one fixed points at x = r and f(0) = ηd1λ > 0. Furthermore, f([0, r]) ⊆
[0, r]. Indeed, the function is increasing on the interval [0, xmax] with xmax = (
1 ηd1γK
) 1 γ−1
. If
η < 1 d1γ(Kλγ−1) 1 γ then r < xmax, implying that
0 < ηd1λ = f(0) ≤ f(x) ≤ f(r) = r, ∀ x ∈ [0, r].
It follows (by induction) that 0 ≤ θ(t)1 ≤ r ∀ t ≥ 0 and θ (t) 1 ≤ θ (t+1) 1 .
From Lemma 25 we can conclude that θ(t)1 → r and∣∣∣θ(t)1 − r∣∣∣ = ∣∣∣θ(t−1)1 − r − ηd1K((θ(t−1)1 )γ − rγ)∣∣∣ = ∣∣∣θ(t−1)1 − r∣∣∣ ∣∣∣∣∣∣1− ηd1K γ−1∑ j=0 ( θ (t−1) 1 )γ−1−j rj ∣∣∣∣∣∣ ≤ ∣∣∣θ(t−1)1 − r∣∣∣ ∣∣1− ηd1Kγrγ−1∣∣ = ∣∣∣θ(t−1)1 − r∣∣∣ ( 1− ηd1Kγ ( λ
K
) γ−1 γ )
≤ r ( 1− ηd1γ ( Kλγ−1 ) 1 γ )t . (138)
Finally, we have the following linear rate of convergence∣∣∣θ(t)L . . . θ(t)1 − λ∣∣∣ = ∣∣∣K(θ(t)1 )γ − λ∣∣∣ ≤ Cλ ∣∣∣θ(t)1 − r∣∣∣ ≤ C̃λ (1− ηd1γ (Kλγ−1) 1γ )t . (139)

heading: D.2 DISCRETIZATION UNDER THE EULER FORWARD METHOD
text: By discretizing the FA dynamics via the standard Euler forward method, we obtain the following recurrence relation for the iterates:
xt+1 = xt + ηd (λ− xtyt) (140) yt+1 = yt + ηxt (λ− xtyt) (141)
with initial condition x0 = y0 = 0. If we use the same substitution strategy as before (λ− xtyt = xt+1−xt
ηd ), we get
yt+1 = yt + ηxt xt+1 − xt
ηd = yt + xt d (xt+1 − xt)
= yt + 1
2d [(xt+1 + xt)− (xx+1 − xt)] (xt+1 − xt) = yt +
1
2d
( x2t+1 − x2t ) − 1
2d (xt+1 − xt)2
= y0 − x20 2d︸ ︷︷ ︸
=0
+ 1
2d x2t+1 −
1
2d t∑ j=1 (xj − xj−1)2, (142)
i.e.
yt = 1
2d x2t −
1
2d St ∀ t ≥ 0, (143)
where St := ∑t i=1 (xi − xi−1)
2 and S0 = 0. Notice that 0 ≤ St ≤ St+1 ∀ t ≥ 0, by definition. Plugging (143) back into equation (140), we get a highly nonlinear recurrence relation where all the past terms of the sequence (up to time t) are involved in determining the subsequent term xt+1:
xt+1 = xt + η
2
( 2dλ− x3t + xtSt ) . (144)
We will now provide a proof of Theorem 8. The proof is quite long and it is divided into four parts for readability. We will first rewrite the statement in the following way:
Theorem 26. Consider the sequence {(xt, yt)} defined recursively in (143)-(144) with zero initialization x0 = y0 = 0. For any λ, d > 0, set the constant step size as
η < min
{ 2
3(S∗ + 1)2 ,
2
max(x,s)∈R P (x, S)
} , (145)
where P (x, S) = 2dλ− x3 + xS, S∗ is the unique positive solution of P (x, x) = 0, and R is the following compact, convex set
R := { (x, S) ∈ R2≥0 | S ≤ x, P (x, S) ≥ 0 } . (146)
Then, the product xtyt linearly converges to the true signal λ:
|xtyt − λ| ≤ Cqt. (147)
for some C > 0 and 0 < q < 1.
Proof of convergence - part 1. In order to prove convergence of the sequence {xt}, we will first consider a more general sequence defined by the following recurrence relation
xt+1 = xt + η
2
( 2dλ− x3t + xtαt ) (148)
where {αt} is an auxiliary sequence that we assume to be increasing, positive and bounded (therefore convergent: αt → α∞ for some α∞ ∈ R>0). Consider the following fixed-point problem:
x = f(x;α) with f(x;α) = x+ η
2
( 2dλ− x3 + αx ) (149)
with x ≥ 0 and α ∈ [0, α∞]. Notice that ∀ α ≥ 0, we have that f(0;α) = ηdλ > 0 and f(x;α)→ −∞ as x→ +∞. Therefore, ∀ α ≥ 0, there exists a unique fixed point `α > 0 for the function f (uniqueness follows by directly inspecting the cubic polynomial f(·;α) with α fixed).
From ∂f∂α = η 2x ≥ 0 for x ≥ 0, it follows that f(x;α) ≤ f(x;α ′) ∀ α ≤ α′. Furthermore, `α ≤ `α′ : indeed, ∀ α ≥ 0, `α is the unique zero of the cubic polynomial Pα(x) = −x3 + αx + 2dλ on the positive real line (Pα(0) = 2dλ > 0, P ′α(0) = α > 0 and Pα(x) → −∞ as x → −∞); its partial derivative ∂Pα∂α = x ≥ 0, for x ≥ 0, implies that as α increases the zero `α shifts rightwards on R. On the other hand, from
∂f ∂x = 1 + η 2 α− 3 2 x2 ≥ 0 (150)
we obtain that the function f is increasing on the interval [ 0, √
2 3η + α 3
] . In order to ensure mono-
tonicity on the interval [0, `α], we tune the parameter η such that `α ≤ √ 2
3η + α 3 . (151)
By taking the supremum over α ∈ [0, α∞] on the left hand side and the infimum on the right hand side, we have
`α∞ ≤ √ 2
3η ; (152)
such a bound is satisfied if
η ≤ 2 3`2α∞ . (153)
As it will be clear in the next parts of the proof (see Remark 29), the value of `α∞ is smaller than S ∗, therefore (145) guarantees that the above bound (153) is satisfied. This implies that
1. f(x;α) ≤ f(x′;α) ∀ x, x′ ∈ [0, `α∞ ], x ≤ x′
2. f(x;α) maps [0, `S∞ ] into itself ∀ α ∈ [0, α∞].
In conclusion, we have
f(x;α) ≤ f(x′;α′) ∀ x ≤ x′,∀ α ≤ α′. (154)
We are now ready to prove the following lemma.
Lemma 27. For η ≤ 23`2α∞ , the sequence (148) (with zero initial condition) is positive, increasing and bounded: ∀ t ≥ 0
0 ≤ xt ≤ xt+1 ≤ `α∞ (155)
Proof. The proof is by induction: for t = 0 we have 0 = x0 ≤ x1 = ηdλ ≤ `α∞ and assuming that xt−1 ≤ xt, it follows that
0 ≤ xt = f(xt−1;αt−1) ≤ f(xt;αt) = xt+1 ≤ f(`α∞ ;α∞) = `α∞ (156)
Therefore the sequence (148) converges: xt → ω for some ω ∈ R>0. By construction, it is clear that ω = `α∞ . Indeed, given
xt+1 = xt + η
2
( 2dλ− x3t + xtαt ) (157)
and taking the limit as t→∞ on both sides of the equation, we obtain
ω = ω + η
2
( 2dλ− ω3 + ωα∞ ) , (158)
i.e. ω is a (positive) solution of the fixed point problem x = f(x;α∞) whose unique solution is `α∞ . A major problem arises when we consider the case of αt = St = ∑t i=1(xi−xi−1)2, i.e. the auxiliary sequence {αt} depends directly on the primary sequence {xt}. Clearly, 0 ≤ αt ≤ αt+1 for all times t ≥ 0, but the boundedness property is not straightforward.
Boundedness of the sequence of partial sums {St}. Since St+1 = St + (xt+1 − xxt) 2, consider the recurrence relations:
xt+1 = xt + η
2
( 2dλ− x3t + xtSt ) (159)
St+1 = St + η2
4
( 2dλ− x3t + xtSt )2 (160)
Lemma 28. Given the convex, compact set R := { (x, S) ∈ R2≥0 | S ≤ x, P (x, S) ≥ 0 } , (161)
with P (x, S) = 2dλ− x3 + xS, and
η < min
{ 2
3(S∗ + 1)2 ,
2
max(x,s)∈R P (x, S)
} , (162)
where S∗ is the unique positive solution to P (x, x) = 0, and with initial values x0 = S0 = 0, the sequence {(xt, St)} defined recursively in (159)-(160) is bounded and lies in the setR.
See Figure 6 for a sketch of the regionR.
Proof. The set of fixed points of the system is the set of zeros of the function P (x, S),{ x = x+ η2P (x, S)
S = S + η 2
4 P (x, S) 2
⇔ P (x, S) = 2dλ− x3 + xS = 0, (163)
which is an algebraic curve in R2, with solutions in (x, S) ∈ R2≥0
S = x2 − 2dλ x ; (164)
furthermore, there exists a unique solution (S∗, S∗) ∈ R2>0 such that P (S∗, S∗) = 0: indeed, S∗ is the unique solution to the equation x3 − x2 + 2dλ = 0 and it’s explicit expression can be recovered from Cardano’s formula (if needed).
We will prove the theorem by induction. For t = 0, trivially 0 = S0 = x0 and P (x0, S0) = 2dλ > 0. For t = 1,
S1 = η 2d2λ2 ≤ x1 = ηdλ if η ≤
1
dλ (165)
P (x1, S1) = 2dλ− (ηdλ)3 + ηdλ · (ηdλ2) = 2dλ > 0. (166)
It is easy to notice that if (162) is satisfied, then η ≤ 1dλ as required in (165). SinceR is compact, the value of the polynomial P (x, S) is bounded and positive for all (x, S) ∈ R:
0 < 2dλ = P (x0, S0) ≤ max (x,s)∈R P (x, S) < +∞.
Assume that ∀ s = 1, . . . , t we have
Ss ≤ xs, P (xs, Ss) = 2dλ− x3s + xsSs ≥ 0 (167)
then,
St+1 = St︸︷︷︸ ≤xt + η2 4 P (xt, St) 2 ≤ xt + η 2 P (xt, St) · η 2 P (xt, St)︸ ︷︷ ︸ ≤1 ≤ xt + η 2 P (xt, St) = xt+1,
(168)
thanks to the induction step and (162). Additionally, P (xt+1, St+1) = P (xt, St) +∇P (x̄, S̄) · [ η 2P (xt, St) η2
4 P (xt, St) 2
] (169)
for some (x̄, S̄) belonging to the line segment connecting the point (xt, St) to the point (xt+1, St+1), by the Mean Value Theorem;
P (xt+1, St+1) = P (xt, St) + η
2 P (xt, St)
( −3x̄+ S̄ ) + η2
4 x̄P (xt, St)
2
= P (xt, St)
[ 1 + η
2 (−3x̄2 + S̄) + η
2
4 x̄P (xt, St) ] ≥ P (xt, St) [ 1− 3η 2 x̄2 ] . (170)
Notice now that 0 ≤ x̄ = xt + δ(xt+1 − xt) = xt + δ η2P (xt, St), for some δ ∈ [0, 1]; since (xt, St) ∈ R by induction,
x̄ ≤ S∗ + 1 · η 2 max (x,S)∈R
P (x, S)︸ ︷︷ ︸ ≤1
≤ S∗ + 1. (171)
Therefore, P (xt+1, St+1) ≥ P (xt, St) [ 1− 3η 2 x̄2 ] ≥ P (xt, St) [ 1− 3η 2 (S∗ + 1) 2 ] ≥ 0, (172)
by (162). The result follows.
Proof of convergence - part 2. Theorem 28 above implies that the sequence of partial sums {St = ∑t i=1(xi − xi−1)2} converges: St → S∞ := ∑∞ i=0(xi − xi−1)2 < +∞ and we can resort to the arguments in the first part of the proof to conclude that xt → `S∞ .
Remark 29. Note that `S∞ is the unique positive solution to the equation P (x, S∞) = 0 and, since the sequence {(xt, St)} ⊂ R, `S∞ lies on the boundary ofR. This implies that `S∞ ≤ S∗.
Finally, the product xtyt → λ as t → +∞. Indeed, since xt → `S∞ =: L, then, by construction yt will also converge to some value L̂ (because of (143)). Moreover, L is the fixed point of the recurrence equation (144):
L = L+ η 2
( 2dλ− L3 + LS∞ ) , (173)
which can be alternatively written as L = L+ ηd ( λ− LL̂ ) ; (174)
L is fixed point if and only if LL̂ = λ, implying xtyt → λ.
Linear rates. The last piece of results that needs to be proven is the linear rate of convergence. Consider the difference `St − xt, where `St is the unique positive solution to P (x, St) = 0 and by construction xt ≤ St:
`St+1 − xt+1− = `St+1 − xt − η
2 P (xt, St)
= `St+1 − xt − η
2 (`St − xt)
( x2t + `Stxt + 2dλ
`St ) ≤ (`St − xt) [ 1− η
2
( x2t + xt`St + 2dλ
`t )] ︸ ︷︷ ︸
=:B1
+ ( `St+1 − `St )︸ ︷︷ ︸ =:B2
(175)
We will bound each term B1 and B2 separately. The first term can be easily estimated:
B1 ≤ 1− η
2
( (S∗)2 + S∗`S∞ + 2dλ
`0
) = 1− ηM (176)
where we defined
M := (S∗)2 + S∗`S∞ + (2dλ) 2 3
2 (177)
(recall that `S0 = 3 √ 2dλ).
In order to estimate B2, we will first prove a series of useful lemmas that illustrate some properties of the sequence {`St}. Lemma 30. ∀ t ≥ 0,
0 < 2dλ
`S∞ ≤ `2St − St ≤
2dλ `S0 . (178)
Proof. By definition of `St ,
P (`t, St) = 2dλ− `3St + `StSt = 0 (179)
i.e.
`2St − St = 2dλ
`St > 0. (180)
Using the fact that {`St} is increasing and bounded (`S0 ≤ `St ≤ `S∞ , ∀ t ≥ 0), the result follows.
Lemma 31. ∀ t ≥ 0,
xt+1 − xt ≤ ηM (`St − xt) (181)
with M defined in (177).
Proof. The proof easily follows from the recurrence relation of the (increasing) sequence {xt}:
xt+1 − xt = η
2
( 2dλ− x3t + xtSt ) ≤ (`St − xt) · η
2
( (S∗)2 + S∗`S∞ + (2dλ) 2 3 ) ︸ ︷︷ ︸
=ηM
(182)
Lemma 32. ∀ t ≥ 0,
`St+1 − `St ≤ C∞ (xt+1 − xt) 2 (183)
for some constant C∞ ∈ R>0.
Proof. Using again the definition of `St , we have
`3St+1 − ` 3 St = `St+1St+1 − `StSt = ( `St+1 − `St ) St+1 + `St (St+1 − St)
= ( `St+1 − `St ) St+1 + `St (xt+1 − xt) 2 ; (184)
on the other hand,
`3St+1 − ` 3 St = ( `St+1 − `St ) ( `2St+1 + `St+1`St + ` 2 St ) ; (185)
combining (184) and (185) and using Lemma 30, we obtain
`St+1 − `St = `St (xt+1 − xt)
2
`2St+1 + `St+1`St + ` 2 St − St+1
≤ `S∞ (xt+1 − xt) 2
2`2S0 + 2dλ `S∞
= `S∞
2(2dλ) 2 3 + 2dλ`S∞︸ ︷︷ ︸
=:C∞
(xt+1 − xt)2 . (186)
Collecting the results from Lemmas 31 and 32, we have
B2 ≤ C∞(xt+1 − xt)2 ≤ η2MC∞ (`St − xt) 2 (187)
In conclusion, for η < 1M we have
`St+1 − xt+1 ≤ (`St − xt) (1− ηM) + η2M2C∞ (`St − xt) 2 = (`St − xt) ( 1− ηM + η2M̃ )
≤ 3 √ 2dλ ( 1− ηM + η2M̃ )t
(188)
where M̃ := 2`S∞M 2C∞ > 0 (we used the fact that {`St} and {xt} are bounded); we stress that the quantity 1 − ηM − η2M̃ is smaller than 1 (thus being indeed a convergence rate) for η small enough.

heading: E EXPERIMENTS

heading: E.1 SOLUTION OF THE FA ODE SYSTEM
text: Figures 1, 3, 4, 5, plotting solutions (and products of solutions) of the FA system
θ̇L = (λ− θL . . . θ1)θL−1 . . . θ1 (189) θ̇L−1 = dL−1(λ− θL . . . θ1)θL−2 . . . θ1 (190) ... (191)
θ̇2 = d2(λ− θL . . . θ1)θ1 (192) θ̇1 = d1(λ− θL . . . θ1) (193)
in the continuous setting (t ∈ R≥0) are obtained using the standard ODE solver ode45 on Matlab (version R2020b).

heading: E.2 LINEAR AUTOENCODERS
text: We set up two linear autoencoders
ŷ = W2W1x (2-layer NN) (194) ŷ = W3W2W1x (3-layer NN (195)
whereW` ∈ Rd×d, d = 20 and the data are synthetically generated as
xi = Azi + i, zi ∼ N (0, Ih), i ∼ 10−3N (0, Id) (196)
h = 5 andA ∈ Rd×h is a fixed matrix that we sampled with entriesAij ∼ U([0, 1]). For the experiments showed in Figure 2, we set the step size η = 0.01 and we initialized the weight matrices asW (0)1 ,W (0) 2 ,W (0) 3 such thatW (0) ` ∼ 10−5U([0, 1]).
The FA matrices M` are generated as M`;ij ∼ U([0, 1]). We repeat the FA training for 15 times, sampling a different set of FA matrices each time (but same initial conditionsW (0)` and same matrix A): in Figure 2 we then report the average of the trace norm and the reconstruction error with error bars calculated as mean ± 2·standard_deviation (the factor 2 is to make the error bars more visible in the final plot). We complement these experiments with the corresponding GD training.

Review summaries:
This paper theoretically analyzes the Feedback Alignment (FA, Eqn 2) algorithm in the optimization of deep linear networks. * Under the assumption that the data matrices, FA matrices and initial weight matrices can be diagonalized simultaneously, the optimization can be divided into several one-dimensional problems, and the paper proves the convergence of both continuous and discrete dynamics. * After the diagonalization, each eigenvalue corresponds to a one-dimensional dynamics. The authors construct two types of initialization to show that the convergence for the larger eigenvalues can be either faster (implicit regularization) or slower (implicit anti-regularization) than the smaller eigenvalues. * Numerical comparison between FA and GD with 2- and 3-layer linear networks and random initialization. FA converges faster than GD and shares similar implicit regularization behavior.

The authors study feedback alignment (FA, Lillicrap '14), an alternative algorithm to the standard backpropagation algorithm to train neural networks. The key idea is to replace the weight matrices of the network with fixed, random "feedback matrices" when computing the weight updates while training the network. A precise understanding how such a learning rule can lead to learning remains an open question, and has recently attracted some interest (see related works)  Here, the authors study the dynamics of learning with FA in (deep) linear neural networks. The dynamics of linear neural networks trained with backpropagation has been studied extensively, thus providing an ideal test bed for this study. The authors derive a set of continuous-time equations governing the dynamics of FA for linear networks, discussing in detail their assumptions etc. They use these equations to discover an interesting implicit bias of FA: for certain initialisations, features of the data are learnt in the *inverse* order of importance, quantified by the singular value associated with each mode of the data. For other initialisations, features are learnt in the order of importance, as expected. The sensitivity of this bias to the initialisation doesn't appear in backprop. Finally, the authors also study linear auto-encoders trained with FA, and find that not FA does not only recover (a rotation of) the eigenvectors, but that it also speeds up training significantly. They also study the discrete-step dynamics of FA.

This paper studies the Feedback Alignment (FA) algorithm on deep linear networks. In continuous time, FA algorithm is an alternative to Gradient Flow (GF) In particular, to compute the time derivative of the weights, FA algorithm starts with the gradient of the weights at each layer, and replace the part that depends on the weights of succeeding layers by a fixed matrix.  The FA algorithm on deep linear networks under spectral initialization is studied regarding both convergence and implicit bias, and the analysis are provided for both continuous-time and discrete-time FA algorithm.