{'fraction_non_alphanumeric': 0.04269341904483621, 'fraction_numerical': 0.01480471413265803, 'mean_word_length': 4.690375023092555, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Biological credit assignment through dynamic inversion of feedforward networks

Learning depends on changes in synaptic connections deep inside the brain. In multilayer networks, these changes are triggered by error signals fed back from the output, generally through a stepwise inversion of the feedforward processing steps. The gold standard for this process -- backpropagation -- works well in artificial neural networks, but is biologically implausible. Several recent proposals have emerged to address this problem, but many of these biologically-plausible schemes are based on learning an independent set of feedback connections. This complicates the assignment of errors to each synapse by making it dependent upon a second learning problem, and by fitting inversions rather than guaranteeing them. Here, we show that feedforward network transformations can be effectively inverted through dynamics. We derive this dynamic inversion from the perspective of feedback control, where the forward transformation is reused and dynamically interacts with fixed or random feedback to propagate error signals during the backward pass. Importantly, this scheme does not rely upon a second learning problem for feedback because accurate inversion is guaranteed through the network dynamics. We map these dynamics onto generic feedforward networks, and show that the resulting algorithm performs well on several supervised and unsupervised datasets. We also link this dynamic inversion to Gauss-Newton optimization, suggesting a biologically-plausible approximation to second-order learning. Overall, our work introduces an alternative perspective on credit assignment in the brain, and proposes a special role for temporal dynamics and feedback control during learning.


Introduction
Synaptic credit assignment refers to the difficult task of relating a motor or behavioral output of the brain to the many neurons and synapses that produced it (Roelfsema and Holtmaat, 2018) -a problem which must be solved in order for effective learning to occur. While credit is assigned in artificial neural networks through the backpropagation of error gradients (Rumelhart et al., 1986;LeCun et al., 2015), a direct mapping of this algorithm to biology leads to several characteristics that are either in conflict with what is currently known about neural circuits, or that violate harder physical constraints, such as the local nature of synaptic plasticity (Grossberg, 1987;Crick, 1989).
Many biologically-plausible modifications to backpropagation have been proposed over the years (Whittington and Bogacz, 2019), with several recent studies focusing on one issue in particular, the fact that error is fed back using an exact copy of the forward weights (the "weight transport" or "weight symmetry" problem, Lillicrap et al. (2020)). Recently, it was discovered that random feedback weights are sufficient to train deep networks on modest supervised learning problems (Lillicrap et al., 2016). However, this method appears to have shortcomings in scaled-up tasks, as well as in convolutional and bottleneck architectures (Bartunov et al., 2018;Moskovitz et al., 2018). Several studies have therefore aimed to identify the necessary precision of feedback (Nøkland, 2016;Xiao et al., 2018), and others have proposed to learn separate feedback connections (Bengio, 2014;Lee et al., 2015;Akrout et al., 2019;Lansdell et al., 2019). While it is indeed plausible that feedback weights are updated alongside forward ones, these schemes complicate credit assignment by making error backpropagation dependent upon an additional learning problem (with uncertain accuracy), and by potentially introducing more learning phases.
One important characteristic of biological neural circuits is their dynamic nature, which has been harnessed in many previous learning models (Hinton et al., 1995;O'Reilly, 1996;Rao and Ballard, 1999). Here, we take inspiration from this dynamical perspective, and propose a model of error backpropagation as a feedback control problem -during the backward pass, feedback connections are used in concert with forward connections to dynamically invert the forward transformation, thereby enabling errors to flow backward. Importantly, this inversion works with arbitrary fixed feedback weights, and avoids introducing a second learning problem for the feedback. In the following, we derive this dynamic inversion, map it onto deep feedforward networks, and demonstrate its performance on several supervised tasks, as well as an autoencoder task. We also make a link between the algorithm and second-order learning. Finally, we discuss the biological implications of this perspective, and the relation to previous dynamic algorithms for credit assignment.

Notation and forward transformation
We consider nonlinear feedforward networks with L layers. The forward pass (forward transformation; Fig. 1a) from one layer to the next is where h l ∈ R d l is the activity of layer l, g(·) is an arbitrary element-wise nonlinearity, a l ∈ R d l is the "pre-activation" activity of layer l, and W l ∈ R d l ×d l−1 denotes the weight matrix from layers l − 1 to l (including bias). The input data, network output, and supervised target are denoted x = h 0 , y = h L , and t, respectively. The error is denoted e = y − t, and the loss function is L(x, t).

Error backpropagation and inversion of the forward transformation
For such networks, each layer's weights are commonly optimized using gradient descent: with the backpropagated error δ l = ∂L/∂a l ∈ R d l . We write δ l in a generalized recursive form where • is the Hadamard (element-wise) product, D l = diag(g (a l )), and M l = W T l ∈ R d l−1 ×d l is the feedback weight matrix (the source of the weight transport problem).
As mentioned above, learning can sometimes be achieved with a fixed random feedback matrix, a strategy termed feedback alignment (FA), in part due to the observed alignment between the forward weights and the pseudoinverse of the feedback weights during training (Lillicrap et al., 2016). The authors of this study also describe a biologically-implausible idealization of this algorithm, pseudobackprop (PBP), which propagates errors through the pseudoinverse of the current feedforward weights. These results, as well as other studies proposing to learn feedback as an inverted forward transformation (e.g., target prop, Bengio (2014); Lee et al. (2015)), motivate the perspective that the goal of credit assignment is to invert the feedforward transformation of the network.
We summarize these variants of backpropagation as different choices for M l in Eq.
where B l is a fixed random matrix and W + l is the Moore-Penrose pseudoinverse of W l . Backpropagation formulated as an integral control problem -error between the forward transformation (δ l ) and the target error value (δ l ) is integrated and fed back to produce a new target error δ l−1 . c: Dynamic inversion during the backward pass implements this control problem.

Dynamic inversion as feedback control
We now introduce a simple recurrent architecture (Fig. 1b,c) which dynamically and implicitly performs inversions similar to those explicitly performed by pseudobackprop and target prop as outlined above. Considering the backward pass of a linear feedforward network (Eq. (1), with g(x) = x), the error from the l-th layer, δ l , should be transformed into an error for the (l − 1)-th layer, δ l−1 . From a linear feedback control perspective, we can let the l-th layer feed a control signal, u(t) ∈ R d l , into the (l − 1)-th layer, such that the state of this layer, δ l−1 , when propagated through the feedforward transformation of the network, reproduces, as close as possible, the target error vector, δ l , of layer l. We define this as a linear control problem of the following form: where δ l−1 (t) ∈ R d l−1 is the system state of layer l − 1,δ l (t) ∈ R d l is the readout or forward transformation of this system, u l (t) ∈ R d l is the control signal fed back from layer l, and B l ∈ R d l−1 ×d l is a matrix of arbitrary feedback weights. We define a fixed, target error value for the readout, δ l , and a separate controller error, e l (t) =δ l (t) − δ l .

Leaky integral control
A standard approach in designing a controller is to use a proportional-integral-derivative (PID) formulation (Åström and Murray, 2010) that acts on the controller error e l (t), with dynamicṡ where K p , K i , and K d are coefficient matrices for the proportional, integral, and derivative components, respectively, along with an additional leak component with coefficients K u . For mathematical simplicity and biological plausibility, we only consider the integral and leak components (see Discussion for interpretation of other terms), setting their coefficients to K i = I l , and K u = −αI l . These components have a direct interpretation in rate networks (Dayan and Abbott, 2001), and have been used in other neuroscience and biological contexts (Miller and Wang, 2006;Somvanshi et al., 2015). We thus obtain the following leaky integral-only controlleṙ which acts on Eq. (5). For a fixed target δ l , this controller has the steady-state equality which suggests that the steady state of δ l−1 approximates the target δ l through the forward transformation (for small α). For α > 0, δ l−1 can be written as (using Eq. (5) in the steady-state), When α = 0, only one of these equalities will hold, depending on the dimensionalities d l and d l−1 . For expository purposes, we also write the solution as a function of the control signal u l : where and W + l is the Moore-Penrose pseudoinverse of the forward matrix W l . We thus see that this system dynamically inverts the forward transformation of the network (for small α), implicitly solving the linear system W l δ l−1 = δ l . For d l ≥ d l−1 (expanding layer), W l has a well-defined left pseudoinverse (or inverse, for d l = d l−1 ), and so the inversion follows directly from Eq. (9). In contrast, for d l < d l−1 (contracting layer), the system may have infinite solutions. The dynamics instead solves the fully-determined system (W l B l − αI)u l = δ l , which is then projected through B l to obtain δ l−1 (i.e., one solution to the desired linear system, constrained by B l ).

Linear stability and and initialization
The dynamic inversion will only be useful if it is stable and fast. Integral-only control may exhibit substantial transient oscillations, which can be mitigated if the system dynamics are fast compared to the controller. Assuming this separation of timescales, we can study the controller dynamics from Eq. (8) when the system is at its steady state (δ l−1 = B l u l from Eq. (5)): Linear stability thus depends on the eigenvalues of (W l B l −αI). Generally, the stability of interacting neural populations (and the eigenvalues of arbitrary matrix products), is an open question and we do not aim to solve it here. We instead propose that clever initialization of B l will provide stability throughout training (in addition to a non-zero leak, α). One easy way to ensure this is to initialize B l = −W T l (0), which makes the matrix product negative semi-definite (zero index indicates the start of training). From Eq. (12), this also means that for d l < d l−1 , dynamic inversion will use the Moore-Penrose pseudoinverse at the start of training. Note that this initialization does not suggest a correspondence between the forward and backward weights throughout training, as they may become unaligned when the forward weights are updated. In the case where d l > d l−1 , the matrix product is singular and requires α > 0 (but see Supplementary Materials for an alternative architecture).

Nonlinearities
Returning to the general case of nonlinear feedforward networks, the controller in this case iṡ where g(·) is an arbitrary nonlinearity such as tanh or ReLU (we keep the backward pass linear for simplicity). Note that we have opted to derive the network from the perspective of the pre-activation variables for consistency with backpropagation. Now, the steady-state for the system is Deriving explicit relationships between δ l−1 and δ l , is trickier, especially with common transfer functions like tanh and ReLU, which do not have well-defined inverses (at least numerically). To gain an intuition, we use somewhat sloppy notation and write an implicit, non-unique inverse g −1 (·), for which g −1 (g(δ l−1 )) ≈ δ l−1 and g(g −1 (δ l )) ≈ δ l . We can then write δ l−1 recursively as with M l from Eq. (12), suggesting that the network is able to approximately invert nonlinear transformations as well. We note that stability is no longer guaranteed, but in practice we find that linear stability still provides a decent indication of stability in the general case. a: Sequential backward pass (right), in which the error is inverted through one layer at a time, with each layer first receiving control signals from the layer above, and then acting as the controller for the layer below. B: Repeat layer backward pass, enabling each layer to both give and receive control, so that the full backward pass converges at once. C: Single loop backward pass features feedback from the output layer to the first hidden layer, which effectively inverts each hidden layer.

Architectures for chained dynamic inversion
The most direct way of mapping multiple dynamic inversions onto a feedforward network is to prescribe that each inversion happens sequentially -from the output to the first hidden layer -with only one pair of layers dynamically interacting at a time (sequential backward pass, Fig. 2a). Such a scheme begins by feeding the output error, δ L , into the output layer, which provides control to the last hidden layer until convergence to the target δ L−1 . Next, this target is held fixed and is re-passed as input back into layer L − 1, which now acts as a controller for layer L − 2, to obtain the target δ L−2 . This is repeated until the first hidden layer converges to its target, δ 1 . This scheme requires a backward pass with multiple steps for networks with more than one hidden layer (L − 1 steps).
The fact that each hidden layer functions as both a recipient of control, and a controller itself, motivates the second architecture, in which the hidden layers have two separate populations, each serving one of these roles (repeat layer backward pass, Fig. 2b). For the forward pass to remain unchanged, these populations (h A l and h B l ) have an identity transformation between them. During the backward pass, each controller receives the target value δ l as it settles, speeding up convergence. The steady state errors will be the same as in the sequential case, but with a single backward pass.
An alternative approach in chaining multiple dynamic inversion control problems together is to turn them into a single problem (single loop backward pass, Fig. 2c). In this scheme, the output layer acts as the controller for the activity of the first hidden layer. This architecture thus uses a single dynamic inversion, with the forward transformation encompassing nearly the entire forward pass. Interestingly, for a linear network, such controls leads each subsequent hidden layer to converge to the same error value as the previous two architectures. However, in general, the forward pass may be highly nonlinear, requiring a similarly complex control signal for successful inversion, and so this architecture could converge to a different solution compared to the other two architectures.

Update rules
We define the backpropagated error signal δ l for dynamic inversion (DI) as the steady state of the feedback control dynamics from Eq. (16), and weight updates as in Eq. (2). Biases are not included in the dynamics of the backward pass, but are updated with the converged error signals similar to standard backprop. Unlike backprop, however, the derivative of the nonlinear transfer function does not need to be applied, as it is implicitly incorporated in the inversion. As a point of comparison, we also implement a non-dynamic inversion (NDI), in which the explicit linear inversion from Eq. (10) is combined with nonlinearity gradients following Eqs. (2) and (3). c: Stability of DI as measured by the maximum real eigenvalues of the system dynamics. d: Angle between the backpropagated error vectors to the hidden layer (δ 1 ) for NDI vs DI. e: Angle between feedback weights and negative transpose of forward weights (shown for DI and NDI). f: Angle between backpropagated error vectors for DI and PBP. g-j: Same as a,c,d, e but for nonlinear regression. f: Angle between the backpropagated error vectors to the hidden layer for NDI vs SDI. Dashed lines refer to output layer, and solid linear refer to middle layer in h-k.

Relation to second-order learning
The inversion of the forward weights in DI suggests a resemblance to second-order learning, which we illustrate here with a specific example. We consider the update rules for the penultimate weight matrix, W L−1 , of a network performing regression with squared-error loss. For simplicity in the comparison, we also assume that the input to each layer is whitened. Following the formulation of Gauss-Newton (GN) optimization for deep networks in Botev et al. (2017) (Supplementary Materials), we can derive the block-diagonal GN update for W L−1 as where D L−1 = diag(g (a L−1 )). As can be seen, such a scheme features an inversion of the forward pass similar to Eq. (16), but with an explicit first-order approximation of the the nonlinearity. Dynamic inversion instead handles nonlinearities implicitly in the dynamics. However, these should be similar for small error values, suggesting that dynamic inversion approximates Gauss-Newton optimization.

Experiments
We tested dynamic inversion (DI) and non-dynamic inversion (NDI) against backpropagation (BP), feedback alignment (FA), and pseudobackprop (PBP) on four modest supervised and unsupervised learning tasks -linear regression, nonlinear regression, MNIST classification, and MNIST autoencoding. We varied the leak values (α) and feedback initializations ("Tr-Init", B l = −W T l ; "R-Init", random stable B l ) for DI and NDI. To impose stability for random initialization, we optimized the feedback matrix B l using smoothed spectral abscissa (SSA) optimization (Vanbiervliet et al., 2009;Hennequin et al., 2014)   . c,d: Examples of input and output digits and a 2D t-SNE representation of the 30-dim. latent space for BP with orthogonal weight initialization and NDI for random weight initialization.

Linear and nonlinear function approximation
Following Lillicrap et al. (2016), we tested the algorithms on a simple linear regression task with a two-layer network (dim. 30-20-10). Training was done with a fixed learning rate (Fig. 3a), or with fixed norm weight updates (Fig. 3b) in order to probe the update directions that each algorithm finds. All algorithms were able to solve this simple task with ease, with DI, NDI, and PBP converging faster than BP and FA in the fixed norm case, providing further evidence for the relation to secondorder learning. With the exception of substantial variability in the first~100 iterations, the dynamic inversion remained stable throughout training for all examples shown, with updates well-aligned to the non-dynamic version (Fig. 3c,d). Furthermore, the alignment between the feedback and the negative transpose of the forward weights settled to around 45 degrees for all DI and NDI models, which also produced alignment with the PBP updates (Fig. 3e,f).
Next, we tested performance on a small nonlinear regression task with a three-layer network (dim. 30-20-10-10, tanh nonlinearities) also following Lillicrap et al. (2016). All inversion algorithms finished with better performance compared to BP and FA (Fig. 3g), but often with slower convergence, which was unexpected considering the approximation to second-order optimization. We speculate that this may be due to inaccuracies in approximating the loss curvature (or overly small, "conservative" updates), which is common for many such algorithms (Martens, 2014). DI dynamics remained stable throughout training, and closely followed NDI updates (except for α = 0, likely due to slow convergence; Fig. 3h,i). Alignment of B varied with the layer and algorithm (Fig. 3j) -some layers settled to~45 degrees, but others remained close to zero -this is intriguing, but may be due to the simplicity of the problem. Finally, we also simulated single-loop dynamic inversion (SDI), which behaved very differently to DI, suggesting that it converges to different steady states (Fig. 3k).

MNIST classification and autoencoder
We next tested dynamic inversion on the MNIST handwritten digit dataset, where we use the standard training and test datasets (LeCun et al., 1998), with a two-layer architecture (dim. 784-1000-10 as in Lillicrap et al. (2016)). Due to the close alignment of the updates for DI and NDI (with the exception of α = 0), and the computational complexity of simulating DI for large-scale problems, we chose to only simulate NDI and use it as a proxy for DI performance. All algorithms performed similarly (Fig.  4a), with PBP and NDI having a slightly worse test error (BP, 1.9%; FA, 1.7%; PBP, all NDI, 2.7%), possibly due to slower convergence again.
Finally, we trained an autoencoder network on the MNIST dataset (dim. 784-500-250-30-250-500-784; nonlinearities tanh-tanh-linear-tanh-tanh-linear, a reduced version of Hinton and Salakhutdinov (2006)) with mini-batch training (100 examples per batch). Notably, first-order optimization algorithms have trouble dealing with the "pathological curvature" of such problems and often have very slow learning (Martens, 2010) (especially FA, Lansdell et al. (2019)). We trained BP, FA, PBP and NDI with two weight initializations -random Gaussian, and random orthogonal, which has been shown to speed up learning (Saxe et al., 2013). We found that BP only learns successfully with orthogonal weight initialization, whereas PBP and NDI perform decently in either case, further suggesting they use second-order information. Notably, PBP and NDI performance is slower with orthogonal initialization, where second-order information is not useful (but this might be mitigated by having non-orthogonal feedback weights). Furthermore, FA performed poorly in both cases, providing evidence that dynamic inversion is superior to random feedback. BP with orthogonal initialization and NDI with random initialization result in similar performance (Fig. 4c,d).

Discussion
Several previous studies have proposed biologically-plausible solutions to credit assignment involving dynamics of some kind, but most contain conceptual differences with our framework -e.g., learning using differences in activity over time or phase (contrastive Hebbian learning, O'Reilly (1996); Scellier and Bengio (2017)), using explicit error-encoding neurons (Whittington and Bogacz, 2017), or incorporating dendritic compartments (Guerguiev et al., 2017;Sacramento et al., 2018;Payeur et al., 2020). Furthermore, the aim of most of these models is to approximate error gradients, whereas dynamic inversion is unique in its link to second-order optimization. Most relevant to our work is a recent paper which also proposes to re-use forward weights in order to propagate errors through a feedback loop (Kohan et al., 2018). However, the authors do not formulate this as an inversion of the forward pass, and they use a contrastive learning scheme that differs substantially from what we do here. Finally, target propagation also aims to learn (non-dynamic) inversions (Bengio, 2014;Lee et al., 2015), making it conceptually similar -in principle dynamic inversion can also be used to propagate targets, which may afford some biological plausibility or other benefits.

Biological implications
Unlike many other biologically-plausible algorithms for credit assignment, dynamic inversion does not require precise feedback weights. This is a crucial distinction, as it not only relaxes the assumptions on feedback wiring, but also allows for feedback to be used concurrently for other roles, such as attention and prediction (Gilbert and Li, 2013). The architectures we proposed for chained dynamic inversion (Fig. 2) suggest different ways of using feedback for learning, and even leaves the possibility for direct feedback to much lower areas which is known to exist in the brain (Felleman and Van Essen, 1991). Additionally, our work depends upon the stability and control of recurrent dynamics between interacting populations (or brain areas), which has received recent interest in neuroscience (Joglekar et al., 2018). Stability and fast convergence of dynamic inversion requires slow control (Eq. (13)), suggesting that higher-order areas should be slower than the lower areas they control. Indeed, activity is known to slow down as it moves up the processing hierarchy (Murray et al., 2014).

Limitations and future work
We see two main limitations of dynamic inversion as a model for credit assignment in the brain. First, its dynamic nature means that the backward pass requires time to converge. Fast convergence would be easier to achieve with additional proportional and derivative control, as in Eq. (7) (Åström and Murray, 2010). Spike-based representations could help here since they effectively add a derivative component to the signal (Eliasmith and Anderson, 2004;Boerlin et al., 2013;Abbott et al., 2016).
Second, due to the relationship of our scheme with second-order learning and natural gradient methods, dynamic inversion could share some of the same problems (Martens, 2014;Kunstner et al., 2019). While our method avoids the cost of explicitly calculating and inverting a Hessian or Gauss-Newton matrix, in common with standard second-order methods (Pearlmutter, 1994;Schraudolph, 2002;Martens, 2010), its stability will be enhanced if the recurrent dynamics can be designed to better condition the inverse computations.
The true test of dynamic inversion will be whether or not it can be successfully scaled up to larger tasks (Bartunov et al., 2018;Xiao et al., 2018). Even so, it may be useful in other contexts where it is necessary to invert a computation, such as motor control (Wolpert and Kawato, 1998) and sensory perception (Pizlo, 2001). As an example, it was pointed out in a recent paper (Vértes and Sahani, 2019) that the successor representation -used in reinforcement learning and requiring an inverse to calculate explicitly -can be achieved dynamically in a similar way to what we propose here.
7 Supplementary materials 7.1 Alternative control architecture for expanding layers As noted in Section 3.2, the eigenvalues of the matrix (W l B l − αI) provide a good measure of the linear stability of dynamic inversion (true linear stability is measured by the eigenvalues of the block matrix in Eq. (S9) below). This precludes stability for α = 0 and d l > d l−1 (expanding layer), as the matrix product W l B l will be singular. To address this, we propose an alternative control architecture: where now the target error for layer l − 1, δ l−1 , integrates the error between δ l andδ l directly, through the feedback matrix B l . This can be interpreted as proportional feedback control with a fast controller. In this system, stability instead depends on the matrix B l W l − αI (assuming the readout dynamics are fast, and sõ δ l = W l δ l−1 ). Note that this scheme either requires identical feedback weights for the target error δ l and the current estimateδ l , or a separate population which calculates the error between these, propagated back as B l (δ l − δ l ) = B l e l .

Relation of dynamic inversion to Gauss-Newton optimization
Following the derivations in Botev et al. (2017), we can write the block-diagonal sample Gauss-Newton (GN) matrix for a particular layer l as where Q l = h l−1 h T l−1 is the sample input covariance to layer l and G l is the "pre-activation" GN matrix, defined recursively as with D l = diag(g (a l )), and C l is a square-root representation of G l . The GN update to the weight matrix of layer l can be written in vectorized form as ∆vec(W l ) ∝ −G −1 l g, where g is a vectorized version of the standard backprop gradient, as in Eq. (2). In order to avoid vectorization (and thus simplify the comparison with dynamic inversion), we make the assumption that the input to this layer is whitened, making Q l = I l . This allows us to write the GN update in non-vectorized form: Note that the approximate equality is due to the assumption that both (C T l+1 W l+1 D l ) has a left pseudoinverse, and C l+1 has a right pseudoinverse, which depends on the relative dimensionality of layers l and l + 1. Considering the simplest case, optimizing the penultimate set of weights WL−1 for a network solving regression with squared error loss, we have GL = I (and thus CL = I), and δL = e, and so the update becomes ∆WL−1 ∝ −(WLDL−1) + eh T L−1 , which is the same as Eq. (17). As mentioned in the main text, this bears a resemblance to dynamic inversion, which in this case would use the update A more thorough analysis is merited on the relationship between Eqs. (S7) and (S8) (as well as the types of nonlinear inversions found in (S8) and a more general comparison of dynamic inversion and GN optimization) but is out of the scope of this study. We note that layer-wise whitening is performed in a recent model proposing to map natural gradient learning onto feedforward networks (Desjardins et al., 2015), suggesting that the strategic placement of whitening transformations in a network with dynamic inversion may produce a more accurate approximation to Gauss-Newton or natural gradient optimization, while still being biologically plausible.

Stability optimization
In general, the dynamic inversion system dynamics for a particular layer l are not stable when initialized with random matrices W l and B l (R-Init). We follow procedures outlined in Vanbiervliet et al. (2009) and Hennequin et al. (2014) to optimize linear stability by minimizing the smoothed spectral abscissa (SSA; a smooth relaxation of the spectral abscissa, the maximum real eigenvalue). The full system matrix can be written in block form as with the first and second rows corresponding to the dynamics of δ l−1 and u l , respectively, of Eqs. (5) and (10). In brief, we calculate the gradient of the SSA with respect to the matrix B l , and make small steps until the maximum eigenvalue is sufficiently negative. We refer the reader to the references above for details. SSA optimization can be done both on W l and B l , but we chose to optimize only B l in order to have full control on the initialization of W l .