{'fraction_non_alphanumeric': 0.04016309887869521, 'fraction_numerical': 0.01653641408993091, 'mean_word_length': 4.402104747919726, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Target Propagation via Regularized Inversion

Target Propagation (TP) algorithms compute targets instead of gradients along neural networks and propagate them backward in a way that is similar yet different than gradient back-propagation (BP). The idea was first presented as a perturbative alternative to back-propagation that may achieve greater accuracy in gradient evaluation when training multi-layer neural networks (LeCun et al., 1989). However, TP has remained more of a template algorithm with many variations than a well-identified algorithm. Revisiting insights of LeCun et al., (1989) and more recently of Lee et al. (2015), we present a simple version of target propagation based on regularized inversion of network layers, easily implementable in a differentiable programming framework. We compare its computational complexity to the one of BP and delineate the regimes in which TP can be attractive compared to BP. We show how our TP can be used to train recurrent neural networks with long sequences on various sequence modeling problems. The experimental results underscore the importance of regularization in TP in practice.


Introduction
Target propagation algorithms can be seen as perturbative learning alternatives to the gradient back-propagation algorithm, where virtual targets are propagated backward instead of gradients (Le Cun, 1986;Le Cun et al., 1988;Rohwer, 1989;Mirowski and LeCun, 2009;Bengio, 2014;Goodfellow et al., 2016). A high-level summary is presented in Fig. 1: while gradient back-propagation considers storing intermediate gradients in a forward pass, target propagation algorithms proceed by computing and storing approximate inverses. The approximate inverses are then passed on backward along the graph of computations to finally yield a weight update for stochastic learning.
Target propagation aims to take advantage of the availability of approximate inverses to compute better descent directions for the objective at hand. Bengio et al. (2013); Bengio (2020) argued that the approach could be relevant for problems involving multiple compositions such as the training of Recurrent Neural Networks (RNNs), which generally suffer from the phenomenon of exploding or vanishing gradients (Hochreiter, 1998;Bengio et al., 1994;Schmidhuber, 1992). Recently, empirical results indeed showed the potential advantages of target propagation over classical gradient back-propagation for training RNNs on several tasks (Manchev and Spratling, 2020). However, these recent investigations remain built on multiple approximations, which hinder the analysis of the core idea of TP, i.e., using layer inverses.
On the theoretical side, difference target propagation, a modern variant of target propagation, was related to an approximate Gauss-Newton method, suggesting interesting venues to explain the benefits of target propagation (Bengio, 2020;Meulemans et al., 2020). Previous works have considered approximating inverses by adding multiple reverse layers (Manchev and Spratling, 2020;Meulemans et al., 2020;Bengio, 2020). However, it is unclear whether such reverse layers actually learn layer inverses during the training process. Even if they were, the additional cost of computational complexity of learning approximate inverses should be carefully accounted for.
In this work, we propose a simple target propagation approach, revisiting the original insights of Le Cun et al. (1988) on the critical importance of the good conditioning of layer inverses. We define regularized inverses through a variational formulation and we obtain approximate inverses via these regularized inverses. In this spirit, we can also interpret the difference target propagation formula (Lee et al., 2015) as a finite difference approximation of a linearized regularized inverse. We propose a smoother formula that can directly be integrated into a differentiable programming framework.
We detail the computational complexity of the proposed target propagation and compare it to the one of gradient back-propagation, showing that the additional cost of computing inverses can be effectively amortized for very long sequences. Following the benchmark of Manchev and Spratling (2020), we observe that the proposed target propagation can perform better than classical gradient-based methods on several tasks involving RNNs.
The code to reproduce the experiments is provided at https://github.com/vroulet/tpri. The appendix details the implementations of the algorithms and discuss previous interpretations of target propagation as a Gauss-Newton method.
Related work. Many variations of back-propagation algorithms have been explored; see Werbos (1994); Goodfellow et al. (2016) for an extensive bibliography. Closer to target propagation, penalized formulations of the training problem have been considered to decouple the optimization of the weights in a distributed way or using an ADMM approach (Carreira-Perpinan and Wang, 2014;Taylor et al., 2016;Gotmare et al., 2018). Rather than modifying the backward operations in the layers, one can also modify the weight updates for deep forward networks by using a regularized inverse (Frerix et al., 2018). Wiseman et al. (2017) recast target propagation as an ADMM-like algorithm for language modeling and reported disappointing experimental results. Recently, in a careful experimental benchmark evaluation, Manchev and Spratling (2020) explored further target propagation to train RNNs, mapping a sequence to a single final output, in an attempt to understand the benefits of target propagation to capture long-range dependencies, and obtained promising experimental results. Another line of research has considered synthetic gradients that approximate gradients using an additional layer instead of using back-propagated gradients (Jaderberg et al., 2017;Czarnecki et al., 2017) to speed up the training of deep neural networks. Recently, Ahmad et al. (2020); Dalm et al. (2021) considered using analytical inverses to implement target propagation and blend it with what they called a gradientadjusted incremental formula. Yet, an additional orthogonality penalty is critical for their approach to work. Recently, Meulemans et al. (2020) considered using as many reverse layers as forwarding operations. We focus here on the optimization gains of using target propagation that cannot be obtained by adding a prohibitive number of reverse layers. Finally, we do not discuss the biological plausibility of TP since we are unable to comment on this. We refer the interested reader to, e.g., (Bengio, 2020).

Target Propagation with Linearized Regularized Inverses
While target propagation was initially developed for multi-layer neural networks, we focus on its implementation for recurrent neural networks, as we shall follow the benchmark of Manchev and Spratling (2020) in the experiments. Recurrent Neural Networks (RNNs) are also a canonical family of neural networks in which interesting phenomena arise in back-propagation algorithms.
Problem setting. A simple RNN parameterized by θ = (W hh , W xh , b h , W hy , b y ) maps a sequence of inputs x 1:τ = (x 1 , . . . , x τ ) to an outputŷ = g θ (x 1:τ ) by computing hidden states h t ∈ R p corresponding to the inputs x t .
Formally, the outputŷ and the hidden states h t are computed as an output operation following transition operations defined asŷ where s is, e.g., the soft-max function for classification tasks, a is a non-linear operation such as the hyperbolic tangent function, and the initial hidden state is generally fixed as h 0 = 0. Given samples of sequence-output pairs (x 1:τ , y), the RNN is trained to minimize the error (y, g θ (x 1:τ )) of predictingŷ = g θ (x 1:τ ) instead of y.
As one considers longer sequences, RNNs face the challenge of exploding/vanishing gradients ∂g θ (x 1:τ )/∂h t (Bengio and Frasconi, 1995); see Appendix A for more discussion. We acknowledge that specific parameterization-based strategies have been proposed to address this issue of exploding/vanishing gradients, such as orthonormal parameterizations of the weights (Arjovsky et al., 2016;Helfrich et al., 2018;Lezcano-Casado and Martınez-Rubio, 2019). The focus here is to simplify and understand target propagation as a backpropagation-type algorithm using RNNs as a workbench. Indeed, training RNNs is an optimization problem involving multiple compositions for which approximate inverses can easily be available. The framework could also be potentially applied to, e.g., time-series or control models (Roulet et al., 2019).
Given the parameters W hh , W xh , b h of the transition operations, we can get approximate inverses of f θ,t (h t−1 ) for all t ∈ {1, . . . , τ }, that yield optimization surrogates that can be better performing than the ones corresponding to regular gradients. We present below a simple version of target propagation based on regularized inverses and inverse linearizations.
Back-propagating targets. The idea of target propagation is to compute virtual targets v t for each layer t = τ, . . . , 1 such that if the layers were able to match their corresponding target at time t, i.e., f θ,t (h t−1 ) ≈ v t , the objective would decrease. The final target v τ is computed as a gradient step on the loss w.r.t. h τ . The targets are then back-propagated using an approximate inverse 1 f −1 θ,t of f θ,t at each time step. Formally, consider an RNN that computed τ states h 1 , . . . , h τ from a sequence x 1 , . . . , x τ with associated output y. For a given stepsize γ h > 0, we propose to back-propagate targets by computing The update rule (2) blends two ideas: i) regularized inversion; ii) linear approximation. We shall describe below that our update (2) allows us to interpret the "magic formula" of difference target propagation in Eq. 15 of Lee et al. (2015) as 0th-order finite difference approximation, while ours is a 1st-order linear approximation. We shall also show that (2) puts in practice an insight from Bengio (2020) suggesting to use the inverse of the gradients in the spirit of a Gauss-Newton method. Once all targets are computed, the parameters of the transition operations are updated such that the outputs of f θ,t at each time step move closer to the given target. Formally, the update consists of a gradient step with stepsize γ θ on the squared error between the targets and the current outputs, i.e., for As for the parameters θ y = (W hy , b y ) of the output operation, they are updated by a simple gradient step on the loss with a stepsize γ θ .

Regularized Inversion
To explore further the original idea of Le Cun et al. (1988), we consider using the variational definition of the inverse, As long as v t belongs to the image f θ,t (R p ) of f θ,t , this definition recovers exactly the inverse of v t by f θ,t . More generally, if v t ∈ f θ,t (R p ), Eq. (4) computes the best approximation of the inverse in the sense of the Euclidean projection. When one considers an activation function a and θ h = (W hh , W xh , b h ), the solution of (4) can easily be computed.
Formally, for the sigmoid, the hyperbolic tangent or the ReLU, their inverse can be obtained analytically for any v t ∈ a(R p ). So for v t ∈ a(R p ) and W hh full rank, we get If v t ∈ a(R p ), the minimizer of (4) is obtained by first projecting v t onto a(R p ), before inverting the linear operation. To account for non-invertible matrices W hh , we also add a regularization in the computation of the inverse. Overall we consider approximating the inverse of the layer by a regularized inverse of the form , with r > 0 and π a projection onto a(R p ).
Regularized inversion vs. parameterized inversion. Bengio (2014); Manchev and Spratling (2020) parameterize the inverse as a reverse layer such that , and learn the parameters θ = (V, c) for this reverse layer to approximate the inverse of the forward computations. The parameterized layer needs to be learned to get a good approximation which involves numerically solving an optimization problem for each layer. These optimization problems come with a computational cost that can be better controlled by using regularized inversions presented earlier.
However, the approach based on parameterized inverses may lack theoretical grounding, as pointed out by Bengio (2020), as we do not know how close the learned inverse is to the actual inverse throughout the training process. In contrast, the regularized inversion (4) is less ad hoc and clearly defined and, as we shall show in the experiments, leads to competitive performance on real datasets.
In any case, the analytic formulation of the inverse gives simple insights on an approach with parameterized inverses. Namely, the analytical formula suggests parameterizing the reverse layer s.t. (i) the reverse activation is defined as the inverse of the activation and not any activation, (ii) the layer uses a non-linear operation followed by a linear one instead of the usual scheme, i.e., a linear operation followed by a non-linear one.

Linearized Inversion
Earlier instances of target propagation used direct inverses of the network layers such that the target propagation update formula would read v t−1 = f −1 θ,t (v t ) in (2). Yet, we are unaware of a successful implementation of TP using directly the inverses. To circumvent this issue, Lee et al. (2015) proposed the difference target propagation formula that back-propagates the targets as . If the inverses were exact, the difference target propagation formula would reduce to v t−1 = f −1 θ,t (v t ). Lee et al. (2015) introduced the difference target propagation formula to mitigate the approximation error of the inverses by parameterized layers. The difference target propagation formula can naturally be interpreted as an approximation of the linearization used in (2), as where ∂ h f −1 θ,t (h t ) denotes the Jacobian of the inverse of the layer at h t .  We show in Appendix D that the first-order approximation we propose (2) leads to slightly better training curves than the finite-difference approximation. Moreover, our interpretation illuminates the "mystery" of this formula, which appeared to be critical to the success of target propagation.

Gradient Back-propagation versus Target Propagation
Graph of computations. Gradient back-propagation and target propagation both compute a descent direction for the objective at hand. The difference lies in the oracles computed and stored in the forward pass, while the graph of computations remains the same. To clarify this view, we reformulate target propagation in terms of displacements λ t := v t − h t such that Eq. (1), (2) and (3) read Target propagation amounts then to computing a descent direction d θ h for the parameters θ h with a graph of computations, illustrated in Fig. 2, analogous to that of gradient-back-propagation illustrated in Appendix A. The difference lies in the use of the Jacobian of the inverse The implementation of TP with the formula (2) can be done in a differentiable programming framework, where, rather than computing the gradient of the layer, one evaluates the inverse and keep the Jacobian of the inverse. With the precise graph of computation of TP and BP, we can compare their computational complexity explicitly and bound the difference of the directions they output.
Arithmetic complexity. Clearly, the space complexities of gradient back-propagation (BP) and our implementation of target propagation (TP) are the same since the Jacobians of the inverse, and the original gradients have the same size. In terms of time complexity, TP appears at first glance to introduce an important overhead since it requires the computation of some inverses. However, a close inspection of the formula of the regularized inverse reveals that a matrix inversion needs to be computed only once for all time steps. Therefore the cost of the inversion may be amortized if the length of the sequence is particularly long.
Formally, the time complexity of the forward-backward pass of gradient back-propagation is essentially driven by matrix-vector products, i.e., where d is the dimension of the input x t , q is the dimension of the parameters θ h , for a function f we denote by T (f ) the time complexity to evaluate f and we consider e.g.
On the other hand, the time complexity of target propagation is is the cost of encoding the inverse, which, in our case, amounts to the cost of encoding g θ : z → (W hh W hh + r I) −1 W hh , such that our regularized inverse can be computed as Encoding g comes at the cost of inverting one matrix of size p. Therefore, the time-complexity of target propagation can be estimated as i.e., for long sequences whose length is larger than the dimension of the hidden states, the cost of TP with regularized inverses is approximately the same as the cost of BP. If a parameterized inverse was used rather than a regularized inverse, the cost of encoding the inverse would correspond to the cost of updating the reverse layers by, e.g., a stochastic gradient descent. This update has a cost similar to BP. However, it is unclear whether these updates get us close to the actual inverses.
Bounding the difference between target propagation and gradient back-propagation. As the computational graphs of BP and TP are the same, we can bound the difference between the oracles returned by both methods. First, note that the updates of the parameters of the output functions are the same since, in TP, gradients steps of the loss are used to update these parameters. The difference between TP and BP lies in the updates with respect to the parameters of the transition operations. For BP, the updates are computed by chain rule as where the term ∂h τ /∂h t decomposes along the time steps as ∂h τ /∂h t = τ s=t+1 ∂ h f θ,s (h s−1 ). The direction computed by TP has the same structure, namely it can be decomposed for γ h = 1 as We can then bound the difference between the directions given by BP or TP as, for any matrix norm · as formally stated in the following lemma.
Lemma 3.1. The difference between the oracle returned by gradient back-propagation ∂ θ h (y, g θ (x 1:τ )) and the oracle returned by target propagation can be bounded as For regularized inverses, we have, For the two oracles to be close, we then need the preactivation u t = W xh x t + W hh h t−1 + b h to lie in the region of the activation function that is close to being linear s.t. ∇a(u t ) ≈ I. We also need (W hh W hh + r I) −1 to be close to the identity which can be the case if, e.g., r = 0 and the weight matrices W h were orthonormal. By initializing the weight matrices as orthonormal matrices, the differences between the two oracles can be closer. However, in the long term, target propagation appears to give better oracles, as shown in the experiments below.
Target propagation as a Gauss-Newton method? Recently target propagation has been interpreted as an approximate Gauss-Newton method, by considering that the difference target propagation formula approximates the linearization of the inverse, which itself is a priori equal to the inverse of the gradients (Bengio, 2020;Meulemans et al., 2020;2021). Namely, provided that f −1 By composing the inverses of the gradients, we get an update similar to the one of Gauss-Newton (GN) method. Namely, recall that if n invertible functions f 1 , . . . , f n were composed to solve a least square problem of the form f n • . . .•f 1 (x)−y 2 2 , a Gauss-Newton update would take the form In other words, GN and TP share the idea of composing the inverse of gradients. However, numerous differences remain as detailed in Appendix C. Note that even if TP was approximating GN, it is unclear whether GN updates are adapted to stochastic problems. In any case, by using an analytical formula for the inverse, we can test this interpretation by using non-regularized inverses, which would amount to directly use the inverses as in a GN method. If the success of TP could be explained by its interpretation as a GN method, we should observe efficient training curves when no regularization is added.

Experiments
In the following, we compare our simple target propagation approach, which we shall refer to as TP, to gradient Back-Propagation referred to as BP. We follow the experimental benchmark of Manchev and Spratling (2020) to which we add results on RNNs on CIFAR and GRUs on FashionMNIST. Additional experimental, details on the initialization and the hyper-parameter selection can be found in Appendix D.
Data. We consider two synthetic datasets generated to present training difficulties for RNNs and several real datasets consisting of scanning images pixel by pixel to classify them (Hochreiter and Schmidhuber, 1997;Le et al., 2015;Manchev and Spratling, 2020). Temporal order problem. A sequence of length T is generated using a set of randomly chosen symbols {a, b, c, d}. Two additional symbols X and Y are added at positions t 1 ∈ [T /10, 2T /10] and t 2 ∈ [4T /10, 5T /10]. The network must predict the correct order of appearance of X and Y out of four possible choices {XX, XY, Y X, Y Y }.
Adding problem. The input consists of two sequences: one is made of randomly chosen numbers from [0, 1], and the other one is a binary sequence full of zeros except at positions t 1 ∈ [1, T /10] and t 2 ∈ [T /10, T /2]. The second position acts as a marker for the time steps t 1 and t 2 . The goal of the network is to output the mean of the two random numbers of the first sequence (X t1 + X t2 )/2. Image classification pixel by pixel. The inputs are images of (i) grayscale handwritten digits given in the database MNIST (LeCun and Cortes, 1998), (ii) colored objects from the database CIFAR10 (Krizhevsky, 2009) or (iii) grayscale images of clothes from the database FashionMNIST (Xiao et al., 2017). The images are scanned pixel by pixel and channel by channel for CIFAR10, and fed to a sequential network such as a simple RNN or a GRU network (Cho et al., 2014). The inputs are then sequences of 28 × 28 = 784 pixels for MNIST or FashionMNIST and 32 × 32 × 3 = 3072 pixels for CIFAR with a very long-range dependency problem. We also consider permuting the images of MNIST by a fixed permutation before feeding them into the network, which gives potentially longer dependencies in the sequential data.
Model. In both synthetic settings, we consider randomly generated mini-batches of size 20, a simple RNN with hidden states of dimension 100, and hyperbolic tangent activation. For the temporal order problem, the last layer uses a soft-max function on top of a linear operation, and the loss is the cross-entropy. For the adding problem, the last layer is linear, the loss is the mean-squared error, and a sample is considered to be accurately predicted if the mean squared error is less than 0.04 as done by (Manchev and Spratling, 2020). For the classification of images with sequential networks, we consider mini-batches of size 16 and a cross-entropy loss. For MNIST and CIFAR, we consider a simple RNN with hidden states of dimension 100, hyperbolic tangent activation, and a softmax output. For FashionMNIST, we consider a GRU network and adapted our implementation of target propagation to that case while using hidden states of dimension 100 and a softmax output.
Target propagation can tackle long sequences better than gradient back-propagation. In Fig. 3, we observe that TP performs better than BP on the temporal ordering problem: it is able to reach 100% accuracy in fewer iterations than BP for sequences of length 60 and, for sequences of length 120, it is still able to reach 100% accuracy in fewer than 40 000 iterations while BP is not. On the other hand, for the adding problem, TP performs less well than BP. The contrast in performance between the two synthetic tasks was also observed by (Manchev and Spratling, 2020) using difference target propagation with parameterized inverses. The main difference between these tasks is the different nature of the outputs, which are binary for the temporal problem and continuous for the adding problem.
In Fig. 4, we observe that TP generally performs better than BP for image classification tasks. For the MNIST dataset, it reaches around 74% accuracy after 4 · 10 4 iterations. This phenomenon is also observed with permuted images, where the optimization appears smoother, and TP obtains around 86% accuracy after 4 · 10 4 iterations and is still faster than BP. On the CIFAR dataset, no algorithms appear to reach a significant accuracy, though TP is still faster. On the FashionMNIST dataset, where a GRU network is used, our implementation of TP performs on par with BP, which shows that our approach can be generalized to more complex networks than a simple RNN.
Target propagation requires a non-zero regularization term. As mentioned in Sec. 3, by using an analytical formula to compute the inverse of the layers, we can question the interpretation of TP as a Gauss-Newton method, which would amount to TP without regularization. To understand the effect of the regularization term, we computed the area under the training loss curve of TP for 400 iterations on a log 10 grid of varying step-sizes γ θ and regularizations r for a fixed γ h = 10 −3 . The results are presented in Fig. 6a, where the smaller the area, the brighter the point and the absence of dots in the grid mean that the algorithm diverged. Fig. 6a shows that without regularization we were not able to obtain convergence of the algorithm. Simply using the gradients of the inverse as in a Gauss-Newton method may not directly work for RNNs. Additional modifications of the method could be added to make target propagation closer to Gauss-Newton, such as inverting the layers with respect to their parameters as proposed by Bengio (2020). For now, the regularization appears to successfully handle the rationale of target propagation.
Target propagation is adapted for long sequences. In Fig. 6b, we compare the performance of BP and TP in terms of accuracy after 400 iterations on the MNIST problem for various widths determined by the size of the hidden states and various lengths determined by the size of the inputs (i.e., we feed the RNN with k pixels at a time, which gives a length 784/k). Fig 6b shows that TP is generally appropriate for long sequences, while BP remains more efficient for short sequences. TP can then be seen as an interesting alternative for dynamical problems which involve many discretization steps as in RNNs and related architectures.

Conclusion
We proposed a simple target propagation approach grounded in two important computational components, regularized inversion, and linearized propagation. The proposed approach also sheds light on previous insights and successful rules for target propagation. The code is available to facilitate the reproduction of the results. We have used target propagation within a stochastic gradient outer loop to train neural networks for a fair comparison to stochastic gradient using gradient backpropagation. Developing adaptive stochastic gradient algorithms in the spirit of Adam that lead to boosts in performance when using target propagation instead of gradient backpropagation is an interesting avenue for future work. Continuous counterparts of target propagation in a neural ODE spirit is also an interesting avenue for future work. Given differentiable activation functions a, the training of recurrent neural networks is amenable to optimization by gradient descent. The gradients can be computed by gradient back-propagation implemented in modern differentiable programming software (Rumelhart et al., 1986;Werbos, 1994;Paszke et al., 2019). The gradient back-propagation algorithm is illustrated in Fig. 6. Formally, the gradients are computed by the chain rule such that, for a sample (y, x 1:τ ) and θ h = (W hh , W xh , b h ), The term ∂h τ /∂h t decomposes along the time steps as As τ grows, the norm of the term ∂h τ /∂h t can then either increase to infinity (exploding gradients) or exponentially decrease to 0 (vanishing gradients). This phenomenon may prevent the RNN from learning from dependencies between temporally distant events (Hochreiter, 1998). Several solutions were proposed to tackle this issue, including changing the network architecture (Hochreiter and Schmidhuber, 1997), Hessian-free optimization (Sutskever et al., 2011), gradient clipping and regularization (Pascanu et al., 2012), or orthonormal parametrizations (Arjovsky et al., 2016;Helfrich et al., 2018;Lezcano-Casado and Martınez-Rubio, 2019). We consider here propagating targets instead of gradients as first presented by LeCun and co-workers (Le Cun, 1986;Le Cun et al., 1988) and recently revisited by Bengio and co-workers (Bengio, 2014;Lee et al., 2015).

B.1 Target Propagation for RNNs
As detailed in Sec. 3, target propagation with linearized regularized inverses amounts to move along a descent direction computed by a forward-backward algorithm akin to gradient propagation. The iterations of linearized target propagation are then summarized in Algo. 1. The iterations of Algo. 1 make calls to any algorithm providing a descent direction which is computed by Algo. 2. In the implementation of the regularized inverses, since the inverse of activation functions such as the sigmoid or the tangent hyperbolic is numerically unstable, we consider projecting on a subset of a(R p ). For the hyperbolic tangent, we clip the target to [−1 + ε, 1 − ε] for ε = 10 −3 . Concretely, for an hyperbolic tangent activation function, the projection is then π(x) = (min(max( To read Algo. 2, we recall our notations for θ = (W hh , W xh , b h , W hy , b y ): Note that Algo. 2 can also be used for mini-batches of sequence-output pairs since all operations are either elementwise or linear with respect to the sample of sequence-output pair.   5) and (6), stepsize γ θ , total number of iterations K 2: for k = 1 . . . K do 3: Draw a sample or a mini-batch of sequences-output pairs (x 1:τ , y).

B.2.1 Formulation
Starting from h 0 = 0, given an input sequence x 1 , . . . , x τ , the GRU network (as implemented in Pytorch 2 (Paszke et al., 2019)), iterates for t = 1, . . . , τ , where is the Hadamard product, σ is a sigmoid. In the following, we will denote simply θ = (θ m , θ z , θ n ) the parameters of the network with The output of the network is e.g. a soft-max operation on the hidden state computed at the last step (if applied to an image scanned pixel by pixel for example). See the main paper for the expression of the output in that case.

B.2.2 Modifying the chain rule
The underlying idea of our implementation of target propagation in a differentiable programming framework is to mix classical gradients and Jacobians of the inverse of the functions. Denote for a given output loss L computed on a given mini-batch with the current parameters θ,∂L/∂h t the direction back-propagated by our implementation of target propagation until the step h t . The directions for the parameters of the network can be output aŝ The main task is to define∂L/∂h t−1 given∂L/∂h t and appropriate regularized inverses. For that, we start with the chain rule for ∂h t /∂h t−1 and we will replace some of the gradients by Jacobians of regularized inverses at some places.
Classical chain rule. We have Now for ∂n t /∂h t−1 , we further decompose the function f n,t (h t−1 ) as with ∇g t (u) = diag(tanh (W in x t + b in + u)).
Inverses. Now, the variables z t , m t and a t are functions of h t that incorporate a linear operation and that can be inverted. Namely, we can define the following regularized inverses We can then do the following substitutions in Eq. (12) and (14) ∂m to define the quantity back-propagated by target propagation. Note that by taking the gradient of the inverse we can ignore the biases and the inputs. Namely, we have for Summary. Combined together, we get, denoting d t =∂ L ∂ht , ∂L This provides a rule to propagate targets through linearized regularized inverses.

C Target Propagation vs Gradient or Gauss-Newton descent C.1 Gradient Back-Propagation vs Target Propagation
Lemma 3.1. The difference between the oracle returned by gradient back-propagation ∂ θ h (y, g θ (x 1:τ )) and the oracle returned by target propagation can be bounded as For regularized inverses, we have, Proof. The first claim is a direct application of Lemma C.1 and the second claim follows from the formulation of the regularized inverse, using that ∇a −1 (h t ) = ∇a(a −1 (h t )) −1 = ∇a(u t ) −1 .
Lemma C.1. Given A 1 , . . . , A n , B 1 , . . . , B n ∈ R n×n , for any matrix norm · , and any 1 ≤ t ≤ n, A convergence to a stationary point for TP can be derived from classical results on an approximate gradient descent detailed below (the proof is akin to the results of Devolder et al. (2014)).
Corollary C.2 (Corollary of Lemma C.3). Denote ε k a bound on the difference between the oracle returned by gradient back-propagation and by target-propagation both applied to the whole dataset. Provided that the objective is L-smooth and the stepsizes of TP are chosen such that γ = γ h γ y < 1/2L, after k iterations, we get min i∈{0,...,k−1} the output of the RNN on a sample x i and the chosen loss. Lemma C.3. Let f : R d → R be a L-smooth function. Consider an ε-approximate gradient descent on f with stepsize γ ≤ 1/(2L), i.e., x k+1 = x k − γ ∇f (x k ), where ∇f (x k ) − ∇f (x k ) 2 ≤ ε k . After k iterations, this method satisfies, for c 1 , c 2 two universal constants, min i∈{0,...,k−1} Proof. Denote g k = ∇f (x k ) − ∇f (x k ) for all k ≥ 0. By L-smoothness of the objective, the iterations of the approximate gradient descent satisfy where in the last inequality we bounded the absolute value of the last term and used that γL ≤ 1. Now we use that for any a, b ∈ R and θ > 0, 2ab ≤ θa 2 + θ −1 b 2 , which gives for θ > 0, a = γ(1 − Lγ)/2 ∇f (x k ) 2 and b = γ(1 − Lγ)/2 g k 2 , for L a squared loss. Namely, we keep the inversion of the gradients of the intermediate functions.
Our objective here is to question whether viewing TP as a GN step with the approximations explained above is meaningful or not.
Does the original TP formulation approximate GN? Meulemans et al. (2020) start by considering the original TP formulation, i.e., targets computed as v t = ψ t (v t+1 ) for ψ t an approximate inverse of f t and with v L = x L − η∂ x L (y, x L ). Meulemans et al. (2020, Lemma 1) show then that, provided that we use the exact inverse, (Meulemans et al., 2020, Theorem 2) conclude that (i) for mini-batches of size 1, (ii) for a squared loss, (iii) for invertible f t , as η → 0, TP uses a Gauss-Newton optimization with block diagonal approximation to compute the targets in the sense that as η → 0, . As the stepsize of any optimization algorithm tends to 0, they all are the same, since the update would be 0 in all cases. An optimization algorithm aims not to have infinitesimal stepsizes. To make the claim of Meulemans et al.
(2020) more precise, the constants hidden in O(η 2 ) need to be detailed in order to understand in which regimes of the stepsize the approximation is meaningful. Assuming the inverses ψ t to be Lipschitz continuous and L-smooth (i.e. with L-Lipschitz continuous gradients), a quick look at the proof of Lemma 1 of Meulemans et al. (2020) The above bound shows that ξ t 2 grows w.r.t. the stepsize η as a polynomial with leading term η 2 L−t . So unless the stepsize is extremely small, it seems unclear whether the original TP formulation approximates GN in this case. Though the above bound may be pessimistic, it captures correctly the dependency of the error w.r.t. η. True, as η → 0, the leading term is in O(η 2 ), but generally stepsizes are computed for η not infinitesimally small.
Finally, if the similarity of TP with GN could explain its efficiency, then by the reasoning of Meulemans et al. (2020), the original TP formulation should be efficient. Yet, the original TP formulation has never been shown to produce satisfying results.
Does TP with the difference target propagation approximate GN? Meulemans et al. (2020) make a similar claim for TP with the Difference Target Propagation formula, i.e., v t = x t + ψ t (v t+1 ) − ψ t (x t+1 ). Namely, Meulemans et al. (2020, Lemma 3) show that Once again, for the claim to be meaningful beyond infinitesimal stepsizes, the terms in O(η 2 ) need to be detailed. A quick look at the proof of Meulemans et al. (2020, Lemma 3) shows that under appropriate smoothness assumptions the error can be bounded as a polynomial in η with a leading term η 2 L−t . So again, unless we consider infinitesimal stepsizes, it is unclear whether this approximation is useful.
Linearized target propagation and GN. If we use a linearized version of the difference target propagation formula as presented in (2) and the idea that TP could be seen as an approximate GN method may be pursued in a meaningful way. However the error of approximation of the inverse of the gradients must be taken into account in order to understand the validity of the approach.
Propagating the approximation error of the gradient inverses. We compute the approximation error incurred by composing gradients of the inverse instead of inverses of gradients. Formally, the approximation error for one layer can be estimated under the assumption that with e an ε-Lipschitz continuous function and the assumption that the minimal singular value σ of ∂ xt f t (x t−1 ) is positive.
The function e a priori depends on θ; we ignore this dependency and simply consider e to be ε-Lipschitz continuous for all θ. For a function e, we define its Lipschitz continuity constant as ε = sup x sup λ 2≤1 ∂ x e(x) λ 2 = sup x ∂ x e(x) , where · denotes the spectral norm. By differentiating both sides of Eq. (15), we get By assuming the minimal singular value σ of ∂ x f t (x t−1 ) to be positive, we get that ∂ x f t (x t−1 ) is invertible and so
Testing the hypothesis that TP could be interpreted as using GN updates directions. Here we come back to the setting of RNNs presented in the paper. In this case the length of the compositions of layers is τ and according to the previous discussion, the error of approximation of the product of the inverse of the gradients by the product of the gradients of the approximate inverses could easily diverge as τ grows (long sequences). Nevertheless, by using analytical formulas for the inverses, we can ensure that the approximation error is zero, which would correspond then to the ideal setting where TP uses GN update directions for the hidden states.
The updates of TP using the formula (2) would then be exactly the ones of a GN update direction, i.e., So by considering our implementation without regularization, we can test whether the interpretation of TP as an approximate GN method is meaningful in terms of optimization convergence. As shown in Fig. 6b, it appears that regularizing the inverses is necessary to obtain convergence, hence the interpretation of TP as GN may not be sufficient to explain why TP can converge.

D.1 Initialization and hyper-parameters
Initialization and data generation. In all experiments, the weights of the RNN are initialized as random orthogonal matrices, and the biases are initialized as 0 as presented by Le et al. (2015) and Manchev and Spratling (2020). For all experiments, the data was not normalized, as done by Manchev and Spratling (2020). We kept a setting as similar as possible as the one of Manchev and Spratling (2020) to be able to compare target propagation with regularized or parameterized inverses.
Hyper-parameters. In the synthetic tasks, for BP we used a momentum of 0.9 with Nesterov accelerated gradient scheme as done by Manchev and Spratling (2020). Otherwise, we did not use any momentum for the experiment on MNIST pixel by pixel presented in the main paper. The learning rates of BP and the parameters of TP were found by a grid-search on a log 10 basis and are presented in Table 1. We did not add a regularization term in the training of the RNNs. For the Fig. 6b, we used batch sizes of size 512 and performed a grid search for the stepsizes of BP and for the stepsizes γ h of TP while keeping the same regularization r and stepsize γ θ to the parameters found for the length 784.
Software. We used Python 3.8 and PyTorch 1.6. The RNN was coded using the cuDNN implementation available in PyTorch that is highly optimized for computing forward passes on the network or gradient back-propagation.
Hardware. All experiments were performed on GPUs using Nvidia GeForce GTX 1080 Ti (12G memory). Each experiment only used one gpu at a time (clock speed 1.5 Ghz).
Time evaluation. On our GPU, we observed that for the MNIST pixel by pixel experiment, 200 iterations (each iteration considering 16 samples) were taking approximately 60s for BP and 800s for TP. Note that with larger batchsizes the cost of the regularized inversion would be amortized by the fact that more samples are treated simultaneously. We kept the setting of Manchev and Spartling (Manchev and Spratling, 2020) for ease of comparison.

D.2 Additional experiments
The overhead of TP can be worth its performance. To account for the additional cost of inversion for each minibatch, we consider the convergence of the algorithms in time rather than in iterations. We found that, on average, 1 iteration of BP takes approximately 13 times less time than one iteration of TP in our implementation (note that BP benefits from highly optimized implementations for GPU machines, and TP could potentially also benefit from the same optimized implementations). Therefore we ran BP for 13 times more iterations than TP and multiplied the number of iterations by the approximate time needed for each iteration for all algorithms. In the right panel of Fig. 7, we observe that in time too, TP performs better than BP, which stays stuck at an accuracy of approximately 22 Regularized inverses outperform parameterized inverses. We evaluate the impact of using regularized inverses as opposed to parameterized inverse and linearized propagation as opposed to finite-difference-based propagation. The variant of target propagation with parameterized inverse and finite-difference propagation corresponds to the approach of Lee et al. (2015) recently implemented by Manchev and Spratling (2020) and referred to in the figure above as DTP-PI. The variant of target propagation with regularized inverse and finite-difference propagation is referred to in the figure above as DTP-RI. Recall that our approach involves regularized inverses and linearized propagation, referred as TP. In Fig. 7, we observe that both TP and DTP-RI outperform DTP-PI, demonstrating the benefits of using regularized inverses. On the other hand, both TP and DTP-RI perform on par overall, with the former being slightly better for the given parameters.
Target propagation is robust to the choice of the target stepsize γ h . In Fig. 6a, we observed how the convergence could be affected by the choice of the regularization and the step-size γ θ . In the left panel of Fig. 8, we observe that varying γ h does not lead to significant changes in the convergence behavior.
Target propagation does not benefit from momentum techniques. Numerous methods have been proposed to enhance the performance of a classical stochastic gradient descent by using, e.g., a momentum term akin to Nesterov's accelerated gradient formula (Sutskever et al., 2013). Since TP also produces a priori a descent direction, we can wonder whether an additional momentum provides faster convergence. In the right panel of Fig. 8, we observe that adding a momentum on TP is possible but does not seem to provide significantly faster convergence while being less stable. On the other hand, on this same figure, the momentum seems to help the gradient descent. Our preliminary experiments using Adam with TP did not conclude; namely, we were not able to obtain a convergence similar to the one illustrated in Fig. (4) that simply used Algo. 2. We leave for future work the implementation of appropriate adaptive stepsizes strategies for TP.