{'fraction_non_alphanumeric': 0.031639013917399574, 'fraction_numerical': 0.007664927776417635, 'mean_word_length': 4.264723378941106, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Deriving Differential Target Propagation from Iterating Approximate Inverses

We show that a particular form of target propagation, i.e., relying on learned inverses of each layer, which is differential, i.e., where the target is a small perturbation of the forward propagation, gives rise to an update rule which corresponds to an approximate Gauss-Newton gradient-based optimization, without requiring the manipulation or inversion of large matrices. What is interesting is that this is more biologically plausible than back-propagation yet may turn out to implicitly provide a stronger optimization procedure. Extending difference target propagation, we consider several iterative calculations based on local auto-encoders at each layer in order to achieve more precise inversions for more accurate target propagation and we show that these iterative procedures converge exponentially fast if the auto-encoding function minus the identity function has a Lipschitz constant smaller than one, i.e., the auto-encoder is coarsely succeeding at performing an inversion. We also propose a way to normalize the changes at each layer to take into account the relative influence of each layer on the output, so that larger weight changes are done on more influential layers, like would happen in ordinary back-propagation with gradient descent.


Introduction and Earlier Work
In recent years there has been a revival of interest in exploring analogues of backpropagation which could be more biologically plausible than back-propagation, for which concerns were voiced early on [Crick, 1989]. This paper introduces ideas and mathematical derivations for a family of such biologically motivated analogues, in line with earlier work on "target propagation" [LeCun, 1986, Lecun, 1987, LeCun et al., 1988, Krogh et al., 1990, Bengio, 2014, Lee et al., 2015b, Lillicrap et al., 2020, in which a feedback path provides targets for each layer, and each layer-to-layer forward and feedback paths form an auto-encoder, respectively with its feedforward path encoder and corresponding feedback path decoder. This paper attempts to provide a solid theoretical foundation for target propagation by studying the case where one can use the auto-encoder to approximately invert the targets of the next layer (l) in order to obtain a target for the previous layer (l − 1).
See Lillicrap et al. [2020] for a review of alternatives or analogues of backpropagation aimed at bridging the gap between backpropgation and biology, and in particular of how target propagation and difference target propagation [Lee et al., 2015b] could potentially be implemented in the brain. The main concerns about a biological implementation of back-propagation can be summarized as follows: 1. The weight transport or weight symmetry problem is that back-propagation requires the backprop calculation (going backwards from outputs towards inputs) to use the exact same weight matrices (transposed) as those used for the forward propagation, but these identical synaptic weights would have to be on different neurons.
2. The non-linearities used in the forward and backward phase of backpropagation are not of the same type and in addition those on the backward path need to be numerically coordinated with the corresponding ones in the forward path (since they should multiply their input by the first derivative of the corresponding forward non-linearities).
3. The brain has recurrent connections (both feedback and lateral ones) which do not disappear while the signals travel on the feedforward paths.

Backpropagation needs a clean separation between examples and between
the forward phase and backward phase, with the forward phase activations somehow being held while the backward phase proceeds because a combination forward phase and backward phase quantities is needed to perform the synaptic weight updates. Various analogues of backpropagation such as difference target propagation and the proposal made here require a form of time-multiplexing and synchronization in order to separate the feedforward and feedback phases, while the brain seems to use feedback paths all the time.
5. Dale's law constrains the signs of the synaptic weights (depending on the neuron type), while they are unconstrained in most artificial neural network implementations, including those aiming at biological plausibility.
6. Neurons communicate with spikes, not with real values.
7. Backpropagating through time (BPTT, used to handle recurrent networks) requires (a) storing all past states of the brain and (b) traversing these states backwards to compute updates (c) waiting for a long time before making these updates (to be able to capture sufficiently long dependencies). Instead, the brain seems to make online updates (or possibly daily updates for longer-term consolidation during sleep).
This paper (and target propagation in general) helps with the problem (1) of weight symmetry or weight transport [Grossberg, 1987, Crick, 1989, Liao et al., 2015, the problem (2) of having to match non-linearities as well, the problem (3) of explaining (to some extent) the feedback and lateral recurrence, and it suggests a way to deal with the problem (4) of synchronization and phases (using neuromodulating gating, e.g., with acetylcholine). This paper does not discuss spiking (but rates can approximately be obtained in various ways by averaging over time and groups of neurons), Dale's law nor BPTT. Many other approaches have been proposed to deal with the weight transport problem, notably feedback alignment [Lillicrap et al., 2016], by showing exact transposes are not always needed (although the approach is still challenged by more complex tasks like ImageNet object recognition [Bartunov et al., 2018]). This is in addition to target propagation approaches [Lecun, 1987, Bengio, 2014, Lee et al., 2015b, Ororbia and Mali, 2019, Lillicrap et al., 2020, which train auto-encoders to obtain feedback weights. Like many other approaches with a feedback path, the approach discussed here assumes that the feedback path mainly influences the computation of targets for the feedforward path. However, variations proposed here use the recurrent computation to settle to approximate inverses of the output target at each layer (so that moving the feedforward activations to these targets would produce the desired output on the final layer). Like Equilibrium Propagation [Scellier and Bengio, 2017], this paper addresses the potential mismatch of non-linearities between the forward and backward phase, here by appropriately training or iterating over the feedback paths (whereas Equilibrium Propagation takes advantage of the propagation of perturbations in the dynamics of a recurrent network with a static input). Like Equilibrium Propagation, the proposed iterative procedure for target propagation is a relaxation, but one that provably converges exponentially fast, with a rate depending on the reconstruction ability of the auto-encoder (hence the need to train the feedback weights), and it avoids the weight transport problem associated with Equilibrium Propagation. Like Podlaski and Machens [2020], this paper proposes to use iterative refinements involving the encoder and decoder at each layer to obtain and then propagate approximate inverses, although Podlaski and Machens [2020] use fixed random feedback weights whereas we use a mechanism that requires trained feedback weights (minimizing a reconstruction error) in order to achieve fast convergence to a good inverse. Finally, this paper borrows the idea from GAIT-prop [Ahmad et al., 2020] that the targets propagated in target propagation should be only small perturbations of the corresponding feedforward activations to facilitate an analysis of the algorithm based on a linearization of the forward propagation around the forward phase computation. In comparison to GAIT-prop, this paper avoids a potentially exponential growth of required learning rates as one considers layers further from the output, and the proposed update rule is justified and derived in a different way, although also starting from the premise of approximate inverses being computable (in our case via auto-encoder iterations extending the original idea in difference target propagation).

Layer-wise approximate inverse
We focus for now on a multi-layer feedforward model and assume that all the layers have the same size, to make the layer-wise inversion feasible and simple. See Podlaski and Machens [2020] for ways to generalize target-propagation methods to the case of unequal layer sizes. See also Section 9.3 to deal with skip connections and architectures which are not structured as a chain. Each layer's forward phase activation h l is computed with its encoder f l : where f l is implicitly parametrized by W l , the parameters (weights and biases) used for computing the output of layer l, h 0 = x is the input vector and h L is the output layer, with loss L(h L , y). The second notation is convenient to describe compositions of these layer-wise functions, so In addition to learning the feedforward weights W , we consider neural networks with a separate feedback path with top-down layers g l and corresponding feedback parameters Ω l to compute an approximate inverse of f l : Note that inverses compose, so and simple target propagation starting at τ L (normally chosen close to h L ) gives rise to or equivalently, applying a change τ l − h l on layer activations h l (to give new values τ l ) would lead to which would change h L in the direction of the gradient of the loss on the output. Note that the inverse condition is achieved when we minimize reconstruction error in the output space (applying the decoder and then the encoder), which is the opposite of the usual auto-encoder reconstruction error (where we apply the encoder and then the decoder). However, the iterative inversion scheme of Section 6 uses the regular auto-encoder property of g l • f l being close to the identity, and the version in Section 10 minimizes the reconstruction error of g l • f l but obtains τ l−1 → f −1 l (τ l ) iteratively. We will assume that g has enough capacity for approximating the inverse of f . For example it may need an extra hidden layer or use powerful nonlinearities assumed to exist in dendritic computation in order to obtain the required flexibility, or the use of iterative refinements to correct the input or output of the decoder in order to approximate the inverse well, as per Sections 5 and 6. A training objective for enforcing the above also needs to be defined (such as squared error on reconstructing pre-activation values).

Differential Target-Prop when a Single Layer is Changing
We will first study the case where a single layer l is to have its weights updated in order to make the output h L match the target output τ L . We will also define differential targets τ l for each layer, i.e., small perturbations of the forward activations of each layer obtained by applying and composing the approximate inverses g l on the differential output target τ L : where 0 < β ≪ 1 is a small positive constant: small enough to linearize g and f and large enough to make computations with limited precision feasible (which possibly may require several neurons to represent a single quantity with enough precision, in the brain). Because of the availability of an approximate inverse, we can propagate targets and obtain τ l for layer l using Eq. 7. It has the property that if h l moves to τ l , we obtain τ L as output: for β small enough (it acts like a learning rate). We can move parameters W l so that h l approaches τ l by considering τ l as a local target for h l , e.g., using what we now call the single layer update Differential Target-Prop (DTP1) update, either by moving activations or with parameter updates (where the update differs from difference target propagation and other proposals which do not have the normalization factor below): The inverse in the above expression is a normalization factor which goes beyond the usual gradient of the squared error (τ l − h l ) 2 with respect to W l , computed by the rest of the expression. The role of this inverse matrix is to make sure that when we apply ∆ DTP1 W l to W l , we obtain the change (τ l − h l ). Note that the matrix is normally block-diagonal with a block per neuron, because the output of one neuron only depends of its incoming synaptic weights and not (directly) of the synaptic weights of the other neurons. Now let us see what the matrix inverse factor provides by computing the change ∆h l (∆ DTP1 W l ) which would result if ∆ DTP1 W l was added to W l . Let us write h l as a function of W l (keeping h l−1 fixed): Substituting Eq. 10 in the above, we obtain that the change in h l due to the DTP1 change in W l is which was the objective. In other words, by using the normalization factor, we can recover the desired target τ l on layer l to make h L go to τ L . Whereas we use an explicit feedback path to propagate the targets down across layers (such that a change in h l leads the desired change in h L ), we use the weight update formula (Eq. 10) to propagate the target change τ l − h l into a corresponding target change in W l , such that applying that parameter change would yield the desired result on h l , and thus on h L . For example, with f l (h l−1 ) = W σ(h l−1 ) for some activation non-linearity σ(), we would get the following update (shown below, see Eq. 14): which is similar to the delta rule except for the normalization of the activation vector of the previous layer. Note also that (τ l − h l ) is already proportional to β and that β acts like a global learning rate. One could however also have local adaptation equivalent to a local learning rate, which could scale or modify the above equation (e.g. using a momentum method). In Section 7 we derive a perlayer learning rate to make the updates proportional to the influence of W l on L.
Let us now derive Eq. 13 from Eq. 10 for the special case f l (h l−1 ) = W σ(h l−1 ), using the Einstein notation to simplify the tensor calculations: where δ i,j = 1 i=j , and for the weights W l,i,j into unit i, as claimed in Eq. 13. The above suggests that if we can move parameters so that h l approaches τ l , we would expect the resulting network to lead a smaller loss on (x, y). We prove a cleaner and stronger result below about how this rule, if the inverses are well approximated, leads to a change τ L − h L = −β ∂L(hL,y) ∂hL on the output layer and how this relates to a regular gradient step through the implicit inverse of a Gauss-Newton matrix.

DTP1 Updates Approximate Gauss-Newton
Let us try to analyze the ∆ DTP1 learning rule in the limit of β → 0 which linearizes the encoder and decoder computations.
We will start from the per-layer approximate inverse of Eq. 3 and its consequence Eq. 6 on the whole path from an intermediate layer h l to the output layer h L . Let us rewrite Eq. 6 as follows: where ∆ DTP1 h l = τ l − h l . If we do a first-order Taylor expansion of the above equation around h l (considering ∆ DT P 1 h l is small when β → 0) we get Let us write the shorthand for the Jacobian matrix for the mapping from h l to h L , so Now let us consider the true gradient of the loss with respect to h l , scaled by −β (which would be the SGD step if we were able to perform one in the space of hidden activities h l ): where for the second line we use the chain rule, for the third line the definition of τ L in Eq. 7, and for the fourth line we replace τ L − h L as per Eq. 19, so we have Now let us compare the above with the change in h l which the DTP1 update (Eq. 10) is trying to achieve: where we make use of Eq. 20. Hence what we see is that DTP1 is pushing the weights so that they move h l not in the direction of the gradient − ∂L(hL,y) ∂h l but instead in the direction corresponding to the Gauss-Newton direction for updating h l . What this provides is a form of invariance with respect to the joint effect of the parametrization of the forward path and the representation h l : we want that the change in the resulting function (i.e., in h L ) be the same if we transform both h l and the forward path while keeping the same h L . This is a general property of target propagation since it focuses on finding the change (in h l or in W l ) which achieves a target τ L , irrespective of the specific parametrization in between. Indeed, if we were to instead perform a regular gradient descent step with respect to h l , h l would move in the direction ∆ SGD h l and the outcome would be a change in h L which is not aligned with τ using the third line of Eq. 20. Note that the main property enjoyed by Gauss-Newton updates, i.e., independence of the parametrization, is also obtained with the parameter update ∆ DTP1 W l and not just by the pseudo-update of the activations ∆ DTP1 h l , because by construction the former cause the latter, as shown in Eq. 12.

Iterative Refinement of an Approximate Inverse by Updating a Correction on the Input of the Decoder
The main concern one may have about the above derivations is that any reasonable approximate inverse, especially if it is learned (e.g., with a form of per-layer auto-encoder objective) is still going to make errors, and that these errors can be composed and grow across layers. What if the error made by the approximate inverse g l is comparable to the targeted change in h l ? We consider two approaches to reduce these potential challenges, which can be combined. First, in this section, we see how any approximate inverse can be iteratively improved at an exponential rate by introducing a gradually improved correction at the input of the decoder. Second, in the next section, we consider how to iteratively obtain a correction at the output of the decoder, again showing exponential convergence at a rate which depends on how good the decoder is at inverting the encoder.
To simplify the notation, let us consider in this and the next section a simple setting with a forward function f with input x and approximate inverse g, with f (g(y)) ≈ y. Given a target value τ y for the output of f , we are looking for an input target τ x such that Let us introduce a modified argument to the decoder, u, which may initially be τ y but which we will iteratively modify such that τ x = g(u) leads to i.e., We update u as follows: where r is the map which transforms u iteratively. If this iteration converges, it means that we have τ y = f (g(u)) as per Eq. 26 in order to obtain the desired inverse, i.e., Eq. 25.
Let us now analyze this relaxation under the assumption that the autoencoder only deviates slightly from a perfect inverse, meaning that f (g(u)) − u has a Lipschitz constant α which is smaller 1 (ideally much smaller than 1), i.e., that its Jacobian is close to the identity, with E a small matrix, i.e., with its largest eigenvalue ρ(E) smaller than a constant α which is also the Lipschitz constant of f (g(u)) − u: Using the Lipschitz continuity of r with constant α, we obtain using Eq. 27 (with iterates indexed by t) from which we can conclude that u t converges exponentially fast to 0 in α t , where α is closer to 0 when we have a better inverse approximation and lower reconstruction errors. Note that convergence implies a perfect inverse is obtained. This analysis also tells us the condition for convergence, i.e., that the spectral radius of the auto-encoder Jacobian deviates from the identity at most by an α < 1. In practice, this condition is easy to obtain with auto-encoders trained to minimize a form of reconstruction error. To see this, note that α = 1 would mean that the reconstruction error is of the same order of magnitude as the reconstruction itself, whereas α = 0.1 would roughly correspond to a 10% error. Although we developed the iterative updates with an intuition based on the feedback weights already doing a good job, the above analysis shows that it only need to do a reasonable job. As soon as reconstruction error gets to be somewhat minimized, we should be able to achieve the convergence condition easily. From a biological perspective, initializing the connectivity pattern to obtain approximate matching of feedforward paths by corresponding feedback paths might already achieve a sufficient condition for the iterations to converge even before training kicks up. Furthermore, minimizing reconstruction error is something which can be done offline, even in the absence of real data (e.g. in the womb).
6 Iterative Refinement of an Approximate Inverse by Updating a Correction on the Output of the Decoder Instead of correcting g by changing its input, we can correct it by changing its output, as discussed here.

Single-Step Refinement
Once we have obtained a pair (x ′ , y ′ ) with y ′ = f l (x ′ ) such that y ′ ≈ τ l (e.g., using the procedure in the previous section) we can make an extra correction which focuses on the error in the space of h l−1 , and in the next section we generalize this idea to obtain an iterative refinement of the inverse, making a correction on the output side of the decoder.
To lighten notation, let us again focus on a single layer with an input x, output y, encoder f , decoder g, and output target τ y near y. The optimal additively adjusted input target τ x is such that f (τ x ) = τ y . If g misses the optimal τ x by some error and g is smooth, then we can assume approximately the same error at a nearby point From this estimated local auto-encoder error vector, we can approximately solve for f −1 (τ y ): Note that this works to the extent that f and g are smooth enough (and for both of them to be smooth, we want their Jacobians to have singular values close to 1, which is good for back-propagation as well). If we had obtained τ x by iterating on the decoder input correction δ t as in the last section, then, instead of estimating the error g(τ y ) − f −1 (τ y ) at τ y using as a proxy the error ) − x, we can use the pair (x ′ , y ′ ) with y ′ closest to τ y we found (normally at the last iteration). Using the notation from the previous sections, this would give us the following corrected update: where t is the last step of the iterative refinement and δ l,t the last increment obtained when trying to invert f l at τ l .

Iterative Correction of the Decoder Output
We can iterate the type of correction proposed in the previous subsection in order to exponentially converge to an inverse, assuming that the regular auto-encoder minus the identity is sufficiently contracting: where we defined the iterate map m and we can initialize τ l−1 ← g l (τ l ). At any step, τ l−1 represents a value of the activations at layer l − 1 which was found to be the best to make f l (τ l−1 ) close to τ l . We could also combine the above formula with the one in the previous section and initialize τ l ← g(τ l + δ l,t ) for the last step t in the iterations of the previous section. We can analyze the convergence of Eq. 35 by using ideas similar to those in the previous section. The trick is to look for a quantity whose iterates have a rapidly converging (i.e. strongly contractive) map. We again assume that the auto-encoder does a good job, but this time we need to consider the g • f regular auto-encoder rather than the f • g reverse auto-encoder. We will assume that the Jacobian of g • f is near the identity, or more precisely that m(u) has a Lipschitz constant α < 1 (which is also a bound on the spectral radius of m ′ (u)), noting that m ′ includes the identity minus the Jacobian of g • f , which only leaves the Jacobian of the reconstruction error: and we correspondingly assume that the Lipschitz constant of m is no more than α. With this Lipschitz condition and α < 1, we obtain (after noting the iterates of τ l−1 at step t as u t ) that Hence the series of increments (u t+1 − u t ) of τ l−1 converges and it converges exponentially fast, in α t . Note that at convergence, g(f (u t )) → g(τ y ), which means (if g has an inverse) that f (u t ) → τ y . This analysis suggests that if we use this iteration to find a good target for the layer below, we should make sure to minimize reconstruction error of the regular auto-encoder g • f .

DTP for Multiple Layers Updated in Parallel
Because β is small, we can consider small variations of h l to operate on a linearized forward propagation computation, and the combined effects of all the increments ∆ DTP1 h l = τ l − h l applied to each layer l on the network output h L simply add up to yield a change on the output layer. As we have seen in Section 3 if h l moves to τ l this will move h L to τ L . We also saw that the DTP1 weight update would move W l such that h l moves to τ l , and thus that h L moves to τ L . Hence doing a DTP1 weight update in parallel over all the layers yields a change L(τ L − h L ) on the output. We could simply scale down each layer update by 1/L but this would give equal influence on h L to all the layers. What should be the relative impact of each layer on the output (and the loss)? This is the question we study in this section.
If we were to do ordinary gradient descent, the layers with the largest gradient would move proportionally to the norm of their gradient: layers with a small gradient (i.e., less influence on the loss) would move less. Instead, if we apply DTP1 (even scaled by a constant 1/L) to all layers in parallel, the layers with the most influence on the loss (which have a smaller target change τ l − h l sufficient to induce the same change τ L − h L in the output) would be changed the least! We do want to keep the invariance property associated with Gauss-Newton within each layer, but if we do it separately for each layer we might be doing something wrong. Also note that any convex combination (and not just 1/L) of the DTP1 updates would produce the desired target on h L . What is the right weight to put on each of several changes each of which can individually achieve the goal? We propose putting more weight on those layers which have more influence on the loss, just like in gradient descent.
To elucidate this in a simple case, consider two adjustable input variables x 1 and x 2 (for now think of them as scalars but later we will think of them as vectors) and a scalar (loss) function L(x 1 , x 2 ). What is the relative magnitude of the loss gradients with respect to each of them? The gradient with respect to x i is the change in the loss ∆L(∆x i ) induced by a change ∆x i , divided by ∆x i , in the limit when the magnitude of ∆x i goes to 0, and the magnitude of this ratio seems to be a good measure of the influence of x i on L. To be able to compare both limits (for i = 1 and i = 2) at once in terms of their effect on the loss, let us consider a scaled input change κ∆xi ||∆xi|| where we obtain the comparable limits by letting κ go to 0 simultaneously for both changes in x 1 and x 2 . Now we obtain same size changes in both x 1 and x 2 , so that the magnitude of their relative gradient (i.e., their relative influence on the loss) can be read off in the corresponding magnitude of the change in loss. As we show below more formally, the ratio of the influences of the two components should be proportional to the ratio of the induced changes in loss for a same-size change in x i .
Given some directions of change ∆x i , one should ask how should we scale them into ∆x i ∝ ∆x i so that we can update all the x i in those directions but with the magnitude of the change, || ∆x i || being proportional to their "influence" on the loss, i.e., a term we formalize here as a measure of the influence of x i on L, which is | ∂L ∂xi | in the scalar case. We thus get for very small changes (scaled by a small β). The influence(∆x i , L) function here measures the magnitude of the effect of changing component x i on the loss L. It is the absolute value of the gradient ∂L ∂xi in the case of scalar x i but otherwise is defined by the above equations. We need to introduce a β to scale ∆xi ||∆xi|| or otherwise the change in x i is too large for the resulting change in L to indicate only first-order effects. To compensate we thus need to also divide by β outside.
Hence the relative influence of the different components x 1 and x 2 on the loss is Now, in our case, ∆x i can either be interpreted as a layer change (under DTP1) or a weight change (also under DTP1). In both cases, the corresponding output change is the same for all l (Eq. 8 and Eq. 12), moving h L to τ L . To change h l or W l proportionally to their influence on L, let us consider what is that influence: where we see that the influence is indeed inversely proportional to ||∆ DT P 1 h l || but that it also includes a factor || ∂L ∂hL || 2 which is the same for all layers for a given example (x, y) but which varies across examples. In the same way that we prefer to make larger changes on layers with a larger gradient, we prefer to make larger changes on examples with a larger gradient. Hence, this analysis is also telling us how to scale the parameter updates to normalize properly across layers and across examples.
To obtain this, we propose that changes in h l or W l be scaled to have as magnitude their influence, times the global learning rate β: where we have introduced ∆ DT P , the differential target-prop update taking into account the issue of normalizing across layers. Note how the only thing which is needed to be shared across layers it the norm of the global output loss gradient ||τ L − h L || 2 = β 2 || ∂L ∂hL || 2 . Similarly, we obtain the target parameter change for layer l: Now, when we update all the layers in parallel with DTP, they each contribute in proportion to their influence on the loss. The special case f l (h l−1 ) = W l σ(h l−1 ) gives the update At this point, a natural question is the following: what target should layer l propagate upstream? The target corresponding to ∆ DTP1 h l which guarantees that lower layers are also trying to make h L reach τ L , or the scaled down target change ∆ DTP h l which is weighted to make it proportional to the influence of h l on the loss? Doing the calculation shows that if we propagate ∆ DTP h l as a target for lower layers and not account for the scaling factor ||τL−hL|| 2 ||τ l −h l || 2 , this will result in updates that have this extra factor in them. It would thus require passing these scaling factors down along with the feedback path in order to obtain the correct updates in lower layers. See Section 9.4 for a related discussion.

Iterating over Multiple Layers in Parallel
Up to now we have only considered iterating the inverse calculation on a perlayer basis. This would require waiting for each inverse to be performed on layer l to obtain a target for layer l − 1 and start iterations for layer l − 1. However, this is probably neither efficient nor biologically plausible. What we propose instead is the simultaneous optimization of targets at all layers in parallel. Of course, at least L steps need to be performed for the information to propagate from τ L down to τ 1 . Nonetheless, this procedure should converge faster, since middle layers can start searching for better targets before the upper layers have converged. This makes sense in the context that the proposed iterations are rapidly converging and gradually provide a more and more appropriate target change at each layer. Also note that from a biological perspective, it is plausible that local iterations can be performed quickly because they involve short local paths. This general idea is illustrated in Algorithm 1.
Algorithm 1: Example of parallel multi-layer iterative Differential Target Propagation and parameter update scheme, for one example (x, y). The network has L layers, each with encoder f l , producing feedforward activation vector h l = W l σ(h l−1 ) with activation non-linearity σ. Biases can be implemented by augmenting the vector σ(x) with an additional element with value 1. Each layer also has a decoder g l (u) = Ω l σ(u), used in propagating targets and iteratively obtain τ l−1 for the layer below, from target τ l for layer l. It is advisable to use an invertible non-linearity such as the leaky ReLU, to facilitate the task of the decoder. The targets are iteratively updated in parallel until the lower layer converges to a given StoppingPrecision hyper-parameter. The overall (input,output) example is (x, y) (note that the non-linearity σ is going to be applied to x so some preprocessing may be appropriate) and there is a loss function L(h L , y) whose gradient with respect to h L is used to initialize the output target τ L , with a global learning rate β ≪ 1. This algorithm shows the weights updated at the end of the relaxation phase but they can also be updated on-the-fly with differential equations as per Section 10.1. The ∀l statements can be executed in parallel.

Training the Decoders and Architecture Issues
Up to now we have not discussed an important piece of the puzzle, i.e., the training of the decoder, as well as its architecture.

Expressive Power of the Decoder
Note that the proposed algorithm requires a separate path with separate neurons for computing the inverses, with these inverses then providing targets for the corresponding feedforward neurons. We have noted that the inverse of a single non-linear layer may require more expressive power than just the same parametrization upside-down. It is thus important that the chosen architecture for the decoder g l at each level has enough expressive power to approximately invert f l , for any value of the forward parameters W l . For example, the decoder neurons may have non-linearities in its dendrites which enable to compute the equivalent of an MLP inside a single neuron (i.e., no need to have an external backprop, but one should then imagine show such gradients could be computed inside the non-linear dendrites). On the other hand, if we use one of the iterative procedures for performing the inverse, the decoder itself may not need to have the capacity to invert the encoder precisely, so it could also be that a simple mirror of the feedforward architecture is sufficient, in that case: there is a trade-off between the number of iterations needed for useful inversion and the expressive power of the decoder. An interesting extension of the proposal in this paper is to make the decoder have its own hidden layer. This could easily be obtained as follows, in the context of encoder-input correction. Imagine we had another group of neurons which predict u l or the correction u l − τ l as a function of τ l . One way to think of this extra hidden layer is as a way to initialize u l (given τ l ). Another way to think about it is that we now have a 2-layer decoder with a training signal for both layers of the decoder. The advantage of this second decoder layer could be even faster convergence of the inverse, so it could be an interesting trade-off between computation time and computational resources. Note also how this principle of iterating over intermediate layers to achieve a global consistency (here the inversion coherence) might provide a more general way to obtain hidden layer targets, and similar ideas have been explored in the deep learning literature [Carreira-Perpinan and Wang, 2014, Lee et al., 2015a, Taylor et al., 2016, Choromanska et al., 2019.

Reconstruction Objective
An important question is whether the decoder weights should be updated to optimize a regular auto-encoder reconstruction error (g • f should match the identity) or the reverse auto-encoder reconstruction error (f •g should match the identity), or both. Depending on the iterative inverse method chosen, it seems that one or the other of these requirements is more appropriate. The regular auto-encoder error is easier to optimize with respect to the decoder g l since the target for the decoder is readily available in the form of the corresponding feedforward activation h l−1 . Each time f l is computed on h l (which e.g. happens at least in the feedforward phase), there is an opportunity to compute and optimize this reconstruction error. This is the simplest and most biologically plausible solution so we propose it as the first thing to explore.
Using the reverse auto-encoder objective is more tricky because it requires doing the equivalent of back-propagating through the encoder to get a target on the decoder. Any of the target-propagation methods discussed here could be used, noting that in that configuration, g l also acts like an approximate inverse for the encoder f l and so we can use the same procedure for target propagation. What is more tricky is to imagine the timing of the different phases of these operations having to overlap with the computations being performed for the overall target propagation.

Skip Connections and Multiple Paths to Output
What if an inverse is not possible, e.g., because the next layer has more units than the previous one? Similarly, what if the architecture has skip connections (as in ResNets for example)? In those cases, there is no exact inverse, so the best we can hope is to find a target τ l−1 (still close to h l−1 ) which minimizes the reconstruction error of the next-level target, when projected down through g. This is already what Eq. 35 does.
Let us now consider skip connections and more generally the architectures which are not in the form of a linear chain. For example, layer A may feed both layers B and C (and they in turn influence the output layer). The targets τ B and τ C each individually guarantee that τ L can be obtained by h L , if either h A or h B were changed to their respective target. We have an encoder f B which maps the inputs of B (including A) to h B , and we have a decoder g BA which maps a target τ B into a target τ BA which would be sufficient to make h B changed to τ B if h A was set to τ BA = g BA (τ B ) or to the result of an iterative process involving f B and g BA . We have a similar decoder g CA which produces (iteratively or not) a target τ CA for h A sufficient to bring h C to τ C . Since either picking τ A = τ AB or τ A = τ CA is sufficient to make h L reach τ L , and since β is small enough to make the forward maps linearized, it means that setting τ A = κτ BA +(1−κ)τ CA also works for any κ (but we will restrict ourselves to convex combinations, i.e., 0 ≤ κ ≤ 1). What would be a good choice of κ? This situation should remind us of the question of how much weight to give on the different layers (in the linear chain scenario) when changing all the layers in parallel (Section 7). This leads to the following proposal: set κ in the above situation to so as to give a weight to each branch which is proportional to its "influence" on the loss. This principle can be generalized to more than two descendants by using as convex weights ∝ 1/||τ XA − h A || 2 for the feedback from descendant X, with the proportionality constant chosen so the sum of the weights is 1.

Numerical Stability
One concerning aspect of the inverse (iterative or not) and of the Gauss-Newton approximation is that they may lead to numerical instabilities, as we are dividing by quantities (like the activation layer squared norm) which may be small, pushing τ l far from h l or making weights change too drastically. One possibility to avoid such problems is to bound the maximum deviation of τ l from h l in the computation of τ l , and to bound the maximum change in the synaptic weights, but this would possibly introduce biases. An alternative we favor is to scale the target changes at each layer to maintain a uniform magnitude of these changes across layers but keep track of the scaling required to make the weight updates conform to the DTP update. This can be done by propagating both the target (via the feedback weight) as well as a scalar quantity, such that their product corresponds to the DTP target changes. That scalar would then be used as an local learning rate for each layer in order to recuperate the correct DTP weight changes.

Biological Implementation
We have not discussed much how the ideas in this paper could be implemented by biological neurons. See in particular the review from Lillicrap et al. [2020] for a general discussion, especially around target propagation's biological implementation. In terms of architecture, we can assume a similar structure as in that review paper, which also raises the question of how to schedule the different phases, while there is no biological evidence that feedback paths are inactive during recognition, for example. In fact, it looks like feedback paths are used not just for learning but also for forms of attention and gain control. An interesting possibility explored here is that instead of having a sharp distinction between a feedforward phase and a feedback phase (with different paths becoming activated), both are always active but with possibly varying degrees of influence on the feedforward computation.

Updating Weights all the Time
In terms of biological plausibility, it would make sense that the proposed synaptic learning rule be applied all the time and not just at specific moments as is typically the case in a software implementation of learning algorithms. On the other hand, in a software implementation where parallel computation is not completely free, it does make sense to wait for the right moment to perform the weight updates. Another desired feature of a good biologically plausible online updating scheme is that at any point during the settling of propagated targets, it should give good results, as good as the number of iterations already done allows, making the algorithm an anytime algorithm, i.e., producing valid answers whenever it is stopped.
To do this online updating, we need to satisfy two conditions. First, we want the rule to lead to no update in the feedforward propagation phase. This is easy to do by letting the output target be τ L ← h L (so the layer-wise targets are automatically converging to τ l ← h l ) during this phase. Second, while the iterative target computation and target propagation proceed, there will be intermediate values of τ l which are suboptimal.
Since the normalized input n l−1 = σ(h l−1 ) ||σ(h l−1 )|| 2 to layer l does not change while the τ l relax, it does not need to be recomputed in a software implementation and the network can be recomputing always the same thing in a biologically plausible implementation where x is fixed during the whole interval dedicated to processing that example. Alternatively, if network input x keeps changing with time, it would be sufficient to keep a delayed trace of n l−1 in the presynaptic areas, with a delay equivalent to the time for the combined forward propagation / target relaxation and propagation. Note that normalization can be achieved with lateral connections and inhibitory neurons which sum up the squared activity coming from each layer. The weight updates can then be performed as follows: whereτ l is the error signal computed by the feedback machinery, with τ l initialized with h l and constant during the forward phase (during which there is therefore no effective update of the weights) and updated with differential equations according to the chosen relaxing scheme. As a result, the total update to W l is the integral of the above, which ends up being the last value of τ l minus the initial value of τ l , times n T l−1 . See also Proposition 5 for the calculation.

CDTP: a Continuously-Running DTP
To make DTP more biologically plausible, we need to address a number of issues.
One of them is that we would like the same "algorithm" to be run at each unit of the system, all the time. If some phases or global coordination are needed, we need to be explicit about them and identify the signals which need to be communicated across the network to achieve them.
In particular, with DTP as formulated, the same "hardware", the encoder f , has to be used twice but with different inputs: the feedforward input on one hand, in the feedforward phase, and the iterated targets from the lower layer on the other hand, in the feedback phase. There is evidence for gating modulation in the brain between feedforward and feedback signals which could be consistent with such a mechanism [Hasselmo, 2006], and we would like to be more explicit about it here. The results from Hasselmo [2006] suggest that acetylcholine would be involved in the bottom-up vs top-down modulation as well as correlated with synaptic weight changes and its cyclic variations linked to the theta rythm oscillations, which in our interpretation may be linked to the control of feedforward vs feedback phases. We will denote γ for that quantity, with the convention that γ = 1 corresponds to the feedforward phase and γ = 0 corresponds to the feedback phase, with γ continuously oscillating between 0 and 1 as we go from one "example" to the next. γ = 1 must either be sustained long enough for the whole feedforward phase to complete, or more likely from a biological perspective, γ may be different for each layer and a wave of γ = 1 would travel in the feedforward direction while the feedforward activation values (a l below) are computed and stored.

Continuous Update Equations
We propose the following continuous update equations (for now assuming a global γ although a layer-wise travelling wave following the initial feedforward propagation would also work). We write down the equations in discrete time but they could as well be written in continuous time. We will analyze these equations next.
∆W l = 1 α<1 ∆τ l−1 n T l−1 ∆ 1 Ω l = (τ l−1 − g l (u l ))σ l (u l ) T (51) or ∆ 2 Ω l = (a l−1 − g l (u l ))σ l (u l ) T where we have considered two variants of the feedback weights update, ∆τ l is the temporal variation of τ l and n l is the normalized value of σ(a l ) (using lateral connections to compute the normalization of the layer l), i.e., and we have chosen the standard linear+nonlinear parametrization with σ() an arbitrary non-linearity (and a bias can be included by having σ() be vector-valued and produce an extra output equal to 1), as in earlier sections of this paper, although it would be straightforward to define other parametrizations for the feedforward and/or feedback paths.

Analysis of the Continuous DTP Equations
We can verify by inspection of Eqs. 46 and 47 that a l acts as a temporary buffer for the activations of the last feedforward phase (obtained while γ = 1). Conversely, when γ = 0, we obtain a purely feedback phase where the encoder f l ignores a l−1 and instead is used in an auto-encoder loop which composes f l and g l and propagates targets.
using the Lipschitz property in the last line. Finally as γ → 0, from Eq. 49 we get τ l−1 → g l (u l ).
Proposition 5. Ignoring the time during which the feedforward activities are computed, the feedforward weight update (Eq. 50) is equivalent to the following total update throughout the combined feedforward-feedback dymamics: which is the same as ∆ DT P W l from Eq. 43. The initial updates (while γ = 1) are 0.
Proof. We first note that τ l−1 → a l−1 at convergence of u l . The sum of temporal changes in Eq. 50 telescopes because The initial updates (while γ = 1) are 0 because τ l does not move (Eq. 49), being fixed to the value τ l−1 (0) = a l−1 if the feedforward computation has been performed.
While the feedforward activations are settling to their value a l from their old value (for the previous example), the weight change (if performed) would aim at predicting the new value of a l given its previous value, which seems to be a useful feature in a temporal world. However this could also be hindrance (where there is no temporal order) so in Eq. 50 we inserted a 1 α<1 factor to avoid updates until the feedforward phase as settled (note they would be zero once a l 's have settled, so long as α = 1). This analysis is also compatible with the idea that τ l becomes a corrected value for a l , towards which a l could move slowly and slightly (not too fast, to make sure the sequence of τ l 's from the top τ L has time to converge).
Proposition 6. The first variant of the feedback weight update ∆ 1 Ω l performs an SGD step with respect to reconstruction error 1 2 ||a l−1 − g l (f l (a l−1 ))|| 2 while γ = 1.
Proposition 7. When u l has converged, the first variant of the feedback weight update ∆ 1 Ω l performs an SGD step with respect to the inversion error 1 2 ||f −1 l (τ l ) − g l (u l )|| 2 while γ < 1, i.e., the feedback weight change is trying to hasten the convergence of g l (u l ) to a value such that g l (u l ) = f −1 l (τ l ). When γ → 0, ∆ 1 Ω l → 0 as well.
We note that for γ varying from 1 to 0, τ l−1 varies linearly from a l−1 to g l (u l ), making the update interpolate linearly between minimization of the reconstruction error and helping to hasten the convergence of g l (u l ) to make g l (u l ) achieve the optimal target propagation f −1 l (τ l ). Proposition 8. For any γ and whether or not u l has converged, Proof.
Hence the two updates are equivalent when γ = 1, and whereas ∆ 1 goes to 0 when γ → 0, it is not the case for ∆ 2 . The additional pressure Ω l receives under ∆ 2 can thus best be understood by considering γ → 0. In that case, the feedback units are pushed by ∆ 2 towards a l−1 and not just towards f −1 (τ l ). This could be seen as a regularizer, especially when there are many inverses: we would prefer to pick among all the inverses ones which are closer to the feedforward activations (i.e. require the least changes in the weights). Another more farfetched consideration is that in the case where the loss can be interpreted like an energy function and the feedback path as a stochastic generative model (if noise is injected in its computation), then we want that generative path to match the feedforward path, as in the Wake-Sleep algorithm [Hinton et al., 1995].

Conclusion and Future Work
This is only a proposal, with many uncertain avenues left to explore and the need to validate these ideas through simulations and confrontation with biological knowledge (and possibly biological experiments) and machine learning experiments. Many variations on the themes introduced here are possible and will need simulations to disentangle experimentally and understand better. This proposal opens the door to a style of biologically plausible analogues to backprop which is a variant of target propagation while sharing some properties of error propagation schemes (the targets being close to the activations) and can potentially offer optimization advantages compared to SGD with backpropagation, because of the Gauss-Newton approximation implicitly being computed (without having to represent and invert high-dimensional matrices).
One generalization which would be interesting to consider is one in which lower layers would not just try to minimize the output error but also make the inversion easier for upper layers, which could be achieved by extending the kind of feedback loop we analyzed here between consecutive layers to a global feedback look. The general idea is that we are trying to make all latent layers move away from the current feedforward activations in a direction which provides less output error but also a better match between the "prior" top-down information and the "data-driven" bottom-up signal. An important principle exploited in this paper is that we only want to move slightly away from the datadriven inference and towards one which is more globally coherent with other signals (like the target), in order to obtain an update training signal for all the weights. Like in the Wake-Sleep algorithm, we want to make the generative part of the model (the feedback loops) slightly more consistent with the data, and we want to data-driven inference to be slightly more consistent with the global loss.
An interesting longer-term question is whether this principle could also be used to implement a limited analogue to backprop through time. If the brain has a memory of key past events, as well as predictive computations to link them (by a form of association) both forward and backward in time, then all the ingredients are present to achieve credit assignment through time, in a style similar to Sparse Attentive Backtracking [Ke et al., 2018, Kerg et al., 2020. In the spirit of the Consciousness Prior [Bengio, 2017], these credit assignment calculations through time may need not be for all the units in some circuit: it would be enough that the forward and backward associations involve only a few relevant high-level semantic variables. The learned back-through-time target propagation could use attention mechanisms to choose aspects of past moments (events) which are relevant and need to be "fixed" to achieve some current targets. Attention could also be used in a top-down fashion when calculating targets, to focus the target change in the most relevant neurons. This is especially relevant when the lower layer for which a target is computed has more neurons than its upper layer.