Learning fixed points of recurrent neural networks by reparameterizing
  the network model
In computational neuroscience, fixed points of recurrent neural networks are
commonly used to model neural responses to static or slowly changing stimuli.
These applications raise the question of how to train the weights in a
recurrent neural network to minimize a loss function evaluated on fixed points.
A natural approach is to use gradient descent on the Euclidean space of
synaptic weights. We show that this approach can lead to poor learning
performance due, in part, to singularities that arise in the loss surface. We
use a reparameterization of the recurrent network model to derive two
alternative learning rules that produces more robust learning dynamics. We show
that these learning rules can be interpreted as steepest descent and gradient
descent, respectively, under a non-Euclidean metric on the space of recurrent
weights. Our results question the common, implicit assumption that learning in
the brain should be expected to follow the negative Euclidean gradient of
synaptic weights.
Vicky Zhu, Robert Rosenbaum
# Learning fixed points of recurrent neural networks by reparameterizing the network model

###### Abstract

In computational neuroscience, fixed points of recurrent neural networks are commonly used to model neural responses to static or slowly changing stimuli. These applications raise the question of how to train the weights in a recurrent neural network to minimize a loss function evaluated on fixed points. A natural approach is to use gradient descent on the Euclidean space of synaptic weights. We show that this approach can lead to poor learning performance due, in part, to singularities that arise in the loss surface. We use a reparameterization of the recurrent network model to derive two alternative learning rules that produce more robust learning dynamics. We show that these learning rules can be interpreted as steepest descent and gradient descent, respectively, under a non-Euclidean metric on the space of recurrent weights. Our results question the common, implicit assumption that learning in the brain should be expected to follow the negative Euclidean gradient of synaptic weights.

## 1 Introduction

Recurrent neural network models (RNNs) are widely used in machine learning and in computational neuroscience. In machine learning, they are typically used to learn dynamical responses to time series inputs. In computational neuroscience, RNNs are sometimes used to model dynamical responses of neurons to dynamical stimuli [1, 2], but are also often used to model stationary, fixed point neural responses to static inputs. For example, many phenomena observed in visual cortical circuits, _e.g._, surround suppression, are widely modeled by stationary states of computational models in which recurrent connections model lateral, intralaminar connectivity [3, 4, 5, 6, 7, 8].

A natural approach to learning fixed points of RNNs is to use direct gradient descent on the recurrent weight matrix after the network has converged toward a fixed point. A direct application of this approach, called "truncated backpropagation through time," can be computationally expensive because it requires the application of backpropagation on a computational graph unrolled over many time steps. Moreover, backpropagation through time is difficult to implement or approximate with biologically plausible models of learning [9].

Alternative approached use the implicit equation for fixed points to derive the exact gradients of the loss with respect to the weight matrix at the fixed point, or some approximations to this quantity [10, 11, 12, 13, 14]. These approaches can also be computationally expensive and difficult to implement in biologically plausible models because the gradient derived from the implicit equation involves matrix inverses, which either need to be computed directly or approximated using, for example, iterative methods. In this work, we additionally show that gradient descent on the recurrent weight matrix can lead to poor learning performance because the associated loss landscape has singularities and implicit biases that make it poorly conditioned for gradient-based learning.
When mentioning "gradient descent" above, we were implicitly referring to the Euclidean gradient on weights, which is standard practice. However, several authors have argued that the default use of the Euclidean gradient in gradient descent is not necessarily optimal for studying artificial or biological learning. In machine learning applications, non-Euclidean gradients informed by information theory, such as the natural gradient, are superior in some settings [15, 16, 17]. In computational neuroscience, the use of a Euclidean gradient implicitly assumes a specific choice of units in a biological model and, more generally, assumes a specific parameterization of the model [18, 19, 20]. Different units or different parameterizations of a biological model will yield different gradients and ultimately different learning dynamics. Hence, gradient descent using the Euclidean gradient of the loss with respect to synaptic weights under a specific choice of parameterization might not capture learning dynamics or learned representations in biological neuronal networks.

In this work, we derive two new learning rules for fixed points of recurrent neural networks by reparameterizing the network model. The first learning rule can be viewed as steepest descent with respect to a non-Euclidean metric. The second rule approximates the first one, but it is more efficient and it can be interpreted as gradient descent with a non-Euclidean gradient. We demonstrate empirically that these learning rules exhibit more robust and efficient learning dynamics than standard, Euclidean gradient descent. We also find that the parameter updates produced by these rules point in substantially different directions in parameter space than the negative Euclidean gradient. In addition to providing new, robust learning rules for learning fixed points in recurrent networks, our results question the common, implicit assumption in computational neuroscience that learning should follow the negative Euclidean gradient of synaptic weights.

Code to apply the proposed learning rules and produce all figures in the manuscript can be found at

[https://github.com/RobertRosenbaum/LearningFixedPointsInRNNs](https://github.com/RobertRosenbaum/LearningFixedPointsInRNNs)

## 2 Background and theory

### Model description

We consider a recurrent neural network (RNN) model of the form [21, 22, 1, 2]

\[\tau\frac{d\mathbf{r}}{dt}=-\mathbf{r}+f(W\mathbf{r}+\mathbf{x}) \tag{1}\]

where \(\mathbf{r}(t)\in\mathbb{R}^{N}\) is a vector of model firing rates, \(\tau>0\) is a time constant, \(W\in\mathbb{R}^{N\times N}\) is a recurrent connectivity matrix, \(\mathbf{x}\in\mathbb{R}^{N}\) models external input to the network, and \(f:\mathbb{R}\rightarrow\mathbb{R}\) is a non-negative, non-decreasing activation function or "f-I curve", which is applied pointwise. For a time-constant input, \(\mathbf{x}(t)=\mathbf{x}\), fixed point firing rates satisfy

\[\mathbf{r}=f(W\mathbf{r}+\mathbf{x}). \tag{2}\]

The stability of fixed point firing rates from Eq. (1) is determined by the eigenvalues of the Jacobian matrix,

\[\mathcal{J}=\frac{1}{\tau}[-I+GW] \tag{3}\]

where \(G=\text{diag}(f^{\prime}(\mathbf{z}))\) is a diagonal matrix with entries

\[G_{jj}=f^{\prime}(\mathbf{z}_{j})\]

and \(\mathbf{z}=[W\mathbf{r}+\mathbf{x}]\) is the vector of neural inputs or pre-activations evaluated at their fixed points. Specifically, a fixed point is hyperbolically stable if all eigenvalues of \(\mathcal{J}\) have negative real part.
We can alternatively consider an recurrent neural network model in discrete time of the form

\[\mathbf{r}(n+1)=f(W\mathbf{r}(n)+\mathbf{x}(n)). \tag{4}\]

Eq. (1) is more common in computational neuroscience while Eq. (4) is more common in machine learning, but they are closely related. Eq. (4) has the same fixed points as Eq. (2), but hyperbolic stability is obtained when eigenvalues of \(GW\) have magnitude less than \(1\). Hence, if a fixed point is stable for Eq. (4), it is also stable for Eq. (1), but the converse is not true. In this work, we focus on the continuous system in Eq. (1), but our approach and learning rules can also be applied to the discrete system in Eq. (4).

In machine learning applications, RNNs are often used to learn mappings from input time series, \(\mathbf{x}(t)\), to output time series, \(\mathbf{r}(t)\), and they are often trained using backpropagation through time. In computational neuroscience, RNNs of the form in Eq. (1) are often studied for their fixed point properties, for example to study orientation selectivity and surround suppression among other phenomena [3, 4, 5, 6, 7, 8], but the weights in these studies are often chosen by hand, not learned. As a combination of these perspectives, we are interested in _learning_ mappings from static inputs, \(\mathbf{x}(t)=\mathbf{x}\), to their associated fixed points, \(\mathbf{r}\), given by Eq. (2).

Specifically, consider a supervised learning task with a cost function of the form

\[J(W)=\frac{1}{m}\sum_{i=1}^{m}L(\mathbf{r}^{i},\mathbf{y}^{i})\]

where \(\mathbf{x}^{i}\) is an input, \(\mathbf{y}^{i}\) is a label, \(L\) is a loss function, and \(\mathbf{r}^{i}=f(W\mathbf{r}^{i}+\mathbf{x}^{i})\) is the fixed point that the network converges to under input \(\mathbf{x}^{i}\). This learning task presents unique challenges because fixed points are defined _implicitly_ by Eq. (2) instead of explicitly as a function of \(\mathbf{x}^{i}\), and also because we only wish to learn _stable_ fixed points. The data set \(\{(\mathbf{x}^{i},\mathbf{y}^{i})\}_{i}\) can be the entire data set in the case of full-batch learning, or a mini-batch in the case of stochastic learning. Updates to \(W\) during learning can be written as

\[W\gets W+\Delta W\]

where

\[\Delta W=\frac{1}{m}\sum_{i=1}^{m}\Delta W^{i}\]

Here, \(\Delta W^{i}\) is an update rule that can depend on \(\mathbf{x}^{i}\), \(\mathbf{y}^{i}\), \(W^{i}\), and \(\mathbf{r}^{i}\). Below, we derive and compare three different update rules, \(\Delta W^{i}_{1}\), \(\Delta W^{i}_{2}\), and \(\Delta W^{i}_{3}\), for minimizing \(J\).

### Gradient descent on the recurrent weight matrix.

The first learning rule we consider is direct gradient descent of the loss surface with respect to \(W\) using the Euclidean gradient,

\[\Delta W^{i}_{1}=-\eta_{W}\nabla_{W}L(\mathbf{r}^{i},\mathbf{y}^{i}) \tag{5}\]

where \(\eta_{W}>0\) is a learning rate and \(\nabla_{W}\) refers to the standard, Euclidean gradient with respect to \(W\). If the fixed point, \(\mathbf{r}^{i}\), is hyperbolically stable, then the Jacobian matrix from Eq. (3) has eigenvalues with negative real part, so \(I-G^{i}W=-\tau\mathcal{J}\) is invertible and we have (see Appendix A.1)

\[\Delta W^{i}_{1}=-\eta_{W}G^{i}\left[I-G^{i}W\right]^{-T}\left(\nabla_{\mathbf{r}^ {i}}L\right)\left(\mathbf{r}^{i}\right)^{T}. \tag{6}\]

where \(G^{i}=\text{diag}(f^{\prime}(\mathbf{x}^{i}))\) evaluated at the fixed point and \(U^{-T}\) denotes the inverse transpose of a matrix, \(U\). If \(G^{i}_{jj}\neq 0\) for all \(j\), then \(G^{i}\) is invertible so Eq. (6) can be simplified to get

\[\Delta W^{i}_{1}=-\eta_{W}\left[\left[G^{i}\right]^{-1}-W\right]^{-T}\left( \nabla_{\mathbf{r}^{i}}L\right)\left(\mathbf{r}^{i}\right)^{T} \tag{7}\]
Evaluating Eqs. (6) and (7) directly is computationally expensive because they require the calculation of a matrix inverse. Truncated backpropagation through time and other methods provide alternative approaches to approximating \(\Delta W_{1}\)[10, 11, 12, 13, 14], but note that truncated backpropagation through time requires the storage of a large computational graph, making it memory inefficient. Moreover, we show in examples below that using \(\Delta W_{1}\) to update weights can lead to poor learning performance. We next propose an alternative update rule based on a nonlinear reparameterization of the model.

### A new learning rule from reparameterizing the RNN

To motivate the reparameterized model, first consider the special case of a linear network defined by

\[f(z)=z\]

In this case, \(G=I\) is the identity matrix and Eq. (2) for the fixed point can be written as

\[\mathbf{r}=[I-W]^{-1}\mathbf{x}.\]

This is a linear model in the sense that \(\mathbf{r}\) is a linear function of \(\mathbf{x}\), but the nonlinear dependence of the cost on \(W\) (especially a nonlinearity involving matrix inverses) produces complicated and computationally expensive update from Eq. (6).

Instead of performing gradient descent with respect to \(W\), we propose instead to first apply a nonlinear change of coordinates to obtain new parameters,

\[A=F(W):=[I-W]^{-1}. \tag{8}\]

If we parameterize the model in terms of \(A\) instead of \(W\), then fixed points satisfy the standard linear model

\[\mathbf{r}=A\mathbf{x} \tag{9}\]

which is linear in the input, \(\mathbf{x}\), _and_ the parameters, \(A\). Gradient descent of the loss with respect to \(A\) gives the standard update rule for a linear, single-layer neural network

\[\begin{split}\Delta A^{i}&=-\eta_{A}\nabla_{A}L( \boldsymbol{r}^{i},\boldsymbol{y}^{i})\\ &=-\eta_{A}\left(\nabla_{\boldsymbol{r}^{i}}L\right)\left( \boldsymbol{x}^{i}\right)^{T}\\ &=-\eta_{A}\left(\nabla_{\boldsymbol{r}^{i}}L\right)\left( \boldsymbol{r}^{i}\right)^{T}A^{-T}\end{split} \tag{10}\]

where we distinguish between the learning rate, \(\eta_{A}\), used for the reparameterized model and the learning rate, \(\eta_{W}\), used for the original parameterization. Eq. (10) gives a gradient-based update to the new parameter, \(A\), but our original RNN model is parameterized by \(W\). To update our original parameters, we need to change the \(\Delta A\) from Eq. (10) back to \(W\) coordinates. To do this, note that we want to find a value for \(\Delta W\) that satisfies \(A+\Delta A=F(W+\Delta W)\) whenever \(A=F(W)\) and \(\Delta A\) comes from Eq. (10). In other words, the update to \(W\) is given by

\[\begin{split}\Delta W_{2}^{i}&=F^{-1}(F(W)+\Delta A ^{i}))-W\\ &=-\left[[I-W]^{-1}-\eta_{A}\left(\nabla_{\boldsymbol{r}^{i}}L \right)\left(\boldsymbol{r}^{i}\right)^{T}[I-W]^{T}\right]^{-1}+I-W\end{split} \tag{11}\]

where \(F^{-1}(A)=I-A^{-1}\) is the inverse of \(F(W)\).

To summarize this approach, if Eq. (11) is used to update parameters, \(W\), under the linear fixed point model, \(\mathbf{r}=f(W\mathbf{r}+\mathbf{x})\) with \(f(z)=z\), then the learning dynamics will be identical to standard linear regression of parameters, \(A\), on the model \(\mathbf{r}=A\mathbf{x}\).

Since gradient descent with respect to \(A\) in Eq. (10) represents steepest descent of the loss surface in the new parameter space of \(A\) and since Eq. (11) gives the same updates in the original
parameter space of \(W\), the learning rule in Eq. (5) corresponds to steepest descent of the cost, \(J(W)\), using a non-Euclidean metric defined by

\[d(W_{1},W_{2})=\|F(W_{1})-F(W_{2})\| \tag{12}\]

where \(\|B\|=\sqrt{\text{Tr}(BB^{T})}\) is the Euclidean or Frobenius norm on matrices. Note that \(d(\cdot,\cdot)\) is a metric when restricted to the space of all matrices, \(W\), for which \(I-W\) is invertible. Hence, if we restrict to \(W\) that yield hyperbolically stable fixed points, Eq. (5) corresponds to steepest descent with respect to a non-Euclidean metric. However, the metric \(d\) is not necessarily generated by an inner product, so Eq. (11) cannot be called _gradient_ descent since the notion of a gradient requires a metric induced by an inner product. In Section 2.4, we show that an approximation to \(\Delta W_{2}^{i}\) produces gradient descent with a non-Euclidean gradient. Moreover, in Section 3, we present examples showing that \(\Delta W_{2}\) is better suited to learning fixed points than the standard approach to gradient descent represented by \(\Delta W_{1}\). But first, we need to generalize the derivation of \(\Delta W_{2}\) to arbitrary activation functions.

Eq. (11) was derived for the specific case \(f(z)=z\), but we can extend it to a model with arbitrary \(f(z)\). To do so, we first linearize Eq. (2) to obtain a linearized fixed point equation,

\[\mathbf{r}=G[W\mathbf{r}+\mathbf{x}] \tag{13}\]

which has a closed form solution given by

\[\mathbf{r}=[I-GW]^{-1}G\mathbf{x}. \tag{14}\]

Note, again, that \(I-GW\) is invertible whenever \(\mathbf{r}\) is a hyperbolically stable fixed point.

Given Eq. (14), a natural choice of new parameters would be

\[A=[I-GW]^{-1}G, \tag{15}\]

because it would again produce a (linearized) model of the form \(\mathbf{r}=A\mathbf{x}\). Note that under the linear model \(f(z)=z\), we have \(G=I\), and recover the parameterization in Eq. (8), so Eq. (15) is a generalization of Eq. (8). However, the update rule to \(W\) derived from gradient descent on \(A\) from the parameterization in Eq. (15) is susceptible to blowup or singularities when some values of \(G_{jj}=f^{\prime}(\mathbf{z}_{j})\) become small in magnitude or zero. To see why this is the case, suppose \(G_{jj}=\mathcal{O}(\epsilon)\) is small for some \(j\) and consider an update to \(W\) of the form \(W=W+\Delta W\). Then the resulting update to \(\mathbf{r}_{j}\) is, to linear order in \(\epsilon\),

\[\Delta\mathbf{r}_{j} =\sum_{k}G_{jj}\Delta W_{jk}r_{k}\] \[=\mathcal{O}(\epsilon\Delta W).\]

On the other hand, an update of the form \(A=A+\Delta A\) gives

\[\Delta\mathbf{r}_{j} =\sum_{k}\Delta A_{jk}\mathbf{r}_{k}\] \[=\mathcal{O}(\Delta A).\]

Hence, if we want \(\Delta W\) to produce the same change, \(\Delta\mathbf{r}\), produced by \(\Delta A\), then we must have \(\Delta W\sim\mathcal{O}(\Delta A/\epsilon)\). This will cause large changes to \(W\) in response to inputs for which \(G\) has small elements at the fixed point, ultimately undercutting the model's performance (see Appendix A.2 for more details). In the extreme case that \(G_{jj}=0\) for some \(j\), updates to \(W\) do not impact \(\mathbf{r}\) (_i.e._, \(\Delta\mathbf{r}_{j}=0\) for any \(\Delta W\) under the linear approximation \(\mathbf{r}=G[W\mathbf{r}+\mathbf{x}]\)), so we cannot derive a \(\Delta W\) to match a given \(\Delta A\), _i.e._, the reparameterization in Eq. (15) is ill-posed.
To circumvent these problems, we instead take the parameterization

\[A=F(W):=[G-GWG]^{-1} \tag{16}\]

in place of Eq. (15). Under the linearized fixed point equation in Eq. (13), we then obtain the linear model

\[\mathbf{r}=GAG\mathbf{x}\]

which generalizes Eq. (9). This equation is linear in \(\mathbf{x}\) and in the new parameters, \(A\). Hence, learning \(A\) is again a linear regression problem, albeit with the extra \(G\) terms. These extra \(G\) terms prevent singularities and blowup when \(G_{jj}\) terms become small or zero because \(\Delta r_{j}=\mathcal{O}(\epsilon\Delta A)\) is small whenever we make an update of the form \(A=A+\Delta A\) with \(G_{jj}=\mathcal{O}(\epsilon)\) small. Under the simple linear model \(f(z)=z\), we have \(G=I\), and recover the parameterization in Eq. (8), so that Eq. (16) (like Eq. (15)) is a generalization of Eq. (8).

Note that each input (_i.e._, each \(i\)) will potentially have a different gain matrix, \(G^{i}=\text{diag}(f^{\prime}(\mathbf{z}^{\mathbf{i}}))\), so each sample will have a potentially different value of \(A^{i}=[G^{i}-G^{i}WG^{i}]^{-1}\) as well. The gradient-based update of the loss, \(L(\mathbf{r}^{i},\mathbf{y}^{i})\), with respect to \(A^{i}\) for each sample becomes

\[\Delta A^{i} =-\eta_{A}\nabla_{A^{i}}L(\mathbf{r}^{i},\mathbf{y}^{i})\] \[=-\eta_{A}G^{i}\left(\nabla_{\mathbf{r}^{i}}L\right)(\mathbf{r}^{i})^{T}[ G^{i}]^{-1}A^{-T}\]

Using the same approach used to derive Eq. (11) above, we can again derive an update to \(W\) given by

\[\Delta W_{2}^{i} =F^{-1}(F(W)+\Delta A^{i})-W \tag{17}\] \[=-\left[\left[I-G^{i}W\right]^{-1}G^{i}-\eta_{A}\left[G^{i} \right]^{2}\left(\nabla_{\mathbf{r}^{i}}L\right)(\mathbf{r}^{i})^{T}\left[I-G^{i}W \right]^{T}G^{i}\right]^{-1}\] \[\quad+\left[[G^{i}]^{-1}-W\right].\]

This update can only be evaluated directly in the situation where \(G_{jj}^{i}\neq 0\) for all \(j\) so that the inverse of the gain matrix, \(G\), exists. However, note that \([W_{2jk}^{i}\to 0\) as \(G_{jj}^{i}\to 0\), as expected, so in situations where \(G_{jj}^{i}=0\), it is consistent to take \([W_{2}^{i}]_{jk}=0\). Note also that Eq. (17) is equivalent to Eq. (11) whenever \(G=I\), as expected, since Eq. (17) generalizes Eq. (11) to the case of arbitrary \(f\).

### A simpler learning rule from linearizing the reparameterized rule

The reparameterized rule in Eq. (17) is rather a complicated learning rule, and the matrix inverses can be computationally expensive to compute or approximate. If we assume that \(\eta_{A}>0\) is small, then we can approximate Eq. (17) by applying Taylor expansion to linear order in \(\eta_{A}\). This gives the linearized parameterized rule (see Appendix A.3 for details),

\[\Delta W_{3}^{i}=-\eta_{A}\left[I-WG^{i}\right]G^{i}\left(\nabla_{\mathbf{r}^{i}} L\right)(\mathbf{r}^{i})^{T}\left[I-G^{i}W\right]^{T}\left[I-G^{i}W\right] \tag{18}\]

In contrast to Eqs. (6) and (17) for \(\Delta W_{1}^{i}\) and \(\Delta W_{2}^{i}\), Eq. (18) for \(\Delta W_{3}^{i}\) does not require the computation of matrix inverses. Like \(\Delta W_{1}^{i}\) and \(\Delta W_{2}^{i}\), \(\Delta W_{3}^{i}\) satisfies \(\Delta W_{jk}\to 0\) whenever \(G_{jj}\to 0\), but unlike Eq. (17) for \(\Delta W_{2}^{i}\), Eq. (18) for \(\Delta W_{3}^{i}\) can be evaluated directly when \(G_{jj}=0\) for some \(j\).

Notably, \(\Delta W_{3}^{i}\) can be interpreted as gradient descent of the loss function with a non-Euclidean gradient. To see why this is the case, first note that \(\Delta W_{3}^{i}\) is related to \(\Delta W_{1}^{i}\) according to

\[\Delta W_{3}^{i}=B^{i}\Delta W_{1}^{i}C^{i}, \tag{19}\]
where

\[B^{i}=[I-WG^{i}][I-WG^{i}]^{T}\]

and

\[C^{i}=[I-G^{i}W]^{T}[I-G^{i}W].\]

Here and for the remainder of this section, we take \(\eta_{A}=\eta_{W}=\eta\) to highlight the relationship between the two update rules, but constant scalar coefficients do not affect these results.

Using Eq. (19), we may conclude that \(\Delta W^{i}_{3}\) is equivalent to gradient descent of the loss with respect to \(W\) using a non-Euclidean gradient. To explain this statement in more detail, note that the gradient of \(L(\mathbf{r}^{i},\mathbf{y}^{i})\) with respect to \(W\) depends on the choice of metric or geometry [18]. Given an inner product, \(\langle\cdot,\cdot\rangle_{a}\), on \(\mathbb{R}^{N\times N}\) the gradient of a function, \(F:\mathbb{R}^{N\times N}\rightarrow\mathbb{R}\), on the geometry imposed by \(\langle\cdot,\cdot\rangle_{a}\) evaluated at \(W\in\mathbb{R}^{N\times N}\) is the unique matrix \(\nabla^{a}_{W}F\in\mathbb{R}^{N\times N}\) such that for every \(U\in\mathbb{R}^{N\times N}\)[23, 18],

\[\langle\nabla^{a}_{W}F,U\rangle_{a}=\lim_{\epsilon\to 0}\frac{F(W+ \epsilon U)-F(W)}{\epsilon}.\]

The standard Euclidean gradient, \(\nabla=\nabla^{E}\), on matrices is given by taking the geometry produced by the Euclidean or Frobenius inner product,

\[\langle U,V\rangle_{E}=\sum_{jk}U_{jk}V_{jk}=\text{Tr}(UV^{T}).\]

Recall that \(\Delta W^{i}_{1}\) is defined by the Euclidean gradient,

\[\Delta W^{i}_{1}=-\eta\nabla^{E}_{W}L(\mathbf{r}^{i},\mathbf{y}^{i})\]

where \(L(\mathbf{r}^{1},\mathbf{y}^{i})\) is interpreted as a function of \(W\). We claim that

\[\Delta W^{i}_{3}=-\eta\nabla^{B}_{W}L(\mathbf{r}^{i},\mathbf{y}^{i}) \tag{20}\]

where \(\nabla^{B}_{W}\) is the gradient under the geometry defined by the inner product,

\[\langle U,V\rangle_{B} =\text{Tr}(B^{-1}UC^{-1}V^{T})\] \[=\langle B^{-1}U,VC^{-1}\rangle_{E}.\]

Note that we can use the cyclic property of the trace operator to write

\[\langle U,V\rangle_{B} =\text{Tr}(B^{-1}UC^{-1}V^{T})\] \[=\text{Tr}\left([I-WG]^{-T}[I-WG]^{-1}U[I-GW]^{-1}[I-GW]^{-T}V^{T}\right)\] \[=\text{Tr}\left([I-WG]^{-1}U[I-GW]^{-1}[I-GW]^{-T}V^{T}[I-WG]^{-T}\right)\] \[=\langle\mathcal{L}U,\mathcal{L}V\rangle_{E}\]

where \(\mathcal{L}:\mathbb{R}^{N\times N}\rightarrow\mathbb{R}^{N\times N}\) is a linear operator on \(N\times N\) matrices defined by

\[\mathcal{L}(U)=[I-WG]^{-1}U[I-GW]^{-1}.\]

Hence, \(\langle\cdot,\cdot\rangle_{B}\) can be viewed as a Euclidean inner product on linearly transformed coordinates. This confirms that \(\langle\cdot,\cdot\rangle_{B}\) defines an inner product on square matrices whenever \([I-WG]\) and \([I-GW]\) are non-singular. For notational convenience here and below, we do not write the explicit dependence of \(B\), \(C\), or \(\mathcal{L}\) on \(i\), but they _do_ depend on \(i\) through \(G^{i}\). In other words, there are distinct matrices, \(B\) and \(C\), and therefore distinct inner products, \(\langle\cdot,\cdot\rangle_{B}\), at each gradient descent iteration. Given Eq. (19), we can prove Eq. (20) by showing that

\[\nabla^{B}_{W}L=B\left[\nabla^{E}_{W}L\right]C. \tag{21}\]
To show Eq. (21), first define the \(N\times N\) standard basis matrices \(\mathbf{1}^{jk}\in\mathbb{R}^{N\times N}\) entrywise by

\[\mathbf{1}_{j^{\prime}k^{\prime}}^{jk}=\begin{cases}1&j=j^{\prime}\text{ and }k=k^{\prime}\\ 0&\text{otherwise}\end{cases}.\]

for \(j,k=1,\ldots,N\). Now compute the inner product of the gradient with \(\mathbf{1}^{jk}\),

\[\begin{split}\left\langle\left[\nabla^{B}L\right],\mathbf{1}^{ jk}\right\rangle_{B}&=\left\langle B^{-1}\left[\nabla^{B}L\right],\mathbf{1}^{ jk}C^{-1}\right\rangle_{E}\\ &=\text{Tr}\left(B^{-1}\left[\nabla^{B}L\right]C^{-1}\left[ \mathbf{1}^{jk}\right]^{T}\right)\\ &=\sum_{n=1}^{N}\left[B^{-1}\left[\nabla^{B}L\right]C^{-1}\mathbf{1}^ {kj}\right]_{n,n}\\ &=\sum_{n,m=1}^{N}\left[B^{-1}\left[\nabla^{B}L\right]C^{-1} \right]_{n,m}[\mathbf{1}^{kj}]_{m,n}\\ &=\left[B^{-1}\left[\nabla^{B}L\right]C^{-1}\right]_{jk}\end{split} \tag{22}\]

where the last line follows from the fact that \(\mathbf{1}_{n,m}^{kj}=1\) when \(n=k\) and \(m=j\), and it is equal to zero for all other \(j,k\). But we also have, from the definition of a gradient,

\[\left\langle\left[\nabla^{B}L\right],\mathbf{1}^{jk}\right\rangle_{B}=\lim_{ \epsilon\to 0}\frac{J(W+\epsilon\mathbf{1}^{jk})-J(W)}{\epsilon}=\frac{ \partial J}{\partial W_{jk}}=\left[\nabla^{E}L\right]_{jk}. \tag{23}\]

Since Eqs. (22) and (23) apply for all indices \(j,k\), we may conclude that

\[B^{-1}\left[\nabla^{B}L\right]C^{-1}=\left[\nabla^{E}L\right]\]

and therefore that

\[\left[\nabla^{B}L\right]=B\left[\nabla^{E}L\right]C,\]

which concludes our proof.

In summary, if \(W\) is updated according to \(\Delta W_{3}^{i}\) from Eq. (18), then this is equivalent to performing gradient descent on the loss with respect to the weight matrix under the geometry defined by the new inner product, \(\langle U,V\rangle_{B}\). Below, we present examples showing that this geometry is better suited to learning \(W\) than gradient descent with respect to the standard Euclidean geometry. Specifically, \(\Delta W_{3}^{i}\) learns more robustly than \(\Delta W_{1}^{i}\).

## 3 Experiments and results

We next evaluate and interpret each of the learning rules derived above on two different learning tasks.

### Learning fixed points in a linear model.

For demonstrative purposes, we first consider an example of linear regression with mean-squared loss. Specifically, we consider \(f(z)=z\) with

\[L(\mathbf{r},\mathbf{y})=\|\mathbf{r}-\mathbf{y}\|^{2}.\]

where \(\|\cdot\|\) is the Euclidean norm on \(\mathbb{R}^{N}\). Note that \(G=I\) is the identity in this case. We define the \(N\times m\) matrices, \(X=\begin{bmatrix}\mathbf{x}^{1}&\mathbf{x}^{2}&\ldots\mathbf{x}^{m}\end{bmatrix}\), \(Y=\begin{bmatrix}\mathbf{y}^{1}&\mathbf{y}^{2}&\ldots\mathbf{y}^{m}\end{bmatrix}\), and \(R=\begin{bmatrix}\mathbf{r}^{1}&\mathbf{r}^{2}&\ldots\mathbf{r}^{m}\end{bmatrix}= [I-W]^{-1}X\). The cost function can be written as

\[J(W)=\frac{1}{m}\left\|[I-W]^{-1}X-Y\right\|^{2}. \tag{24}\]
It is useful to also write the cost in terms of the parameters \(A=[I-W]^{-1}\) to get a standard quadratic cost function,

\[J_{A}(A)=\frac{1}{m}\left\|AX-Y\right\|^{2}. \tag{25}\]

For this problem, minimizers of \(J(W)\) and \(J_{A}(A)\) can be found explicitly. Before continuing to empirical examples, we derive and discuss these explicit minimizers.

#### 3.1.1 Computing explicit minimizers in a linear model.

In the under-parameterized case (\(N\leq m\) when all matrices full rank), \(J_{A}(A)\) has a unique minimizer defined by

\[A^{*}=YX^{+}\]

where \(X^{+}=X^{T}(XX^{T})^{-1}\) is the Moore-Penrose pseudo-inverse of \(X\) when \(N\leq m\). Therefore, \(J(W)\) has a unique minimizer at

\[W^{*}=I-[A^{*}]^{-1}=XX^{T}(YX^{T})^{-1}\]

under the assumption that \(A^{*}\) is invertible.

The over-parameterized case (\(N>m\) when matrices are full rank) is more relevant and interesting. In this case, there are infinitely many choices of \(W\) and \(A\) for which \(J(W)=0\) and \(J_{A}(A)=0\). The problem of choosing a solution to \(J_{A}(A)=0\) is a standard least squares problem and a common approach is to take

\[A^{*}=YX^{+}\]

where \(X^{+}=(X^{T}X)^{-1}X^{T}\) is the Moore-Penrose pseudo-inverse of \(X\) when \(N>m\). It is tempting to use this value of \(A^{*}\) and then take \(W^{*}=I-[A^{*}]^{-1}\). However, note that \(A^{*}\) is the solution to \(AX=Y\) that minimizes the Frobenius norm of \(A\), _i.e._,

\[A^{*}=\operatorname*{argmin}_{A}\left\|A\right\|\text{ s.t. }\ AX=Y.\]

Therefore, \(W^{*}=I-[A^{*}]^{-1}\) represents a solution, \(W\), that minimizes the norm of \(A=[I-W]^{-1}\). Since the Jacobian matrix is given by \(\mathcal{J}=(-I+W)/\tau=-A^{-1}/\tau\), stability is promoted by \(W\) having a small spectral radius (all eigenvalues of \(W\) must have real part less than \(1\) for stability). Hence, \(W^{*}=I-[A^{*}]^{-1}\) is a poor choice for \(W^{*}\). Minimizing the Frobenius norm of \(A\) will tend to push the eigenvalues of \(A\) toward zero, which can lead to large eigenvalues of \(W=I-A^{-1}\) and \(\mathcal{J}=-A^{-1}/\tau\), promoting unstable fixed points. Instead, to find a good optimizer, \(W^{*}\), we can find solutions that minimize the norm of \(W\) instead of \(A\). To this end, we can solve

\[W^{*}=\operatorname*{argmin}_{W}\left\|W\right\|\text{ s.t. }\ [I-W]^{-1}X=Y.\]

To solve this problem, we re-write it in a more standard form

\[W^{*}=\operatorname*{argmin}_{W}\left\|W\right\|\text{ s.t. }\ WY=Y-X.\]

This problem has the solution

\[W^{*}=[Y-X]Y^{+} \tag{26}\]

where \(Y^{+}=(Y^{T}Y)^{-1}Y^{T}\) is the Moore-Penrose pseudo-inverse of \(Y\) when \(N>m\). This is the solution with minimal Frobenius norm on \(W\) and is therefore more likely than \(I-[A^{*}]^{-1}\) to have a small spectral radius and therefore more likely to give stable fixed points. Hence, Eq. (26) provides a good optimizer in the over-parameterized case (\(N>m\)).
#### 3.1.2 Visualizing the loss landscape of a linear model.

For empirical examples, we first generated inputs, \(\mathbf{x}^{i}\), independently from a Gaussian distribution and generated targets \(\mathbf{y}^{i}\) using a ground truth weight matrix, \(\hat{W}\), and adding noise. Specifically, we define

\[\begin{split} X&\sim\sigma_{x}Z_{N\times m}\\ Y&\sim\left[I-\hat{W}\right]^{-1}X+\sigma_{y}Z_{N \times m}\end{split} \tag{27}\]

where \(\sigma_{x}=0.1\) controls the magnitude of the inputs, \(\sigma_{y}=0.01\), and each \(Z_{N\times m}\) represents an \(N\times m\) matrix of independent, standard, Gaussian random variables. The ground truth weight matrix is generated by

\[\hat{W}\sim\frac{\sigma_{w}}{\sqrt{N}}Z_{N\times N}.\]

Following Girko's circular law, the eigenvalues of \(\hat{W}\) lie approximately within a circle of radius \(\sigma_{w}\) with high probability [24]. Hence, we take \(\sigma_{w}=0.5<1\) to control the spectral radius of the circle to be less than \(1\), so that all eigenvalues of the Jacobian matrix, \(\mathcal{J}=(-I+\hat{W})/\tau\), are negative and fixed point firing rates are stable under the ground truth parameters, \(\hat{W}\).

The cost landscape, \(J(W)\), cannot easily be visualized as a function of \(W\) for \(N>1\) because \(W\) has \(N^{2}\) dimensions, so even \(N=2\) would be difficult to visualize. To help visualize \(J(W)\), we first plotted it on a random line segment passing through \(W^{*}\) in \(\mathbb{R}^{N\times N}\). Specifically, we defined the parameterized function

\[W(t)=W^{*}+\frac{ct}{\sqrt{N}}Z_{N\times N} \tag{28}\]

where \(c=2.5\) scales the maximum magnitude of the perturbation and \(t\) was varied from \(-1\) to \(1\) to create the visualization of \(J(W(t))\) (Figure 1A). This corresponds to plotting \(J(W)\) along a one-dimensional slice of the space \(\mathbb{R}^{N\times N}\) on which \(W\) lives. Note that the true minimizer, \(W=W^{*}\), is sampled when \(t=0\). The values of \(W\) sampled by the process can produce stable or unstable fixed points. Making the approximation \(W^{*}\approx\hat{W}\), we have that \(\rho(W(t))\approx\sqrt{\sigma_{w}^{2}+c^{2}t^{2}}\) and therefore an approximate stability condition is given by \(|t|<\sqrt{1-\sigma_{w}^{2}}/c\approx 0.346\).

Figure 1A shows the resulting cost curve for five random values of \(Z\) with the blue dashed lines marking the approximate stability boundary. The cost is relatively well behaved within the boundary, but poorly conditioned outside of the boundary because of the singularities produced by the matrix inverses in Eq. (24). Specifically, outside of the stability region, the spectral radius of \(W\) is larger than \(1\) so some eigenvalues are near \(1\) in magnitude. As a result, the \([I-W]^{-1}\) in Eq. (24) can lead to very large values of \(J(W)\).

To further visualize the loss landscape, we repeated the procedure above in two dimensions by sampling values of \(W\) from a random _plane_ passing through \(W^{*}\). Specifically, we defined the parameterized function

\[W(t_{1},t_{2})=W^{*}+\frac{c}{\sqrt{N}}(t_{1}Z_{1}+t_{2}Z_{2}) \tag{29}\]

where \(Z_{1},Z_{2}\sim N_{N\times N}(0,1)\), \(t_{1}\) and \(t_{2}\) were each varied from \(-1\) to \(1\) to create the visualization of \(J(W(t))\), and \(c=2.5\) scales the perturbation (Figure 1B). Note that \(W(t_{1},t_{2})=W^{*}\) when \(t_{1}=t_{2}=0\), so the center of the square corresponds to the minimum cost, \(J=0\). The approximate stability condition becomes \(\sqrt{t_{1}^{2}+t_{2}^{2}}<\sqrt{1-\sigma_{w}^{2}}/c=0.346\), so the approximate stability boundary is a circle (Figure 1B, dashed blue curve). Singularities create intricate ridges of large cost outside of the stability boundary (Figure 1B).

In summary, Figure 1 shows that the cost landscape, \(J(W)\), is extremely poorly conditioned outside of the stability region, _i.e._, when \(W\) has a spectral radius larger than \(1\). Note, however, that the effective cost landscape, \(J_{A}(A)\), of the reparameterized model is a simple quadratic
landscape, given by Eq. (25). Gradient-based learning according to \(\Delta W_{1}\) must traverse the poorly conditioned landscape from Figure 1. But the learning dynamics of the reparameterized rule, \(\Delta W_{2}\), are equivalent to those produced by \(A\) traversing a comparatively well-behaved quadratic landscape. We show in empirical examples below that this difference helps \(\Delta W_{2}\) and its linear approximation, \(\Delta W_{3}\), perform more robustly than \(\Delta W_{1}\).

#### 3.1.3 Gradient descent on the recurrent weight matrix in a linear model.

We first perform direct gradient descent on \(J\) with respect to \(W\) using \(\Delta W_{1}\). The gradient-based update rule from Eq. (7) can be written as

\[\begin{split}\Delta W_{1}&=\frac{1}{m}\sum_{i=1}^{ m}\Delta W_{1}^{i}\\ &=\frac{-2\eta_{W}}{m}\left[I-W^{T}\right]^{-1}\left[R-Y\right]R ^{T}.\end{split} \tag{30}\]

Empirical simulations show relatively poor learning performance (Figure 2A). Learning is slow for small learning rates, but larger learning rates fail to converge to good minima. Recall that the true minimum is zero because the model is over-parameterized. We next show that the linearized approximation to \(\Delta W_{2}\) performs similarly well.

#### 3.1.4 Learning using the reparameterized rule in a linear model.

For this linear example, the reparameterized learning rule from Eqs. (11) and (17) can be written as

\[\Delta W_{2}=\left[I-W\right]-\left(\left[I-W\right]^{-1}-\frac{2\eta_{A}}{m} (R-Y)X^{T}\right)^{-1}. \tag{31}\]

Recall that the learning dynamics produced by Eq. (31) are equivalent to those produced by learning the standard quadratic cost function, \(J_{A}(A)\), along with the standard gradient-based

Figure 1: **Visualizing the cost landscape for a linear model.****A)** The cost function \(J(W(t))\) as a function of \(t\) from Eq. (28). This represents the cost evaluated along five random line segments in \(\mathbb{R}^{N\times N}\), each passing through \(W^{*}\) at \(t=0\). Two blue dash lines show the stability boundary, \(|t|=0.346\). The vertical axis is cutoff at \(J=1000\) to better visualize the curves. Blue and black circles show stable and unstable initial conditions used for learning. **B)** The cost function \(J(W(t_{1},t_{2}))\) from Eq. (29). This represents the cost evaluated on a randomly oriented square with center at \(W^{*}\). The color axis is cutoff at \(J=1000\).
Figure 2: **Performance of three different learning rules for a linear regression problem.****A)** The cost function, \(J(W)\), from Eq. (24) for five different learning rates, \(\eta=\eta_{W}\), using direct gradient descent on the weight matrix, \(\Delta W_{1}^{i}\) from Eq. (30). **B)** Same as A, but for the reparameterized learning rule, \(\Delta W_{2}^{i}\) from Eq. (31) with \(\eta=\eta_{A}\). **C)** Same as B, but for the linearized rule, \(\Delta W_{3}^{i}\) from Eq. (33).
update rule,

\[\Delta A=-\frac{2\eta_{A}}{m}(AX-Y)X^{T}, \tag{32}\]

which is often called the "delta rule."

The behavior of the learning dynamics under Eq. (32) in the overparameterized case is well understood [25, 26, 27, 28]. Specifically, \(A\) tends toward solutions to \(AX=Y\) that minimize the distance, \(\|A-A_{0}\|\), of \(A\) from its initial condition under the Frobenius norm. As a result, \(\Delta W_{2}\) finds solutions, \(W\), to \([I-W]^{-1}X=Y\) that minimize the distance, \(d(W,W_{0})\), of \(W\) from its initial condition under the metric, \(d\), defined in Eq. (12). In addition, since the Jacobian matrix is given by \(\mathcal{J}=-A^{-1}/\tau\), we may conclude that \(\Delta W_{2}\) finds solutions that minimize the distance, \(\left\|\mathcal{J}_{0}^{-1}-\mathcal{J}^{-1}\right\|\), between the inverse Jacobian and its initial condition under the Frobenius norm.

In comparison to the gradient-based update, \(\Delta W_{1}\), from Eq. (30) (Figure 2A), we see that \(\Delta W_{2}\) from Eq. (31) performs much more robustly (Figure 2B). The cost reliably converges toward zero with increasing rates of convergence at larger learning rates.

#### 3.1.5 Learning using the linearized reparameterized rule in a linear model.

The linearized, reparameterized update rule from Eq. (18) for this linear model can be written as

\[\Delta W_{3}=-\frac{2\eta_{A}}{m}\left[I-W\right](R-Y)X^{T}\left[I-W\right]. \tag{33}\]

This approximated learning rule gives a simpler equation that is more efficient to compute, but still shows excellent agreement with the reparameterized rule, \(\Delta W_{2}\) (Figure 2C, compare to B).

Note that Eq. (33) does not require any explicit computation of matrix inverses with the exception of the computation of firing rates \(R=[I-W]^{-1}X\). However, since \(R\) is left-multiplied by \([I-W]\) in Eq. (33), we can get rid of this matrix inverse and write \(\Delta W_{3}\) in the form

\[\Delta W_{3}=-\frac{2\eta_{A}}{m}\bigg{(}XX^{T}\left[I-W\right]-\left[I-W \right]YX^{T}\left[I-W\right]\bigg{)}. \tag{34}\]

Note that the equation \(R=[I-W]^{-1}X\) and Eq. (34) are specific to the linear case, \(f(z)=z\). When using a nonlinear \(f(z)\), fixed point firing rates, \(R\), cannot generally be computed in closed form, but must be approximated by directly simulating Eq. (1) until convergence.

#### 3.1.6 Comparing the direction of updates.

To check the similarity between the updates from each learning rule, we calculated the angle between the updates at each iteration, defined by

\[\theta_{\alpha\beta}=cos^{-1}\left(\frac{\Delta W_{\alpha}\cdot\Delta W_{ \beta}}{\sqrt{(\Delta W_{\alpha}\cdot\Delta W_{\alpha})(\Delta W_{\beta} \cdot\Delta W_{\beta})}}\right)\]

for \(\alpha,\beta\in\{1,2,3\}\) where \(A\cdot B=\text{Tr}(A^{T}B)\) is the Frobenius inner product. For sufficiently small learning rates, any update, \(\Delta W\), that decreases the cost must satisfy \(\Delta W\cdot\Delta W_{1}>0\) where \(\Delta W_{1}\) is the gradient-based update [29] since the change in cost can be written as

\[\Delta J =\Delta W\cdot\nabla_{W}J+\mathcal{O}(\eta^{2}) \tag{35}\] \[=-\frac{\eta_{W}}{m}\Delta W\cdot\Delta W_{1}+\mathcal{O}(\eta^{ 2}).\]

Additionally, \(\Delta W_{2}\rightarrow\Delta W_{3}\) as \(\eta_{A}\to 0\). Hence, we should expect that \(\theta_{\alpha\beta}<90^{\circ}\) for all pairs, \(\alpha\) and \(\beta\).
Figure 3 shows the angles, \(\theta_{12}\) and \(\theta_{23}\), during learning. The angles, \(\theta_{13}\), were virtually identical to \(\theta_{23}\), so they are not shown. In each example, we used the same update, \(\Delta W_{\beta}\), to update \(W\) throughout learning. Hence, the two updates, \(dW_{\alpha}\) and \(dW_{\beta}\), were compared starting at the same initial \(W\) at each learning step.

The angle, \(\theta_{12}\), between the gradient-based updates and the reparameterized updates is relatively close to \(90^{\circ}\) (Figure 3A), indicating that they point in different directions, nearly as different as possible under the condition that they both decrease the cost. Unsurprisingly, \(\theta_{23}\) is near zero (Figure 3B), indicating that the reparameterized rule is similar to its linearization.

### Training fixed points on a nonlinear categorization task.

So far, for demonstrative purposes, we considered only simple examples of linear regression in which closed equations for optima are known. We next consider an example of image categorization using the MNIST hand-written digit benchmark.

The learning goal is to minimize a cross-entropy loss on \(C=10\) classes using one-hot encoded labels. Specifically,

\[L(\mathbf{y}^{i},\mathbf{s}^{i})=-\mathbf{y}^{i}\cdot\log(\mathbf{s}^{i})\]

where \(\mathbf{y}^{i}\) is a "one-hot" encoded label for digit \(i\),

\[\mathbf{s}^{i}_{l}=\frac{e^{\mathbf{z}^{i}_{l}}}{\sum\limits_{k=1}^{C}e^{\mathbf{z}^{i}_{ k}}},\]

is the softmax output, and \(\mathbf{z}^{i}\in\mathbb{R}^{C}\) is a logit computed from a random projection of fixed point rates of a recurrent network. Specifically,

\[\mathbf{z}^{i}=W_{\text{out}}\mathbf{r}^{i}\]

where \(W_{\text{out}}\in\mathbb{R}^{C\times N}\) is a fixed, random readout matrix and \(\mathbf{r}^{i}=f(W\mathbf{r}^{i}+\mathbf{x}^{i})\) is the fixed point from an \(N\times N\) recurrent network with input \(i\). Inputs are flattened \(28\times 28\) MNIST images, \(\mathbf{p}^{i}\in\mathbb{R}^{M}\), where \(M=28*28=784\) and we multiply them by a fixed, random read-in matrix to form the input to the network,

\[\mathbf{x}^{i}=W_{\text{in}}\mathbf{p}^{i}\]

where \(W_{\text{in}}\in\mathbb{R}^{N\times M}\) and \(N=300\) is the number of neurons in the network. We did not train \(W_{\text{out}}\) or \(W_{\text{in}}\) because we wanted to focus on the effectiveness of learning the recurrent weight

Figure 3: **Angles and correlations between weight updates.****A)** Angle (\(\theta_{12}\)) between the weight updates for the gradient-based and reparameterized learning rules. **B)** Same as A, but comparing the parameterized rule with its linearization.
matrix, \(W\). We used a hyperbolic target activation function, \(f(z)=\tanh(z)\). To compute \(\mathbf{r}\), we simulated Eq. (1) using a forward Euler scheme for \(500\) time steps with \(\tau=100dt\) where \(dt\) is the step size used in the Euler method. The model was trained on \(3\) epochs of the MNIST data set using a batch size of \(512\).

For this learning task, the gradient-based update rule from Eq. (7) can be written as

\[\Delta W_{1}=-\frac{\eta_{W}}{m}\sum_{i=1}^{m}\left[\left[G^{i}\right]^{-1}-W^ {T}\right]^{-1}W_{\text{out}}^{T}\left[\mathbf{s}^{i}-\mathbf{y}^{i}\right]\left[\mathbf{r }^{i}\right]^{T}.\]

We found that this gradient based learning rule performed poorly (Figure 4). Small learning rates learned slowly, as expected, while larger learning rates produced instabilities that caused the loss and accuracy to jump erratically during learning. Indeed, analysis of the Jacobian matrices showed that fixed points became unstable for the two largest learning rates considered in Figure 4.

We next tested the linearized, reparameterized update rule, \(\Delta W_{3}\). We did not include results for \(\Delta W^{2}\) because, as in the linear examples considered above, they are very similar to \(\Delta W_{3}\), and they are computationally more expensive to calculate. For this learning task, \(\Delta W_{3}\) can be written as

\[\Delta W_{3}=-\frac{\eta_{A}}{m}\sum_{i=1}^{m}\left[I-WG^{i}\right]G^{i}W_{ \text{out}}^{T}\left[\mathbf{s}^{i}-\mathbf{y}^{i}\right]\left[\mathbf{r}^{i}\right]^{T} \left[I-G^{i}W^{T}\right]\left[I-G^{i}W\right].\]

Using this linearized, reparameterized update rule significantly improved learning performance (Figure 5). Learning performance improved consistently with increasing learning rates and higher accuracy was achieved without instabilities. Analysis of the Jacobian matrices showed that fixed points were stable for all of the learning rates considered in Figure 5. We conclude that the

Figure 4: **Gradient based learning on a non-linear classification task. Results from training the fixed points of a recurrent network to categorize MNIST digits using the gradient-based update rule, \(\Delta W_{1}\). A,B,C,D) Training and testing losses and accuracies evaluated at each step over the course of \(3\) epochs.**
linearized, reparameterized learning rule can improve the learning of fixed points in non-linear recurrent neural network models.

## 4 Discussion

In summary, we have shown that when learning fixed points of recurrent neural network models, the direct application of gradient descent with respect to the recurrent weight matrix under the Euclidean geometry is computationally expensive and not robust. Badly conditioned loss surfaces can cause ineffective learning. Moreover, matrix inverses in the equations for the gradients are expensive to evaluate or approximate.

We derived two alternative learning rules derived from a reparameterization of the recurrent network model. These learning rules perform more robustly than the standard gradient descent approach. Moreover, one of the two learning rules is simpler and more computationally efficient. The learning rules can be interpreted as steepest descent and gradient descent on the recurrent weight matrix under a non-Euclidean metric. Our results support recent calls to re-consider the default use of Euclidean gradients on parameters in machine learning [15, 16, 17] and computational neuroscience [18, 20].

Recently, authors have argued that the use of Euclidean gradients for modeling learning in the brain is justified because any learning rule that takes small steps and reduces the loss must be positively correlated with the negative Euclidean gradient [29]. Put another way, the angle between the parameter updates and the negative Euclidean gradient must be less than \(90^{\circ}\) (see Eq. (35) and surrounding discussion). While this is true of the learning rules that we studied, the angle is very close to \(90^{\circ}\) in practice, indicating only a weak correlation. Hence, our work shows

Figure 5: **A reparameterized learning rule on a non-linear classification task.** Results from training the fixed points of a recurrent network to categorize MNIST digits using the gradient-based update rule, \(\Delta W_{3}\). **A,B,C,D)** Training and testing losses and accuracies evaluated at each step over the course of \(3\) epochs. Compare to Figure 4.
that the Euclidean gradient is not always strongly correlated with effective learning rules.

We focused on a single, fully connected recurrent layer, which limits the ease with which our model can be applied to larger data sets. Partly for this reason, we only considered the relatively simple MNIST data set as a benchmark. Future work could extend our results to multi-layer recurrent networks in which read-in and read-out matrices are trained and in which at least some fully connected layers are replaced by convolutional connectivity. These extensions will allow our approach to be applied to larger and more challenging datasets.

Fixed points of recurrent neural networks are widely used in computational neuroscience to model static neural responses to static stimuli [3, 4, 5, 6, 7, 8] and our results could be useful for these modeling approaches. On the other hand, recurrent neural networks in machine learning are almost exclusively used for time-varying inputs. Our results rely on the assumption of a time-constant input, \(\mathbf{x}(t)=\mathbf{x}\), which limits their direct application to many machine learning problems. Moreover, even in neuroscience, the assumption of a static stimulus is only an approximation. Natural stimuli are dynamical. However, if fixed points are approached faster than the stimulus changes (_i.e._, \(\tau\) is faster than \(\mathbf{x}(t)\)) then the response, \(\mathbf{r}(t)\), is approximated by the fixed point in Eq. (2) and our results provide an approximation. Moreover, a combination of our fixed point learning rules with dynamical learning rules, such as backpropagation through time, could improve learning in situations where some components of the input are static and others are dynamical. Future work should test whether our learning rules can be combined with backpropagation through time to improve performance on tasks with multiple timescales.

## 5 Acknowledgments

This material is based upon work supported by the Air Force Office of Scientific Research (AFOSR) under award number FA9550-21-1-0223 and National Science Foundation (NSF) award numbers DMS-1654268 and DBI-1707400.
Appendix

### Derivation of the Direct Gradient Descent Update, \(\Delta W_{1}\)

Here, we derive Eq. (6) for direct gradient descent on \(W\). To derive Eq. (6), it is sufficient to show that

\[\nabla_{W}L=\left(\mathbf{r}\left[\nabla_{\mathbf{r}}L\right]^{T}\left[I-GW\right]^{-1}G \right)^{T}.\]

Here we are considering a single input, label, and fixed point - \(\mathbf{x}\), \(\mathbf{y}\), and \(\mathbf{r}\) - so we can omit the \(i\) superscripts that appear in Eq. (6). Note that \(\nabla_{W}L(\mathbf{r}(W))\) is a matrix with elements

\[\frac{\partial L}{\partial W_{jk}}=\left[\nabla_{\mathbf{r}}L(\mathbf{r},\mathbf{y})\right] \cdot\frac{\partial\mathbf{r}}{\partial W_{jk}}. \tag{36}\]

To derive \(\frac{\partial\mathbf{r}}{\partial W_{jk}}\), we first derive the change of firing rate, \(\Delta\mathbf{r}\), to linear order in an update \(\Delta W_{jk}\). Consider an initial \(\mathbf{r}_{0}\) satisfying \(\mathbf{r}_{0}=f(W_{0}\mathbf{r}_{0}+\mathbf{x})\) and an update to \(W\) defined by \(W=W_{0}+\Delta W\) for some \(\Delta W\). The new fixed point satisfies \(\mathbf{r}=f(W\mathbf{r}+\mathbf{x})\) and we wish to compute \(\Delta\mathbf{r}=\mathbf{r}-\mathbf{r}_{0}\) to linear order in \(\Delta W\). Define \(\mathbf{z}_{0}=W_{0}\mathbf{r}_{0}+\mathbf{x}\) and \(\mathbf{z}=W\mathbf{r}+\mathbf{x}\). Then

\[\Delta\mathbf{r} =f(\mathbf{z})-f(\mathbf{z}_{0})\] \[=f(\mathbf{z}_{0})+f^{\prime}(\mathbf{z}_{0})(\mathbf{z}-\mathbf{z}_{0})-f(\mathbf{z }_{0})+O(\mathbf{z}-\mathbf{z}_{0})^{2}\] \[=G(\mathbf{z}-\mathbf{z}_{0})+O(\mathbf{z}-\mathbf{z}_{0})^{2}.\]

To linear order in \(\Delta\mathbf{r}\), we have

\[\Delta\mathbf{r} =G(\mathbf{z}-\mathbf{z}_{0})\] \[=G\left((W\mathbf{r}+\mathbf{x})-(W_{0}\mathbf{r}_{0}+\mathbf{x}_{0})\right)\] \[=G\left((W_{0}+\Delta W)\mathbf{r}-W_{0}\mathbf{r}_{0}\right)\] \[=G(W_{0}\mathbf{r}+\Delta W\mathbf{r}-W_{0}\mathbf{r}_{0})\] \[=G(W_{0}\Delta\mathbf{r}+\Delta W\mathbf{r})\] \[\Delta\mathbf{r}-GW_{0}\Delta\mathbf{r} =G\Delta W\mathbf{r}\] \[[I-GW_{0}]\Delta\mathbf{r} =G\Delta W\mathbf{r}.\]

As a result, we have that

\[\frac{\partial\mathbf{r}}{\partial W_{jk}}=[I-GW]^{-1}G\mathbf{1}^{jk}\mathbf{r}\]

which is interpreted as a column vector. Here, \(\mathbf{1}^{jk}\) is the matrix with all entries equal to zero except for element \((j,k)\), which is equal to \(1\). Eq. (6) then follows from the following Lemma.

**Lemma 1**.: \[[I-GW]^{-1}G\mathbf{1}^{jk}\mathbf{r}=\mathbf{r}_{k}\left[[I-GW]^{-1}G\right]_{(:,j)}\] (37)

_where \(\mathbf{r}_{k}\) is the \(k\)th element of \(\mathbf{r}\) and \(B_{(:,j)}\) denotes the \(j\)th column of a matrix, \(B\)._

Proof.: We first calculate \(\mathbf{1}^{11}\mathbf{r}\), \(\mathbf{1}^{12}\mathbf{r}\), and \(\mathbf{1}^{21}\mathbf{r}\):

\[\mathbf{1}^{11}\mathbf{r}=\begin{bmatrix}1&0&\ldots&0\\ 0&\cdot&\ldots&\cdot\\ \cdot&\cdot&\ldots&\cdot\\ 0&0&\ldots&0\end{bmatrix}\begin{bmatrix}\mathbf{r}_{1}\\ \cdot\\ \cdot\\ \cdot\\ \mathbf{r}_{M}\end{bmatrix}=\begin{bmatrix}\mathbf{r}_{1}\\ 0\\ \cdot\\ \cdot\\ 0\end{bmatrix}=\mathbf{r}_{1}I(:,1)\]
\[\mathbf{1}^{12}\mathbf{r}=\begin{bmatrix}0&1&\dots&0\\ 0&\dots&\dots&\cdot\\ \cdot&\cdot&\dots&\cdot\\ 0&0&\dots&0\end{bmatrix}\begin{bmatrix}\mathbf{r}_{1}\\ \cdot\\ \cdot\\ \cdot\\ \cdot\\ \cdot\\ \cdot\\ \cdot\end{bmatrix}=\begin{bmatrix}\mathbf{r}_{2}\\ 0\\ \cdot\\ \cdot\\ 0\end{bmatrix}=\mathbf{r}_{2}I(:,1)\] \[\mathbf{1}^{21}\mathbf{r}=\begin{bmatrix}0&0&\dots&0\\ 1&\dots&\cdot\\ \cdot&\cdot&\dots&\cdot\\ 0&\dots&\cdot\\ 0&0&\dots&0\end{bmatrix}\begin{bmatrix}\mathbf{r}_{1}\\ \cdot\\ \cdot\\ \cdot\\ \cdot\\ \cdot\\ \cdot\end{bmatrix}=\begin{bmatrix}0\\ \mathbf{r}_{1}\\ \cdot\\ \cdot\\ 0\end{bmatrix}=\mathbf{r}_{1}I(:,2).\]

Denote \(A:=[I-GW]^{-1}G\), so \(A\mathbf{1}^{11}\mathbf{r}=\mathbf{r}_{1}A_{(:,1)}\), \(A\mathbf{1}^{12}\mathbf{r}=\mathbf{r}_{2}A_{(:,1)}\), and \(A\mathbf{1}^{21}\mathbf{r}=\mathbf{r}_{1}A_{(:,2)}\). Notice that they are column vectors. WLOG, \(A\mathbf{1}^{jk}\mathbf{r}=\mathbf{r}_{k}A_{(:,j)}\)

\[LHS=\nabla_{w}L(\mathbf{r}(W))=\begin{bmatrix}\frac{dL}{dW_{1}}& \frac{dL}{dW_{12}}&\dots&\dots&\frac{dL}{dW_{M}}\\ \frac{dL}{dW_{1}}&\frac{dL}{dW_{22}}&\dots&\dots&\frac{dL}{dW_{M}}\\ \frac{dL}{dW_{1}}&\frac{dL}{dW_{22}}&\dots&\dots&\frac{dL}{dW_{M}}\\ \frac{dL}{dW_{1}}&\frac{dL}{dW_{M2}}&\dots&\dots&\frac{dL}{dW_{M}}\\ \end{bmatrix}\] \[=\begin{bmatrix}\mathbf{r}_{1}[\nabla_{\mathbf{r}}L(\mathbf{r})]\cdot A_{(:,1 )}&\mathbf{r}_{2}[\nabla_{r}L(\mathbf{r})]\cdot A_{(:,1)}&\dots&\mathbf{r}_{M}[\nabla_{\bm {r}}L(\mathbf{r})]\cdot A_{(:,1)}\\ \mathbf{r}_{1}[\nabla_{\mathbf{r}}L(\mathbf{r})]\cdot A_{(:,2)}&\mathbf{r}_{2}[\nabla_{\mathbf{r }}L(\mathbf{r})]\cdot A_{(:,2)}&\dots&\mathbf{r}_{M}[\nabla_{\mathbf{r}}L(\mathbf{r})]\cdot A_ {(:,2)}\\ \dots&\dots&\mathbf{r}_{k}[\nabla_{\mathbf{r}}L(\mathbf{r})]\cdot A_{(:,j)}&\dots\\ \dots&\dots&\dots&\mathbf{r}_{M}[\nabla_{\mathbf{r}}L(\mathbf{r})]\cdot A_{(:,M)}\end{bmatrix},\] \[RHS=\big{(}\mathbf{r}[\nabla_{\mathbf{r}}L(\mathbf{r})]^{T}[I-GW]^{-1}G \big{)}^{T}=\big{(}\mathbf{r}[\nabla_{\mathbf{r}}L(\mathbf{r})]^{T}A\big{)}^{T}\]

\[=\begin{bmatrix}\mathbf{r}_{1}[\nabla_{\mathbf{r}}L(\mathbf{r})]^{T}A_{(:,1 )}&\mathbf{r}_{1}[\nabla_{\mathbf{r}}L(\mathbf{r})]^{T}A_{(:,2)}&\dots&\mathbf{r}_{1}[\nabla_{ \mathbf{r}}L(\mathbf{r})]^{T}A_{(:,M)}\\ \mathbf{r}_{2}[\nabla_{\mathbf{r}}L(\mathbf{r})]^{T}A_{(:,1)}&\mathbf{r}_{2}[\nabla_{\mathbf{r}}L( \mathbf{r})]^{T}A_{(:,2)}&\dots&\mathbf{r}_{2}[\nabla_{\mathbf{r}}L(\mathbf{r})]^{T}A_{(:,M)} \\ \dots&\dots&\dots&\dots\\ \dots&\dots&\dots&\dots\\ \mathbf{r}_{M}[\nabla_{\mathbf{r}}L(\mathbf{r})]^{T}A_{(:,1)}&\mathbf{r}_{M}[\nabla_{\mathbf{r}} L(\mathbf{r})]^{T}A_{(:,2)}&\dots&\mathbf{r}_{M}[\nabla_{\mathbf{r}}L(\mathbf{r})]^{T}A_{(:,M)} \end{bmatrix}^{T}\] \[=\begin{bmatrix}\mathbf{r}_{1}[\nabla_{\mathbf{r}}L(\mathbf{r})]^{T}A_{(:,1 )}&\mathbf{r}_{2}[\nabla_{\mathbf{r}}L(\mathbf{r})]^{T}A_{(:,2)}&\dots&\mathbf{r}_{1}[\nabla_{ \mathbf{r}}L(\mathbf{r})]^{T}A_{(:,M)}\\ \mathbf{r}_{2}[\nabla_{\mathbf{r}}L(\mathbf{r})]^{T}A_{(:,1)}&\mathbf{r}_{2}[\nabla_{\mathbf{r}}L( \mathbf{r})]^{T}A_{(:,2)}&\dots&\mathbf{r}_{2}[\nabla_{\mathbf{r}}L(\mathbf{r})]^{T}A_{(:,M)} \\ \dots&\dots&\dots&\dots\\ \mathbf{r}_{M}[\nabla_{\mathbf{r}}L(\mathbf{r})]^{T}A_{(:,1)}&\mathbf{r}_{M}[\nabla_{\mathbf{r}} L(\mathbf{r})]^{T}A_{(:,2)}&\dots&\mathbf{r}_{M}[\nabla_{\mathbf{r}}L(\mathbf{r})]^{T}A_{(:,M)} \end{bmatrix}^{T}\] \[=\begin{bmatrix}\mathbf{r}_{1}[\nabla_{\mathbf{r}}L(\mathbf{r})]\cdot A_{(:, 1)}&\mathbf{r}_{2}[\nabla_{\mathbf{r}}L(\mathbf{r})]\cdot A_{(:,1)}&\dots&\mathbf{r}_{M}[ \nabla_{\mathbf{r}}L(\mathbf{r})]\cdot A_{(:,1)}\\ \mathbf{r}_{1}[\nabla_{\mathbf{r}}L(\mathbf{r})]\cdot A_{(:,2)}&\mathbf{r}_{2}[\nabla_{\mathbf{r }}L(\mathbf{r})]\cdot A_{(:,2)}&\dots&\mathbf{r}_{M}[\nabla_{\mathbf{r}}L(\mathbf{r})]\cdot A _{(:,2)}\\ \dots&\dots&\mathbf{r}_{k}[\nabla_{\mathbf{r}}L(\mathbf{r})]\cdot A_{(:,j)}&\dots\\ \mathbf{r}_{1}[\nabla_{\mathbf{r}}L(\mathbf{r})]\cdot A_{(:,M)}&\mathbf{r}_{2}[\nabla_{\mathbf{r }}L(\mathbf{r})]\cdot A_{(:,M)}&\dots&\dots&\dots\\ \mathbf{r}_{M}[\nabla_{\mathbf{r}}L(\mathbf{r})]\cdot A_{(:,M)}\end{bmatrix}\] \[=LHS.\]

Combining Eq. (36) with Eq. (37) gives

\[\nabla_{W}L=\left(\mathbf{r}\left[\nabla_{\mathbf{r}}L(\mathbf{r})\right]^{T}\left[I-GW \right]^{-1}G\right)^{T}\]

which can be simplified to get Eq. (6) for \(\Delta W_{1}\).
### Analysis of a natural reparameterization and its linear approximation

We now consider the updates given by the reparameterization \(A=\left[G^{-1}-W\right]^{-1}\). The direct reparameterized update, \(\Delta W_{2}\) in this case is given by

\[\Delta W_{2} =-\left[\left(A-\eta_{A}(\nabla_{\boldsymbol{r}}L)(\boldsymbol{x} )^{T}\right)^{-1}-A^{-1}\right]\] \[=-\left[\left(\left[G^{-1}-W\right]^{-1}-\eta_{A}(\nabla_{ \boldsymbol{r}}L)(\boldsymbol{r})^{T}\left[G^{-1}-W^{T}\right]\right)^{-1}-G^{ -1}-W\right].\]

Proof.: Since \(A=\left[G^{-1}-W\right]^{-1}\), we have \(W=G^{-1}-A^{-1}\). Let \(W^{0}\) and \(A^{0}\) represent previous step update before \(W\) and \(A\), then

\[\Delta W =W-W^{0}\] \[=G^{-1}-A^{-1}-\left(\left[G^{0}\right]^{-1}-[A^{0}]^{-1}\right)\] \[=\left(G^{-1}-\left[G^{0}\right]^{-1}\right)-A^{-1}+[A^{0}]^{-1}\] \[=-\left(A^{0}+\Delta A\right)^{-1}+[A^{0}]^{-1}\] \[=-\left(A^{0}-\eta_{A}\left(\nabla_{\boldsymbol{r}}L\right)( \boldsymbol{x})^{T}\right)^{-1}+[A^{0}]^{-1}\] \[=-\left(A^{0}-\eta_{A}\left(\nabla_{\boldsymbol{r}}L\right)( \boldsymbol{r})^{T}\left[A^{0}\right]^{-T}\right)^{-1}+[A^{0}]^{-1}.\]

To get the expression that has only \(G\) and \(W\), we can substitute \(A=[G^{-1}-W]^{-1}\) and \(A^{-1}=G^{-1}-W\), and use \(G=G^{T}\) and \(G^{-T}=G^{-1}\) since \(G\) is a diagonal matrix. This gives

\[\Delta W_{2} =-\left(A-\eta_{A}\left(\nabla_{\boldsymbol{r}}L\right)( \boldsymbol{r})^{T}\,A^{-T}\right)^{-1}+A^{-1}\] \[=-\left(\left[G^{-1}-W\right]^{-1}-\eta_{A}\left(\nabla_{ \boldsymbol{r}}L\right)(\boldsymbol{r})^{T}\left[G^{-1}-W^{T}\right]\right)^{- 1}+\left[G^{-1}-W\right].\]

Note that as \(G_{jj}\to 0\), \(A_{jj}^{-1}=[G_{jj}^{i}]^{-1}-W_{jj}\to\infty\), so this reparameterizatin is poorly behaved in situations where \(G_{jj}^{j}=f^{\prime}(\boldsymbol{z}_{j})\) becomes small or zero because the second term in the sum diverges while the first term does not.

We also show that linearizing this parameterization around \(\eta_{A}=0\) still leads to updates that diverge when elements of \(G\) become small. Following the linearization from Section 2.4, the linearized, reparameterized update is given by

\[\Delta W_{3} =-\eta_{A}A^{-1}(\nabla_{\boldsymbol{r}}L)(\boldsymbol{x})^{T}A^{ -1}\] \[=-\eta_{A}\left[G^{-1}-W\right](\nabla_{\boldsymbol{r}}L)( \boldsymbol{r})^{T}\left[G^{-1}-W^{T}\right]\left[G^{-1}-W\right].\]

Proof.: First note that \(\left.\Delta W_{2}\right|_{\eta_{A}=0}=0\), so we have to linear order in \(\eta_{A}\),

\[\Delta W_{2}=\left.\frac{d\Delta W_{2}}{d\eta_{A}}\right|_{\eta_{A}=0}\eta_{A }+\mathcal{O}(\eta_{A}^{2}) \tag{38}\]

Now let

\[V=A+\Delta A=A-\eta_{A}(\nabla_{\boldsymbol{r}}L)\left(A^{-1}\boldsymbol{r} \right)^{T}\]

then \(\Delta W_{2}=A^{-1}-V^{-1}\) so

\[\frac{d\Delta W_{2}}{d\eta_{A}} =\frac{dA^{-1}}{d\eta_{A}}-\frac{dV^{-1}}{d\eta_{A}}\] \[=V^{-1}\frac{dV}{d\eta_{A}}V^{-1}\]
since \(dA^{-1}/d\eta_{A}=0\). Combining this with Eq. (38) and the definition of \(V\) gives the linearized update,

\[\Delta W_{3} =\left.V^{-1}\frac{dV}{d\eta_{A}}V^{-1}\right|_{\eta_{A}=0}\eta_{A}\] \[=V^{-1}\left(-\left(\nabla_{\boldsymbol{r}}L\right)\left(A^{-1} \boldsymbol{r}\right)^{T}\right)V^{-1}\Bigg{|}_{\eta_{A}=0}\eta_{A}\] \[=-A^{-1}(\nabla_{\boldsymbol{r}}L)\left(\boldsymbol{r}\right)^{T }A^{-T}A^{-1}\eta_{A}\] \[=-\left[G^{-1}-W\right]^{-1}\left(\nabla_{\boldsymbol{r}}L \right)\left(\boldsymbol{r}\right)^{T}\left[G^{-1}-W^{T}\right]\left[G^{-1}-W \right]\eta_{A}.\]

Again, substitute \(A^{-1}=G^{-1}-W\), to get the final expression. Notice that \(\Delta W_{3}=A^{-1}A^{-T}\Delta W_{1}A^{-T}A^{-1}\) so one can let \(B=A^{-1}A^{-T}\) and \(C=A^{-T}A^{-1}\), which are symmetric, and \(\Delta W_{3}=B\Delta W_{1}C\). 

Note, again, that \(\Delta W_{3}\) diverges if elements of \(G\) go to zero. Therefore, the natural reparameterization \(A=[G^{-1}-W]^{-1}\) is not well suited for learning.

### Linearization of the corrected reparameterization

Here, we derive the linearized update, \(\Delta W_{3}\), given in Eq. (18). This update rule is derived by expanding \(\Delta W_{2}\) from Eq. (17) to linear order. Recall that \(\Delta W_{2}\) was derived from the reparameterization \(A=[G-GWG]^{-1}\). Let \(U=[G^{-1}-W]^{-1}\), then we can rewrite Eq. (17) as

\[\Delta W_{2}=-\left[U-\eta_{A}G^{2}\left(\nabla_{\boldsymbol{r}}L\right) \left(\boldsymbol{r}\right)^{T}\left[I-GW\right]^{T}G\right]^{-1}+U^{-1}.\]

Now, denote everything inside of the inverse as \(V\) so

\[V=U-\eta_{A}G^{2}\left(\nabla_{\boldsymbol{r}}L\right)\left(\boldsymbol{r} \right)^{T}\left[I-GW\right]^{T}G.\]

Then Eq. 17 can be further rewritten as

\[\Delta W_{2}=U^{-1}-V^{-1}.\]

Now, following the same approach as Appendix A.2, note that \(\left.\Delta W_{2}\right|_{\eta_{A}=0}=0\), so the linearization of \(\Delta W_{2}\) around \(\eta_{A}=0\) is given by

\[\Delta W_{3} =\left.\frac{d\Delta W_{2}}{d\eta_{A}}\right|_{\eta_{A}=0}\eta_{A}\] \[=\left[\frac{dU^{-1}}{d\eta_{A}}-\left(-V^{-1}\frac{dV}{d\eta_{A} }V^{-1}\right)\right]_{\eta_{A}=0}\eta_{A}\] \[=V^{-1}\left(-G^{2}\left(\nabla_{\boldsymbol{r}}L\right)\left( \boldsymbol{r}\right)^{T}\left[I-GW\right]^{T}G\right)V^{-1}\Bigg{|}_{\eta_{A} =0}\eta_{A}\] \[=U^{-1}\left(-G^{2}\left(\nabla_{\boldsymbol{r}}L\right)\left( \boldsymbol{r}\right)^{T}\left[I-GW\right]^{T}G\right)U^{-1}\eta_{A}\] \[=-\eta_{A}[G^{-1}-W]G^{2}\left(\nabla_{\boldsymbol{r}}L\right) \left(\boldsymbol{r}\right)^{T}\left[I-GW\right]^{T}G[G^{-1}-W]\] \[=-\eta_{A}\left[I-WG\right]G\left(\nabla_{\boldsymbol{r}}L\right) \left(\boldsymbol{r}\right)^{T}\left[I-GW\right]^{T}[I-GW].\]