A Theoretical Framework for Inference Learning

Backpropagation (BP) is the most successful and widely used algorithm in deep learning. However, the computations required by BP are challenging to reconcile with known neurobiology. This difficulty has stimulated interest in more biologically plausible alternatives to BP. One such algorithm is the inference learning algorithm (IL). IL has close connections to neurobiological models of cortical function and has achieved equal performance to BP on supervised learning and auto-associative tasks. In contrast to BP, however, the mathematical foundations of IL are not well-understood. Here, we develop a novel theoretical framework for IL. Our main result is that IL closely approximates an optimization method known as implicit stochastic gradient descent (implicit SGD), which is distinct from the explicit SGD implemented by BP. Our results further show how the standard implementation of IL can be altered to better approximate implicit SGD. Our novel implementation considerably improves the stability of IL across learning rates, which is consistent with our theory, as a key property of implicit SGD is its stability. We provide extensive simulation results that further support our theoretical interpretations and also demonstrate IL achieves quicker convergence when trained with small mini-batches while matching the performance of BP for large mini-batches.

Introduction
Backpropagation (BP) [38], the most successful and ubiquitously used learning algorithm in deep learning, has long been criticized as being incompatible with known neurobiology [11,24]. This criticism has inspired researchers to search for alternatives to BP that may better fit biological constraints. One alternative to BP is the inference learning algorithm (IL) [37, 49,43]. IL is an energy based algorithm [42,50], which first minimizes an energy function, known as free energy, w.r.t. neuron activities, then at convergence updates weights to further minimize energy. Interest in IL has grown in recent years because IL arguably better fits biological constraints than BP, while matching BP's performance on many tasks. For example, BP is sometimes said to compute weight updates using information non-local to the synapse, whereas synapses in the brain depend on local information having to do with the activity of pre-synaptic and post-synaptic neurons [11,24]. IL, on the other hand, uses local learning rules that are similar to Hebbian models of synaptic plasticity in the brain (e.g., [49,50]). Another criticism is that error feedback propagated by BP does not alter feed forward (FF) activities, which is hard to reconcile with the highly interconnected circuitry of FF and feedback signals in the brain [23,26]. Feedback used in IL, on the other hand, necessarily alters FF activity in order to minimize free energy. Further, IL has long been used to train predictive coding (PC) models of circuits in the cortex (e.g. [37,14,44,4,1,8]). PC and IL have even been proposed as canonical computations used throughout the brain [16,15,4]. More recently, IL has achieved comparable accuracy to BP on supervised learning tasks with small and medium scale data-sets (e.g, [49,2,40], see also below). IL has also been used to train PC networks to perform auto-associative memory tasks, and on some tests has achieved a higher recall accuracy than both modern Hopfield networks [18,36] and denoising auto-encoders trained with BP [41,40].
IL's biological connections motivates further analysis into its optimization properties given that such an analysis may lead to a better understanding of the principles behind synaptic plasticity and optimization in the brain. Additionally, IL's competitive performance with BP further motivates an investigation into the question of why IL is competitive with BP. Similar questions, such as How does IL behave differently than BP? What advantages/disadvantages does IL have as compared to BP? can also aid in understanding which applications, if any, are better suited for IL rather than BP and vice versa. The optimization properties of IL and their relations to BP, however, are still not fully understood. This gap in our understanding is the main motivation for this paper.
In this paper, we develop a novel theoretical framework for IL that describes its optimization properties and their similarities and differences from BP. More specifically, 1) we expand beyond current literature by deriving a general form of the IL algorithm, which we call Generalized IL (G-IL). 2) We demonstrate our main result, which is that G-IL closely approximates an optimization method known as implicit stochastic gradient descent (implicit SGD). Implicit SGD is distinct from explicit SGD (the standard form of SGD used in machine learning and the sort executed by BP). 3) We provide theoretical guarantees concerning IL's stability and a connection to Gauss-Newton optimization and BP. 4) We identify a learning rule and variable settings needed in order for IL to equal implicit SGD and use this result to develop a novel implementation of IL called IL-prox. 5) Finally, we present extensive simulation results that support our theoretical findings, and for the first time show certain performance advantages of IL over BP, such as better stability across learning rates and improved convergence speed with biologically realistic streaming data and small mini-batches. These results collectively suggest that, not only does IL better fit biological constraints but is mathematically justified and has performance advantages over BP in more biologically realistic training scenarios.

Related Works
Whittington and Bogacz [49] proved that parameter updates performed by the IL algorithm approach those of BP as the global loss approaches zero and optimized activities approach feed-forward activities. However, this same proof implies that non-zero IL updates essentially never equal those of BP. Additionally, this analysis leaves open the question of how IL is able to minimize the loss in a stable manner early in training when the loss is large and IL updates diverge significantly from those of BP. Other work has shown IL is formally related to variational expectation maximization [29,25], which is a learning algorithm with convergence guarantees [31]. This analysis still leaves open the question of how IL relates to BP and the broader framework of gradient-based methods that are the backbone of deep learning. The focus of the current work is to develop insights into these open questions. More recently, several works have altered IL to more closely approximate SGD and BP [30,43]. These variants significantly change the original IL algorithm, and some the alterations seem hard to reconcile with neurobiology [30,43]. Alternatively, instead of altering IL to better fit SGD and BP, we show the standard implementation of IL already closely approximates a working optimization method known as implicit SGD, and we show the variable settings under which the standard implementation of IL better approximates implicit SGD.
Algorithms similar to IL have been proposed. For example, the alternating minimization algorithm [10], which was developed independently of the predictive coding and IL framework from neuroscience, updates weights to minimize local prediction errors similar to IL (see equation 3 below). However, [10] do not connect their algorithm to implicit SGD. A more recent variant of this algorithm by Qiao et al. [35] show some formal links to proximal operators, an optimization process related to implicit SGD [34] (see appendix B.1). However, they do not interpret the weight updates as performing the proximal update or implicit SGD as we do here with IL. Proximal operators have also recently been incorporated into learning algorithms for deep networks. Frerix et al. [13] developed a BP, proximal algorithm hybrid. Lau et al. [20] developed an algorithm that combined and algorithm known as block coordinate descent with the proximal operator. Other works use proximal operators to reduce noise or perform other tasks, then BP is used to to update at convergence [27,52]. All of these algorithms, however, are distinct from IL, as they use BP or other non-local gradient information to compute weight updates or local targets.
Finally, IL has similarities to target propagation algorithms, in which local target activities are computed and used to compute local error gradients to update weights. Target propagation [5] and difference target propagation [22] utilize approximate Gauss-Newton updates on activities to create local targets [28,6], while other algorithms compute targets using gradient updates on activities [21,33,3,19]. These algorithms have important differences from IL, e.g., they use a different learning rules than IL, which we discuss further below. In sum, previous works have not developed the same mathematical interpretation of IL that we do here and have not developed the same mathematical descriptions of the differences between IL and BP. Further, none of the works mentioned above produced the same empirical findings we do in our simulations below.

Notation
Term Description W n Weight Matrix, pre-synaptic layer n h n Feedforward Activity layer n, h n = W n−1 f (h n−1 ) h n Optimized/Target Activity layer n ∆h n Optimized Activity Change layer n, ∆h n =ĥ n − h n p n Local Prediction layer n, p n = W n−1 f (ĥ n−1 ) e n Local Error, layer n, either e n =ĥ n − h n or e n =ĥ n − p n Table 1: Notation Notation describing a multi-layered FF network (MLP) is summarized in table 1. To simplify notation, we assume bias is stored in an extra column on each weight matrix W n and we assume a 1 is concatenated to the end of the pre-synaptic activity vector (h n orĥ n depending). All the following results should hold with and without biases.

Generalized Backpropagation
Consider a set of parameters θ (b) at training iteration b with global loss function L. An explicit SGD iteration subtracts the gradient of the loss from the parameters: θ (b+1) = θ (b) − α ∂L ∂θ (b) , with step size α. BP [39] is a common algorithm for implementing SGD in deep networks. Here, to allow for easier comparison to IL, we consider an alternative implementation of SGD in deep networks that involves computing local target activities at hidden and output layers, similar to [21], [3], [33]. Specifically, the local target at layer n,ĥ n , is computed as followŝ where local loss l n+1 = 1 2 ĥ n+1 − h n+1 2 = 1 2 e n+1 2 . The target at the output layerĥ N = h N − ∂L ∂h N . Notice that e n+1 =ĥ n+1 − h n+1 = h n+1 + ∂ln+2 ∂hn+1 − h n+1 = ∂ln+2 ∂hn+1 = − ∂L ∂hn+1 . As noted in the above equation, this local target update is equivalent to subtracting the global loss gradient from the FF activity h n (for proof see appendix [3]). To update weight W n one uses the gradient of the local loss: This update uses only local information and is equivalent to subtracting the global loss gradient from the weights (see [3]). This local target formulation of SGD is similar to other target propagation (TP) [5,22] algorithms, which use a similar learning rule but compute targets slightly differently. For example, a well-known variant of TP is difference target propagation (DTP) [22], which computes local targets by approximating a Gauss-Newton update on activities rather than gradient updates [28, 6]:ĥ n ≈ h n − J + N,n e N , where J + N,n is the pseudo-inverse of the Jacobian of the forward network from layer n to N and e N = y − h N . Gauss-Newton updates approximate (are within 90°) of SGD updates [28], and recent attempts to improve DTP generally do so by making DTP updates more similar to SGD updates (e.g. [28,12]). Additionally, variants of BP algorithms like random feedback alignment [23] and sign symmetric feedback alignment [51] approximate SGD. Since all of these methods approximate the (explicit) SGD update over weights, we consolidate them all under the heading of generalized BP (G-BP). The target-based formulation of G-BP is shown in algorithm 1.

Generalized Inference Learning
We now present a general description of the IL algorithm, which we call Generalized Inference Learning (G-IL). A summary of the algorithm is presented in algorithm 2. At each layer n of an MLP, G-IL stores two variables:ĥ n and p n . At initialization these variables are set to feed-forward activity values:ĥ n = p n = h n . Thenĥ n are altered over time to minimize the energy function F : where L is the global loss, 1 2 ĥ n − p n 2 are local losses, γ decay n 1 2 ĥ n 2 is an optional regularization term, and γ are positive scalar weighting terms. The prediction p n is computed as p n = W n−1 f (ĥ n−1 ). (Note p n may also be computed as p n = f (W n−1ĥn−1 ). Our theoretical results use the first formulation, but we find little difference in performance between the two in practice. See simulations below.) G-IL first minimizes F w.r.t. to activitiesĥ (inference phase). While optimizingĥ, predictions also change since p n+1 = W n f (ĥ n ). After argminĥ F is approximated, F is again minimized, now w.r.t. the weight, i.e. argmin W F . Typical IL algorithms use SGD to perform the inference phase and either the LMS rule (below) or Adam optimizers to update weights (e.g. [37,49]). G-IL, on the other hand, is general w.r.t. the optimization processes that approximate argmin W F and argminĥ F . We consider two weight updates in our analysis and simulations below. First, is a local error gradient update, which is equivalent to the least-mean squares (LMS) rule: where e n+1 =ĥ n+1 − p n+1 and local loss l n+1 = 1 2 e n+1 2 . The LMS rule does not solve argmin W F but only approximates it. The next rule, which is equivalent to the normalized least mean squared rule (NLMS) [32] is (5) The NLMS rule performs a gradient update but with an inverse squared 2 norm of the pre-synaptic activities as its learning rate instead of a hyper-parameter. The NLMS rule is the minimum-norm solution to argmin W F in the case where we use mini-batch size 1 (see proposition B.1). The LMS rule closely approximates the NLMS update in the sense the two are proportional. Comparing these rules to the update rule for G-BP (equation 2) we see the main difference is the pre-synaptic term used in the rule: BP uses FF activity h n , while IL uses optimized/target activityĥ n . This difference leads to significantly different stability properties between the algorithms, as we note below. Target activitiesĥ are initialized to FF activities h then updated to minimize energy F . Afterward, weights are updated to further minimize F . Weight updates 'align' withĥ activities soĥ become FF activities given the same data point.ĥ are activities that will both improve the loss L and minimize ∆θ 2 . This process is equivalent to minimizing the proximal loss and approximates implicit SGD.

Implicit Gradient Descent
Let the set of MLP weight parameters where b is the current training iteration. The explicit SGD update uses the loss gradient produced by the current parameters: where α is step size and L is the loss measure. This update is explicit because the gradient can be readily computed given (θ (b) , x (b) , y (b) ). The implicit SGD update uses the loss gradient of the parameters at the next training iteration: This is an implicit update because θ (b+1) shows up on both sides of the equation. Unlike explicit SGD, θ (b+1) cannot be readily computed using available quantities. However, the implicit gradient is equivalent to the solution of the following optimization problem (see appendix equation 13): This update is equivalent to the proximal operator/update [34]. This proximal update changes parameters θ (b) in a way that minimizes the loss L and the magnitude of the update, which helps keep θ (b+1) in the proximity of θ (b) . (For more details on the proximal update and its relation to implicit SGD see appendix B.1).

Main Results
We now present our main result, which is to show G-IL is equivalent to implicit SGD in certain limits that are approximated well in practice. We focus on the case where single data-points (mini-batch size 1) are presented each training iteration. This scenario is more biologically realistic than mini-batching and is the case where IL best approximates implicit SGD. We discuss the case of mini-batching in appendix B.5. We first define a kind of IL algorithm that minimizes the proximal loss (equation 8) w.r.t. neuron activities and show G-IL approaches this algorithm in certain limits.
Here b is the current training iteration and θ * are the parameters updated with the NLMS rule. IL-prox is the same as G-IL except during the inference phase, it minimizes the proximal loss w.r.t. activitieŝ h instead of the energy F . Intuitively, minimizing the proximal loss w.r.t. activities will result inĥ that yields a small loss after weights are updated and small weight update norms. Weight update norm 1 2 θ * − θ (b) 2 is known during the inference phase because the NLMS learning rule is used and known explicitly: where α n = f (ĥ n ) −2 . Thus, it can be minimized w.r.t.ĥ n . As we explain in the appendix (see lemma B.1 and proposition B.1), the loss L(θ * ) is also explicitly known during the inference phase and can be optimized w.r.t.ĥ n . Now, let α be a 'global' learning rate hyper-parameter, and α n be layer-wise learning rates used in the actual weight update ∆W n .
Theorem 4.1. Let α n = f (ĥ n ) −2 and assume mini-batch size 1. In the limit where γ decay Hence, in these limits G-IL is equivalent to IL-prox.
The proof can be found in supplementary materials theorem B.1. This theorem states that when using the NLMS rule, G-IL is increasingly similar to IL-prox, as the γ weighting terms in the free energy (equation 3) approach the scalar values specified in the above theorem. In these limits, argminĥ prox = argminĥ F and G-IL is equivalent to IL-prox. The intuitive explanation of this theorem goes as follows: minimizing the proximal loss (equation 11) requires minimizing L w.r.t. θ * while also minimizing the update norm 1 2α θ * − θ (b) 2 . The inference phase of IL minimizes L w.r.t. θ * by initializing output layer activityĥ N to FF output h N then shifting it toward global target y. Upon a weight update this will yield a θ * that produces a smaller loss L. See figure 4.1 for visualization. The magnitude of a weight update 1 2α αe n+1 f (ĥ n ) T 2 intuitively depends on the magnitude of local prediction error e n+1 . F is a sum of the magnitudes of errors (see equation 3), so minimizing F minimizes the magnitudes of errors and consequently weight update norms 1 2α θ * − θ (b) 2 . Importantly, the γ settings under which theorem 4.1 holds can be computed exactly in practice and are of reasonable magnitude that is easily approximated in practice (e.g., they do not approach ∞ or 0 and typically approach positive scalars). This is opposed to the limits under which IL approximates BP/explicit SGD, which is the limit whereĥ n → h n and thus ∆W n → 0 [49], a limit that does not typically obtain in practice.
IL−prox be the parameters updated by IL-prox (see def. 4.1) and θ This theorem states that the IL-prox update is equivalent to the proximal update (equation 8), and thus implicit SGD. It follows that, under the γ settings in theorem 4.1 the G-IL update is also equivalent to the proximal update/implicit gradient update.
Further theoretical results are developed in the appendix. A novel closed form description of IL targets is developed C.5. The closed form description shows local targets approximate a regularized Gauss-Newton update on FF activities (see section C). We call this closed form approximation IL-GN. We study the stability of IL-GN when using the LMS learning rule, and compare to an analogous BP network, which uses Gauss-Newton updates (BP-GN). We show that IL-GN weight updates push output layer activities down their loss gradients for any positive learning rate, i.e., for any positive learning rateĥ where j is a positive scalar (theorem D.1). We show the same is not true of BP-GN, which only has this property for a small range of learning rates (theorem D.2). We develop a partial explanation of why IL is more stable than BP that is based on differences in their learning rule (for details see section E). These results make the important point that the way IL computes and uses local targets in its weight update significantly contributes to its stability, while the way G-BP computes and uses local targets in its weight update (which is distinct from IL, see section 2) is largely responsible for its instability. The IL learning rule and method of computingĥ is advantageous over that of G-BP in this sense.   In this section, we test and compare BP based algorithms to several IL algorithms. BP algorithms compute and use explicit gradients to update weights, while IL algorithms compute and use approximate implicit gradients to perform updates. We test BP models that either perform simple gradient descent (BP-SGD) or use Adam optimizers (BP-Adam). IL-SGD is a simple implementation of IL that is similar to implementations of [37,49]. It uses the LMS rule (local error gradient update) to update weights and optimizesĥ using SGD to near convergence (25 update steps). IL-Prox is our novel algorithm that is based on definition 4.1. IL-prox uses the NLMS rule to update weights. As theorem 4.1 shows, to minimize the proximal loss, the learning rate is only used to determine how much the output layer activitiesĥ N are pushed toward the target y. As α → ∞,ĥ N → y and as α → 0,ĥ N → h N . IL-Prox optimizesĥ using SGD to near convergence (25 update steps). IL-Prox Fast is the same as the IL-Prox algorithm except it truncates the optimization ofĥ to only 12 iterations. IL-Adam computes the weight gradients using IL-SGD, then uses Adam optimizers to update weights. IL-prox Adam computes the normalized weight gradient using IL-prox and uses Adam optimizers to update weights. All simulations average results over at least 5 seeds of each model type. Best test scores from these averages are shown in tables. All models use ReLU activations at hidden layers. Softmax is used at output layer for classification, while sigmoid is used on the autoencoder task.  Table 3: Best test reconstruction loss after one epoch of training on with CIFAR-10 images, minibatch size 1. BCE was summed across pixels and averaged across data-points. Fully connected network layer sizes 3072-1024-512-100-512-1024-3072. Top two scores highlighted in bold.

Experiments
Supervised and Self-Supervised Learning We begin by training these models on supervised classification tasks. Results for CIFAR-10 are summarized in table 2. Consistent with previous results [49,2], in the case where large mini-batches are used, IL algorithms achieve near equal accuracy as BP-SGD and BP-Adam when using Adam optimizers. But IL-SGD and IL-prox algorithms, which do not use Adam optimizers, do not perform as well as BP algorithms in this scenario. However, in the case where they are trained with mini-batch size 1, IL algorithms improve accuracy significantly faster than BP-SGD for the entire first epoch (50,000 iterations) of training. IL-SGD and IL-prox algorithms perform similar to and even slightly better than BP-Adam, despite not using momentum or stored, parameter-wise adaptive learning rates like Adam. As far as we know, this is a novel result. Previous works have not tested IL in case where small mini-batches are used (e.g., [49,2]). To further test this improvement in optimization speed, we also trained these models on MNIST and fashion-MNIST classification tasks (table 6), and we trained autoencoders on CIFAR-10 images (table 3), with mini-batch size 1. Slight speed-ups in training were seen on MNIST data sets (see figure F.2). On the autoencoder task with CIFAR-10, we found again IL-prox algorithms improve loss as quickly and even slightly quicker than BP-Adam (see figure 2).
Stability Test Implicit SGD is unconditionally stable, which, roughly, means that it is able to discount the loss in a stable manner for nearly any positive learning rate [46]. We compare performance of BP-SGD, IL-SGD, and IL-prox on MNIST (table 5) and CIFAR-10 (table 4) classification task across different learning rates. As a control, we also test a hybrid of BP and IL-prox (BP-prox). BP-prox, like IL-prox, uses the learning rate only to adjust the target at the output layer and uses the NLMS rule to update weights. Unlike IL-prox, BP-prox computes local targets using the G-BP algorithm rather than the IL algorithm. This control helps to show the stability of IL-prox is not only due to normalizing the weight updates and softly clamping the output layer. Rather, the way IL computes and uses local targets in its learning rule also contributes to stability. Results can be found in table 5 and 4. Blank entries are those with accuracies below 12%. IL-SGD is more stable than BP-SGD, and IL-prox is more stable than BP-prox. These results provide evidence that the way IL computes and uses targets contributes to stability. Additionally, the prox algorithms are more stable than SGD algorithms, showing that using the normalized gradient to update weights and α to adjust the output layer target rather than scale weight updates also contributes to stability.

Further Results
We computed the update norms ∆θ during training of each model on the CIFAR-10 classification task. Average magnitudes and standard deviations are shown in appendix F. IL algorithms had similar or smaller update magnitudes than BP-SGD on most of their training runs. This suggests IL algorithms do not minimize loss more quickly because their updates are greater in magnitude than BP-SGD. One possibility is that IL updates take a more direct path to local minima than the path of steepest descent. This is consistent with other experiments in which we analyzed the output layer of BP-SGD and IL-SGD. We measured how the FF output layer values h N changed after updates. IL updates tended to change the output layer values more and in a direction nearer the minimum norm (shortest) path to that target y than BP-SGD updates (figures F.4, 6). . For both of these initialization, the proximal quantity is reduced during inference to near zero. We next initializedĥ randomly (Rand Init. and Rand + Full Clamp). Both networks reduce the proximal loss, but nowhere near the global minima, which may help explain why the network needs to be initialized to FF activities when F is minimized with SGD. These results generally support theorem 4.1, which states minimizing F minimizes the proximal loss.

Limitations
Our theoretical results show that IL algorithms closely approximate implicit SGD in the case where only a single data-point (i.e., mini-batch size 1) is presented each training iteration. This approximation is looser in the case where large mini-batches are used (see appendix B.5). This is no problem from a biological point of view, since the brain does not train with large mini-batches. Additionally, when Adam optimizers are used IL still performs similarly to BP. However, it may still be desirable to further study the optimization properties of IL with mini-batches to improve performance on machine learning tasks where mini-batches are typically used.

Discussion
IL was originally used to train predictive coding models of cortical circuits [37]. The success of subsequent predictive coding models lead some neuroscientists to propose the hypothesis that predictive coding and IL are canonical computations used throughout the brain (e.g., [14,15,4,17]).
In this paper, we provided novel mathematical justification for IL by showing it closely approximates a proper optimization method known as implicit SGD, which is distinct from the explicit SGD executed by BP. This finding, in conjunction with the theory IL and predictive coding are canonical computations in the brain, suggests the interesting possibility that implicit SGD, and not explicit SGD, is the central optimization method utilized by biological neural circuits. We observed two performance advantages that fit with this possibility. First, consistent with our theoretical interpretation, IL algorithms are more stable across learning rates than BP. This property is compatible with the fact that synaptic plasticity in the brain can fluctuate rapidly, e.g. due to neuromodulation, yet maintains stable learning properties. The stability of IL is arguably its main advantage over BP, which is consistent with claims that the main advantage of implicit SGD over standard/explicit SGD is its stability across learning rates (e.g., see [46,48]). Second, we found IL minimized the loss more quickly than BP in the more biologically realistic scenario where single data points, rather than large mini-batches, are presented each training iteration. These simulations provide evidence that IL is more compatible with biologically realistic learning scenarios than BP. Collectively, our findings suggest IL is a promising basis for developing biologically constrained, high performing, local learning algorithms.

Broader Impacts
Our work is largely theoretical so it does not suggest any direct, possibly negative, societal impacts, as far as we can see.

A Preliminaries
The theoretical analyses developed in this appendix focuses on the learning scenario where single data-points are presented to an MLP each training iteration (i.e., mini-batch size 1). We focus on this scenario in large part because it is more biologically realistic than training with large mini-batches, and a central goal of the paper is to analyze IL in order to better understand how optimization and credit assignment may work in the brain. Other analyses of learning algorithms make similar simplifying assumptions (e.g., see [28]). We discuss the case of mini-batches briefly in section B.5, and our empirical simulations in this paper and previous works in the literature show that IL algorithms are able to perform competitively with BP when trained on mini-batches (e.g., see [49,2,40] as well as in the scenario where single data-points are presented. This suggests using IL with large mini-batches is approximated by the case of training with single data points. We leave it to future work to analyze differences in the behavior of IL in the case of small versus large mini-batches in more detail. Our analyses considers the version of IL that uses local predictions computed p n+1 = W n f (ĥ n ) at hidden layers and p 1 = W 0ĥ0 at the input layer. This kind of prediction can be interpreted as a prediction of neuron pre-activations. This is slightly different from predictions of neuron activations: Although we only analyse the pre-activation version of the prediction, we find the two behave very similarly in practice (e.g., our IL-SGD model uses prediction of activations, while our IL-prox models use predictions of pre-activations and both achieve similar performance on a variety of tasks). Thus, our analyses here should apply approximately to IL models that use the activation prediction equation.
Finally, we will sometimes use two slightly different expressions to describe the loss produced by MLP parameter θ. First, the form L(θ (b) , x (b) , y (b) ) refers to the loss produced by MLP parameters θ (b) at training iteration b, given data point x (b) and prediction target y (b) . This loss generally measures the difference between the FF activity at the output layer h N and target y. Second, the form L(ĥ N , y (b) ) refers to the loss between the output layer optimized activitiesĥ N and target ) the loss between output layer FF activity and the target. Note that both h N andĥ N are initialized to W N −1 f (h N −1 ). We generally assume some non-linearity is applied to these activities at the output layer before being inputted into loss. For example, in the case where cross entropy loss is used, one would compute L(σ(h where σ may be sigmoid in the case binary cross entropy or softmax in the case a categorical cross entropy loss. To simplify notation, we do not write this σ below, but assume it is 'built' into the loss as needed. None of the proofs below are affected by this assumption. B G-IL and Implicit SGD B.1 A Brief Introduction to Implicit SGD Let θ (b) be a set of parameters, x (b) input data, and global output target is y (b) , where b is the current training iteration. The explicit SGD update uses the loss gradient produced by the current parameters: where α is step size and L is the function being minimized. This update is explicit because the gradient can be readily computed given The implicit SGD update uses the loss gradient of the parameters at the next training iteration: This is an implicit update because θ (b+1) shows up on both sides of the equation. Unlike explicit SGD, θ (b+1) cannot be readily computed using available quantities. One can still perform an implicit gradient update, however, by computing the solution of the following optimization problem: This update is known as the proximal algorithm or proximal point method [34]. The proximal algorithm turns the general optimization problem of minimizing L across all the data set into a series of sub-optimization problems. This algorithm is a specific application of the more general proximal operator: where f is a function with scalar output. We can see the proximal algorithm changes parameters θ (b) in a way that minimizes the loss L and the magnitude of the update, which helps keep θ (b+1) in the proximity of θ (b) . We can also see the proximal algorithm only requires using the known quantities Equation 11 is closely related to incremental proximal methods [7], which is a stochastic version of the standard proximal algorithm, which usually performs batch updates. The equivalence of the proximal algorithm in the stochastic setting to the implicit SGD update (equation 10) can be shown easily using the fact that at the minimum θ = θ (b+1) and the gradient of 11 is zero: Explicit and implicit SGD have similar convergence guarantees. Analysis of the convergence properties of implicit SGD in linear regression models was done by [46,47], in generalized linear models by [45]. However, explicit and implicit SGD have importantly distinct stability properties. It is well-known that explicit SGD is highly sensitive to learning rate. Implicit SGD, on the other hand, is unconditionally stable in the sense that equation 11 monotonically decreases the loss function ). One can also gain an intuition for the unconditional stability of implicit SGD by noting the differences in the way the learning rate α is used in the explicit SGD update (equation 9) versus the proximal algorithm/implicit update (equation 11). In the explicit SGD update α is used to scale the gradient, which clearly implies that as α → ∞ so does the magnitude of the update. However, with the proximal algorithm α only controls the relative weighting of the two terms in equation 11. As α → ∞ we see the regularization term 1 2α θ − θ (b) 2 is down-weighted to 0. In this case, the updated parameters are simply those that best minimize the loss function L, which lead to an update that reduces the loss given x (b) . For further details on the stability of explicit and implicit SGD in linear models see [46,45,47].

B.2 G-IL Local Targets as Target FF Activities
In this section, we show that when a single datapoint is presented each training iteration, there are two interesting properties concerning the local targets computed by IL and the NLMS rule used by IL. We express these properties in lemma B.1 and proposition B.1. This lemma and proposition are used in proofs below.
First, we show that IL local targetsĥ at the next training iteration given the same data point when weights are updated to fully minimize local prediction errors. In this sense, IL local targets are target FF activities.
. Let the solution parameters at iteration b be θ * = argmin θ F , where F is the energy function defined in equation 3. In the IL algorithm, minimizing F w.r.t. weight W n equates to updating weights in a way that minimizes the local prediction error e n+1 = 1 2 ĥ n+1 − p n+1 2 (see section 3.3). In the case where the MLP is only presented with a single data-point each iteration, this local learning problem can be solved such that there is zero error. More specifically, let W * n be the weight matrix at layer n of θ * . When training with single data-points, this solution weight matrix has the propertyĥ n ) (for an example of one solution with this property see proposition B.1). It follows trivially thatĥ n are equivalent to the FF activities of θ * when given the same data-point x (b) as input: Lemma B.1. Consider a non-linear MLP trained with G-IL at iteration b with mini-batch size 1. Let h * n be the FF activities of solution parameters θ * at layer n given data-point Proof. At initialization,ĥ Next, by definition at the input layerĥ This property holds for remaining layers because by definitionĥ This same procedure can then be repeatedly applied to show the same is true of all remaining layers n.
The solution matrices W * n are non-unique in the case where single-data points are used. However, there are unique minimum-norm solution matrices, which can be computed using the NLMS rule Proof. Here we follow the proof of [9], which proved a similar result concerning linear regression models trained with a variant of NLMS.
Using the method of Lagrangian multipliers we can rewrite argmin Wn where λ is the Lagrangian multiplier and e (b+1) n ). First we compute the gradient of P w.r.t. the weights and set to 0, The partial w.r.t. the Lagrangian is then Rearranging 15 we get Finally, substituting the value for λ into 17 we get which is exactly equal to the NLMS rule used in the IL algorithm (equation 5).
In sum, when training with single data-points, activitiesĥ optimized by IL are target FF activities in the sense that they become FF activities after F is minimized completely w.r.t. weights and local errors go to zero. The NLMS rule yields such a solution. One important implication of this is that the LMS update (equation 4), which is commonly used in practice, is a good approximation of the NLMS solution since the two are proportional.
This lemma and proposition allow us to do something important, which lays the basis for the connection between implicit SGD and G-IL. In particular, this lemma and proposition show us that the network solution θ * , its FF activities h * n , and ∆θ can be expressed implicitly in terms ofĥ. For example, lemma B.1 shows that h * n =ĥ n . This implies that we know what the loss produced by θ * : where L(ĥ N , y (b) ) is the loss measure, which in the case of MLPs is some measure of the difference between prediction target y (b) and FF output layer value, which for θ * is equivalent toĥ N (as implied by lemma B.1).
Additionally, proposition B.1 shows that each ∆W n can be expressed as f (ĥ n ) T , where e n+1 again is can be expressed in terms of optimized activities: e n+1 =ĥ This means This allows for the possibility to minimize the proximal quantity w.r.t.ĥ n . That is, it allows for an IL algorithm, where during the inference phase the algorithm approximates argminĥ(L(θ * ) + 1 2α θ * − θ (b) 2 ), where L(θ * ) and θ * − θ (b) 2 are defined implicitly in terms ofĥ using equations 20 and 21. At each step in the inference phaseĥ are updated to minimize the proximal loss, which results in a change in θ * , since θ * is defined in terms of theĥ. After a minimum is reached weights are updated with the NLMS rule such thatĥ values become the FF values given the same data point. The inference phase, in other words, finds the FF values of a θ * that best minimizes the proximal loss, then uses those FF values to update the actual parameters θ (b) such that θ (b+1) = θ * (see figure 4.1 for visualization).
Let's define such an algorithm as proximal inference learning (IL-prox), since it, like G-IL, computeŝ h n by minimizing an objective w.r.t. activities, then updates weights afterward: Definition B.1. Proximal Inference Learning (IL-prox). An algorithm identical to G-IL (algorithm 2), except activities are optimized to approximate argminĥ n (L(θ * ) + 1 2α θ * − θ (b) 2 ) and weights are updated using the NLMS rule.
Again, the θ * here is defined implicitly in terms of theĥ values using equations 20 and 21. In the next, section we show that under certain variable settings it is the case that argminĥ F = argminĥ(L(θ * ) + 1 2α θ * − θ (b) 2 ).

B.3 Generalized Inference Learning Approximates Proximal Inference Learning
Let α be a 'global' learning rate hyper-parameter for θ, and α n be the layer-wise learning rate used in the weight update ∆W n . Finally, let θ * be defined as the parameters after the NLMS update is applied to each weight matrix, as in the last section.
Hence, in these limits, G-IL is equivalent to IL-prox.
Gradient of F In IL networks, activities are updated with local gradients of F , as follows: At the output layer, ∂F (ĥ N ) and at hidden layers: where f (ĥ n ) = ∂f ∂ĥn . We see the gradient of the loss L term is only taken w.r.t. the local output layer activityĥ N , and the gradient of prediction error at layers n and n + 1 are taken w.r.t. to local layer n. We assume whatever process is used to minimize F similarly only uses local information such that its minima are described by the above gradients when they equal zero.
Gradient of Prox w.r.t. Output Layer We first compute gradients for the output layer. As with G-IL, we assume only local gradients of the proximal loss are used to update the activities. The gradients of the proximal loss w.r.t. output layer activityĥ N is sinceĥ N is local to L and since W N −1 are the only weights local toĥ N .
The gradient ∂P rox(ĥ N ) ∂ĥ N , can be expressed in terms of prox 1 and prox 2 as follows: ∂P rox(ĥ N ) can be computed using the chain rule. First, where , f is an element-wise non-linearity, and e N =ĥ N − p N . Then using the chain rule we have where ∂prox2 Note that here α N −1 = f (ĥ N −1 ) −2 as noted in the theorem. In the limit where γ N → α N −1 α , this gradient comes to ∂L(ĥ N ,y) ∂ĥ N + γ N e N , which is equivalent to the gradient ∂F  As above, the gradient ∂prox2(ĥn) ∂ĥn is computed using the chain rule.
where p n = W n−1 f (ĥ n−1 ) and e n =ĥ n − p n . Using the chain rule where ∂prox2 ∂∆Wn−1 = αn−1 α e n f (ĥ n−1 ) T and ∂∆Wn−1 ∂en ∂en ∂ĥn = α n−1 f (ĥ n−1 ). Multiplying together we get ∂prox 2 (ĥ n ) Now we derive prox1(ĥn) hn : where p n+1 = W n f (ĥ n ) and e n+1 =ĥ n+1 − p n+1 . Using the chain rule Unlike the previous gradients, we now need propagate the gradient through the learning rate, which is what is done by the term ∂prox1 Now we substitute our ∂prox1 ∂ĥn and ∂prox2 ∂ĥn terms back into the gradient of P rox(ĥ n ): Finally, we set this gradient equal to 0 (i.e. at a minimum of Prox) and in the limit where γ decay n → e n+1 2 (1 + −2 α 6 n ) and γ n → α −1 n−1 for all n < N . To simplify notation, we just label this limit as lim When set to zero, α 2 n α cancels. The result is exactly equal to ∂F ∂ĥn (equation 23). It follows that, in these limits, for any hidden layer n it is the case that ∂F ∂ĥn = 0 ⇐⇒ ∂P rox ∂ĥn = 0. Since the same result holds at the output layer, it follows that argminĥ F = argminĥ L(θ * ) + 1 2α ∆θ 2 . Hence, under these γ settings G-IL is equivalent to IL-prox.

B.4 G-IL Approximates Implicit SGD
Here we present the theorem showing the IL-prox, and G-IL with the γ setting noted in the above theorem, are equivalent to the proximal algorithm and thus implicit SGD. An intuitive way to think about the relation between the typical proximal algorithm and IL-prox is that IL-prox does what the typical proximal algorithm would do, but in reverse order. The standard proximal algorithm argmin θ L(θ)+ 1 2α θ−θ (b) 2 minimizes the proximal loss w.r.t. parameters, then the new parameters can be used to compute new FF values given the data-point x (b) . IL prox, on the other hand, first computes the new FF values of the parameters that best minimize the proximal loss given x (b) during the inference phase. Then it updates weights so they become the parameters that best minimize the proximal loss. Since both optimization problems are unconstrained (and thus optimize over the same space of possible parameter values), they yield the same parameter updates in the end.
Proof. The weight update procedure performed by IL-prox can be described as follows: where argmin θ F is computed using the NLMS rule and θ * are the parameters updated with the NLMS rule (equation 15). Lemma B.1 and proposition B.1 imply that theĥ produced by argminĥL(θ * ) + 1 2α θ * − θ (b) 2 ) are the FF activities of θ * . The fact that argmin θ F is computed using the NLMS rule implies that the resulting parameters are the θ * of the optimized activities.ĥ are unbounded, real-valued vectors, which implies that argminĥL(θ * ) + θ * − θ (b) 2 is an unbounded optimization problem and also implies that the possible θ * values are also unbounded (because they can vary over a set of FF values that have unbounded possible values). This implies that argminĥL(θ * ) + θ * − θ (b) 2 outputs the FF values of the parameters θ * that best minimize the proximal loss. Then argmin θ F uses those FF activities to update θ (b) such that it become equal to the θ * that best minimize the proximal loss. This implies θ

B.5 A Note about Mini-batches
Our analysis above focuses on the more biologically realistic scenario where a single data-point is presented each training iteration. In this case, the NLMS rule provides the minimum-norm solution to argmin θ F and the LMS update approximates this solution in the sense it is proportional to the NLMS update. However, in the case where mini-batches of size > 1 are used, the NLMS rule is not a solution to argmin θ F . The gradient update of F w.r.t. W n is where f is an element-wise non-linearity and H n is the mini-batch matrix of neurons activities. Each row of H n is aĥ n for one data-point in the mini-batch.
The solution to the local error minimization problem is found by setting the gradient equal to zero and solving for W n : This update is generally not equal or proportional to the gradient update. Gradient updates are thus poorer approximations of the solution(s) to local least mean squared problem in the case of mini-batches size > 1 than they are in the case with mini-batches size equal to one. This may provide some insight into why the IL models we tested do not show performance advantages over BP-SGD when mini-batches are used as opposed to when single datapoints are used to update weights. Future research could explore alternative learning rules that better approximate the solution above in the mini-batch case.

C A Closed Form Description of G-IL Local Targets
In this section, we derive a closed form description of argminĥF for linear MLPs. In addition to simplifying notation, we focus on linear networks because we aim to use this closed form description to analyze the stability properties of linear networks trained with IL in the next section.

C.1 A Brief Intro to Gauss-Newton Optimization and Ridge Regression
Let θ be a vector of parameters, y a vector of prediction targets, X a data matrix. The linear least squares problem is defined as Its closed form solution is (XX T ) −1 X T y which equals X + y when X has linearly independent columns. X + is the left pseudo-inverse of X. For matrix M the left pseudo-inverse has the property I = M + M . When the regularization term λ θ 2 is added to 40 we get ridge regression whose closed form solution is (XX T + λI) −1 X T y, where I is the identity matrix and λ is a scalar.
The Gauss-Newton method is an iterative method used to solve non-linear least squares problems, which has some similarities to the solutions to linear least squares problem just mentioned. Let e (b) be the residuals of the model at where f is the non-linear model function. The Gauss-Newton (GN) update for θ (b) is where J is the Jacobian of f w.r.t. θ (b) and J + is the left pseudo-inverse of J. Note the similarity to the closed form solution to the linear least squares. The parameters are iteratively updated until convergence or a near convergent state is reached.
It is also common to add a regularization term to the Gauss-Newton update as follows where I is the identity matrix and λ is some scalar. Notice the similarity to the ridge regression solution. The regularization term helps prevent the change to θ from growing too large. Note that as λ → ∞ the update approaches gradient descent:

C.2 A Closed form Description of G-IL Local Targets
Consider a linear MLP trained with G-IL using a squared error global loss. Local targets at hidden layers of such networks are the result of the optimization process that attempts to find where λ = γn γn+1 which represents the relative weighting of the two terms (see equation 3). Here we ignore the optional decay term in equation 3.

Theĥ (b)
n can be expressed in closed form by noting that at the minimum of the expression above h n =ĥ (b) n , and the gradient of the expression equals 0. One can thus compute the gradient, set it equal to zero, and solve forĥ Notice the term on the left (W T n W n + λI n+1 is identical to the closed form solution for ridge regression (see previous section), which is a regularized version of a simple target propagation n+1 . The term on the right (W T n W n + λI) −1 p n ≈ p n , where is a scalar, since (W T n W n + λI) −1 is positive definite and approaches a scalar multiple of I as λ → ∞. Thus, at the minimum h (b) n closely approximates a weighted average between p n and the solution to the regularized squared error at layer n + 1.

C.3 Relation to Gauss-Newton Updates
Proposition C.1. In the limit where λ → 0 ( in equation 44) n+1 , which is a Gauss-Newton update on initial activities h Proof.
More generally, in the case where the influence of the top down error term ĥ n − p (b) n 2 is negligible, G-IL targets converge to Gauss-Newton targets.

C.4 Gauss-Newton G-IL as a Closed Form Approximation of G-IL
In practice, one approximates argminĥ n F by initializing h n =ĥ n = p n then minimizes F using some optimization process. At the beginning of this optimization process the top down error ĥ n − p (b) n 2 is small since it is initialized to zero and grows slowly afterward. Additionally, if optimization is stopped early (which is typically done in practice) this error's effect on activities will remain negligible compared to the larger bottom-up error term. This fact along with proposition C.1 suggest Gauss-Newton updated activities are a good approximation of G-IL local targets. Let's define a network that computesĥ n using Gauss-Newton updates as follows Definition C.1. Gauss-Newton G-IL (IL-GN). IL-GN is a closed-form approximation of G-IL for training MLPs which uses the same weight update as G-IL and computes local targets using a Gauss-Newton update. In linear networks the update isĥ n = h IL-GN approximates G-IL when γ = 1 as shown by proposition 45. A γ < 1, however, better captures the fact that in the true closed form solution (see equation 44),ĥ n is a value in between a regularized GN-target and top-down prediction. Simulations below provide further justification that IL-GN with 0 < γ < 1 is a good approximation of G-IL (see figure F.4.

C.5 Some Properties of IL-GN
There are several lemmas concerning IL-GN that will either be used in subsequent proofs or are relevant for understanding IL generally.
The first lemma says that predictions p n in IL-GN networks are a weighted average of h n andĥ n .
Lemma C.1. After allĥ n are computed, local predictions will lie between local sub-targets and FF activity: p n = (1 − γ)h n + γĥ n .
Proof. First, prediction p n can be described as follows: p n = W n−1 (ĥ n−1 + ∆ĥ n−1 ) = W n−1 (h n−1 + γW + n−1 ∆h n ), which implies It then follows that the error after the targets and predictions are updated, e n , is smaller but proportional to the initial error: Lemma C.2. After allĥ n are computed, local errors e n = (1 − γ)∆h n . Proof.
where the second line is computed using 46. Proof.
where the first and third lines are computed using lemma C.2 and the second line using definition C.1.
These lemmas can be used to show local errors are smaller than but proportional to the global loss GN update w.r.t. hidden layer activities. For example, e n can be expressed as follows: Proof. According to definition C.1,ĥ n = h n + γW + n ∆h n+1 , which implies ∆h n = γW + n ∆h n+1 . Applying this same operation recursively from the output layer backward we get ∆h n = γ N −n (W + n ...W + N −1 )∆h N .
Thus, local errors e n can be seen as scaled down version of the global GN update, which is The scaling is such that the nearer e n is to the input layer the more scaled down the error is.

D Stability of IL versus BP
In this section, we analyze certain stability properties of IL and BP algorithms. We begin by briefly discussing the unconditional stability of IL-prox, and equivalent G-IL algorithms. These algorithms use the NLMS update rule, so we next analyze the case where the LMS rule is used in an IL algorithm to update weights, as this case is more easily compared to G-BP algorithms which perform and LMS update over weights. In particular, we analyze how the FF output layer values h N change after weight updates are applied to a linear network using IL-GN with the LMS rule and using an analogous BP algorithm, which we call BP-GN. We show IL-GN, trained with the LMS rule, pushes output layer activitiesĥ N toward the target y (down it loss gradient) for any positive learning rate, while BP-GN only does this for a small finite range of learning rates. We focus on linear networks because constraining neuron activities with non-linearites may improve the stability of an algorithm by limiting the possible values a neuron may take or, e.g., by decreasing the magnitude of weight updates. We want to separate the effects of learning rules/algorithms on stability from those imposed by architectural constraints, so we focus on linear networks.
These theorems show that the way IL-GN computes and uses targets in its weight updates is a key mechanism that aids stability. In particular, IL-GN computes weight gradients for W n using pre-synaptic targetĥ n , while BP-GN uses pre-synaptic FF activities h n . This is the main difference between the two algorithms. The LMS update in the IL algorithm in linear networks is e n+1ĥ T n = (ĥ n+1 − W nĥ T n )ĥ T n , while the BP update is e n+1 h T n = (ĥ n+1 − W n h n )h T n . IL updates thus 'chain together' the local learning problems such that the local prediction target at one layer, is the input to the next local prediction problem. G-BP updates do not have this property. The input to one local prediction problem is the FF activity, which is computed independently from local prediction targets. These results show that the way IL chains together the local prediction problems plays an important role in its stability. In particular, it shows it plays an important role in ensuring that output layer FF activities h N change in the desired direction across a large range of learning rates. If we update weights like BP and do not chain the local learning problems together in this way, we lose this stability guarantee.
Ensuring h N changes in the desired direction for any learning does not guarantee the loss will be minimized, since, e.g., h N may overshoot the global target significantly. As we explain next, however, one can ensure minimization across any positive learning rate using the strategy of IL-prox, which uses the NLMS rule and uses the learning rate to control the output layer target rather than to scale the weight updates.

D.1 The Unconditional Stability of IL-prox
Implicit SGD is unconditionally stable in the sense that the proximal algorithm (equation 11), which is equivalent to performing an implicit SGD update, monotonically decreases the loss function L for 0 < α. This can be seen from the proximal update, equation 11, from which it trivially follows . Theorem B.2 shows that IL-prox and G-IL, under the γ settings specified in theorem 4.1, are equivalent to an implicit SGD update, and thus are unconditionally stable. Our simulations find that our implementation of IL-prox indeed displays this property of unconditional stability (tables 5 and 4). It is worth, however, briefly discussing the mechanics of how this unconditional stability comes about in the IL-prox algorithm (and equivalent G-IL algorithms). First, theorem 4.1, tells us that the learning rate α of IL-prox (and G-IL under certain γ settings), is only used to determine how muchĥ N is pushed toward global target y during the inference phase. More specifically the gradient of the proximal update w.r.t.ĥ N is ∂P rox where e N =ĥ N − p N and p N = W N −1 f (ĥ N −1 ). Let's assume a linear network and L = 1 2 y −ĥ N 2 . One way to compute the update ofĥ N would be to take the negative of the gradient, set equal to zero and solve forĥ N : Thus, each iteration of the inference phase we updateĥ N to a weighted average between p N and y. The learning rate is only used to determine the weighting between these two terms. As α → ∞, h N → y and as α → 0,ĥ N → p N . According to lemma B.1,ĥ N becomes the FF output layer value after the weights are updated (since the NLMS rule is used by IL-prox). Thus, when α → ∞, IL-prox finds a set of weights that best minimize the loss and as α → 0, weight are unchanged, which is the exact same property the proximal algorithm has (equation 11). Intuitively, we see that in any case the loss is either minimized or remains unchanged (when α = 0). IL-prox, and the proximal algorithm, thus gain their stability from the fact they do not use the learning rate to scale weight updates (as explicit SGD does) but rather uses it to determine the relative weighting of the loss term and regularization term in the proximal operator.

D.2 Gauss-Newton IL Updates are Minimum-Norm and Stable
Here we show that the IL-GN weight updates collectively minimizes global loss along its minimum norm path for any positive learning rate in linear networks. This theorem is proven true under the conditions whereĥ T n p n > 0, h T nĥn > 0, andĥ T nĥn >ĥ T n p n for all n. This means local predictions, p n , and initial activities, h n must be within 90 degrees of sub-targets. This is an easily satisfied condition as long as targets are not moved far from initial activities. Also, the inequalitŷ h T nĥn >ĥ T n p n is generally true since the magnitudes ofĥ n and p n will be similar andĥ n is always more similar to itself than p n . Thus, these conditions generally hold in practice.
Theorem D.1. Assumeĥ T n p n > 0,ĥ T n h n > 0, andĥ T nĥn >ĥ T n p n for all n. For a mini-batch of size 1, the IL-GN weight updates applied to a linear MLP at iteration b collectively push h N down its global loss gradient toward the target for any positive learning rate α: h Proof. We define the linear network after weight updates are applied using a recursive formulation. LetŴ * n =Ŵ * n+1 (W n − ∆W n ) for all n < N − 1, and letŴ * N −1 = W N −1 − ∆W N −1 . We can now express the output of this updated network given input h 0 asŴ * 0 h 0 . We assume a linear gradient update ∆W n = −αe n+1ĥ T n , where α is the learning rate and e n+1 =ĥ n+1 − p n+1 . The feedforward pass of the updated MLP at iteration b + 1 is described as follows: Note thatĥ 0 = h 0 and is the same across all training iterations, butĥ n = h n for n > 0.
Notice that the first termŴ * 1 h 1 can be expanded recursively using the same expansion ofŴ * 0 h 0 in equation 51: The leftmost terms h N + c N e N was computed as follows: Now we need to expandŴ * n for all n in equation 52.
We now apply the same derivation ofŴ * n e (b) n to allŴ in equation 52 to yield the following: Note that ( where j = c N + N −1 i=n c i g i and j > 0 is both c > 0 and g > 0. Hence, the weight updates applied at iteration b will collectively move the FF output layer values h N down its loss gradient h N − je It should be noted that moving h N down its loss gradient will not necessarily minimize the loss e N , given that j may be large and could cause h N to significantly overshootĥ N preventing convergence.

D.3 Gauss-Newton Back Propagation is Minimum-Norm but Only Conditionally Stable
In this section, we show that, unlike IL-GN, weight updates performed by an analogous BP based network, we call BP-GN, are not guaranteed to collectively move the global prediction down its loss gradient each training iteration. Meulemans et al. [28] showed that each individual weight update in a BP-GN network will push the global prediction h N down its loss gradient when its effects on h N are considered independently of other weight updates in the network. Here we show this property does not hold for any learning rate when the weight updates are considered collectively.
Definition D.1. Gauss-Newton backpropagation (BP-GN). Assume e n =ĥ n − h n . The weight update of BP-GN is the general BP weight update: ∆W n = −e n+1 h T n . Local targets are computed at each hidden layer are computed using the following equation: We can now apply the same analysis we did for IL-GN to BP-GN, which proves the following. Theorem D.2. Consider a linear MLP trained with BP-GN. For a mini-batch of size 1, the BP-GN weight updates applied at iteration b only push the global prediction down its global loss gradient h for a finite range of α.
Proof. Here we follow the same notation and analysis as the proof for theorem D.1. LetŴ * n = W * n+1 (W n − ∆W n ) for all n < N − 1, and letŴ * We assume a linear gradient update ∆W n = −αe n+1 h T (as in definition C.1), where α is the learning rate.
First, we can describe the FF pass of the updated MLP trained with GN-TP at iteration b + 1 as we did above for IL-GN (see equation 51): Let c n = αh T n−1 h n−1 such that αŴ * 1 e 1 h T 0 h 0 = c 1Ŵ * 1 e 1 . Because 1 > α > 0 and because h T n h n > 0 for all n, c n > 0. We can see this equation is the same as equation 51, except e 1 and ∆W 0 were defined according to GN-TP rather than IL-GN.
As above, the first termŴ * 1 h 1 can be expanded recursively: Equation 52 shows that for IL-GN, each of these recursively computed terms push the initial prediction down the global loss gradient. However, the same is not true here of GN-TP: W n e n can be rewritten in terms of e n+1 as follows: W n e n = W n W + n e n+1 . Notice that it will often be the case that h T nĥn < h T n h n , since generally h n = h n . Thus, α(h T nĥn − h T n h n ) will often be a negative scalar. This is distinct from the corresponding term in equation for IL-GN α(h T n h n −ĥ T n h n ), which will generally be a positive scalar. The difference is due to the IL-GN update being conditioned on the presynaptic sub-targetĥ n rather than presynaptic FF activity h n . Also notice there is no further dampening term γ(1 − γ) as there is with IL-GN.
Let k n = α(h T nĥn − h T n h n ) and d n = (1 + k n ). For reasons just noted, d n is not guaranteed to be positive for all 1 > α > 0. Continuing the derivation above: Because d n is not guaranteed to be positive, then ( is not guaranteed to be a positive scalar. Thus, the product of all d n may or may not be positive for a given 1 > α > 0. Let g n = ( Also, note that e N = − ∂L ∂h N . We can now rewrite equation 57 in terms of e N : We can see that if enough g n are negative it is possible that j is also negative, especially given that g n are scaled by c n , which will always be positive and may be larger than 1. Thus, GN-TP, unlike IL-GN, does not guarantee that for any positive value of α, the weight updates at iteration b will collectively move h N down it loss gradient at b + 1: GN-TP does not guarantee j > 0 for any α > 0 or even any 1 > α > 0. In order for j > 0, a large number of g n must be positive to offset any negative g n . In order for g n > 0, an even number of d n in g n = ( N −1 i=n d i ) must be negative. The simplest way to guarantee this is the case is to ensure d n > 0 for all n. In order for all d n to be positive in each g n we can see that alpha must be small since d n = 1 + k n = 1 + α(h T nĥn − h T n h n ), and α must be small enough such that k n > −1 for all n. This implies learning rate must be small: α <

E G-IL Avoids Interference Effects
The last section showed the IL-GN will push output layer values toward the global target for any positive learning rate, while BP-GN will only do so for a finite range of learning rates. Here we provide mathematical intuition for why this occurs. Specifically, we describe how BP weight updates necessarily interfere with each other's ability to minimize the global loss and this interference generally grows with learning rate, whereas IL weight updates do not necessarily interference and may actually improve each other's ability to minimize loss. In this sense, IL updates are more compatible with each other than BP updates. In line with this intuitive description, we show that under certain assumptions IL-GN updates (at hidden layers) are often most effective at pushing the output layer values toward the target when applied in conjunction with other updates than when applied alone, while the opposite is true of BP-GN. Figure 6 provides evidence these theoretical results hold true in practice for IL-SGD and BP-SGD networks when small learning rates are used.

E.1 A Mathematical Intuition for Interference Effects
Consider the BP weight update: ∆W n = α n e , where h n = W n−1 f (h n−1 ) and n > 0. Now, consider the effect of this weight update on the FF values at hidden layer n + 1 given the same We can see that the ability of the weight update to drive its output down the error gradient (where e , see section 3.2) depends directly on the similarity between pre-synaptic activities before and after the weight update. In particular, if the two are similar, i.e. f (h ) > 0, then the weight update will affect FF values in the desired direction. However, if the two are orthogonal, i.e. f (h there will be no effect, and if dissimilar, i.e. f (h ) < 0, they will affect FF values in an undesired direction (up rather down the loss gradient).
The challenge for the BP update is that if all weights are updated at iteration b then FF activities will change, which will generally make f (h ) less similar. Additionally, presynaptic activities will generally become less similar as α n is increased because changes to weights, and thus changes to FF activities, at previous layers will increase. Only when learning rates at previous because an α of zero means no change in weights or FF activities. We can describe this as a sort of interference effect: non-zero weight updates elsewhere in the network reduce the ability of ∆W n to drive FF activities at layer n + 1 in the desired direction. 2 Now, consider the G-IL weight update: where n > 0. The effect of this weight update on the FF values at hidden layer n + 1 given the same input, x (b+1) = x (b) , can be described as follows: Unlike with G-BP we can see that the ability of the weight update to drive its output down the error gradient e n+1 depends directly on the similarity between pre-synaptic optimized activities and the FF activities after the weight update. If the two are similar, i.e. f (ĥ ) > 0, then the weight update will affect FF values in the desired direction. If the two are orthogonal there will be no effect, and if they are dissimilar they will affect FF values in an undesired direction.
What's interesting about the G-IL update is that with a properly tuned α n , f (ĥ n ) (b) and f (h n ) (b+1) will get more similar when weights prior to layer n are updated. As lemma B.1 shows, under a specific (typically) non-zero value of α n it is actually the case that f (ĥ ). This value of α n is the true solution to G-IL update argmin W F that the LMS learning rule above attempts to approximate. Intuitively, G-IL largely avoids interference because f (ĥ

E.2 G-IL Weight Updates are Most Effective when Applied Together
In previous sections we compared the direction of the effects of weight updates on a linear MLP output layer in an IL and BP network. We now assess the magnitude of the affects of IL and BP weight updates on the output layer of a linear MLP with a single hidden layer. In sum, we show that under certain assumptions, updates in IL-GN should have greater effects on the output layer when applied together than when applied alone, while the opposite is true for BP-GN. This shows that interference effects not only can alter the direction the output layer changes with weight updates but also the magnitude of the change in output layer activities. 2 Interference might be reduced by using non-linearities that keep activities positive, e.g., ReLU, which However, for non-zero learning rates it would still be the case that updates to weights at layers before n would reduce similarity, i.e., generally f (h be the output layer values after W 0 is updated alone, h (b+1/2), 1 2 be the output layer values after W 1 is updated alone, and h (b+1) 2 be the output layer values after both weight matrices are updated. Let c n and g n be defined as above in the proof of theorem D.1. We assume that g n > 1 which will often (though not always) be true in practice for reasons noted above (see theorem D.1). We also assume c 1 = c 2 as it significantly simplifies the calculations. We test how well these theoretical results hold in practice and find they hold approximately for small learning rates (see figures 6 for details). Proposition E.1. Consider a linear MLP with a single hidden layer trained with IL-GN, mini-batch size 1, at iteration b. Assume g 1 > 1, c 1 = c 2 . Under these assumptions, and those in theorem D.1, each weight update, ∆W 0 and ∆W 1 has a larger effect on the output layer when applied in conjunction with each other than when applied alone: h Proof. The effect of ∆W 0 on the output layer h N = h 2 when applied alone is where the second line is computed using C.3, which says e n = γW + n e n+1 . Here c 1 = αh and is clearly a positive scalar.
The effect of ∆W 1 on the output layer, h 2 , when applied alone is where c 2 = αĥ 1 and is a scalar. Finally, when both updates are applied we get the same expression as that in equation 56 targets are computed and used in the weight update. For example, while BP-prox uses presynaptic FF activity h n in its update of W n , IL-prox uses presynaptic target valueĥ n in its update of W n (compare equations 74 and 76). Clearly, these differences have a significant effect on stability.   Table 6: Test accuracies after 1 epoch of training with mini-batch size 1 (60,000 training iterations). Fully connected networks were used with dimensions 784-2x500-10.
The training runs for are plotted in figure F.2 below. Differences in performance between IL and BP algorithms are smaller with these data-sets than they are with CIFAR-10. IL-prox improved accuracy slightly faster in the first few hundred or thousand iterations on both MNIST and Fashion-MNIST, but BP algorithms soon caught up.

F.3 Convolutional Networks
We also trained a small convolutional networks on CIFAR-10. Large mini-batches (size 64) were used and the Adam based algorithms were used for training (see 2), since IL algorithms worked best with Adam optimizers when large mini-batches were used. IL-prox Adam achieved essentially the same accuracy as BP-Adam, though IL-SGD did not achieve comparable accuracy.  Table 7: Here we train on CIFAR-10 with mini-batches of size 64 for 77,000 iteration (≈ 100 epochs). Runs were averaged over five seeds and the best score is shown. There were three convolutional layers followed by one fully connected layer. We used a similar architecture as [3]. Convolutions were (5x5, 64, 2), (5x5, 128, 2), (3x3, 256, 2). Fully connected layer was 1024x10. ReLU activations were used at hidden layers, and softmax at the output layer.  N using cosine similarity. For small linear networks trained on a regression task with, IL-SGD does indeed produce changes in the output (see F.4 left) very close in angle to the minimum norm path. This supports the claim that IL-GN is a useful approximation of G-IL. MLPs trained with BP-SGD produce changes in the output layer that are significantly less close to the minimum norm path. This suggests IL-SGD takes a shorter, more direct path toward local minimum than that of steepest descent. Similar results were found when ReLU activations were applied at hidden layers (see F.4 middle). Finally, we directly test theorem D.1 and theorem D.2 by exactly implementing IL-GN and BP-GN in a toy linear network. We measure the cosine similarity between the change in output FF values and the minimum norm change, as with the previous tests, and we do this over a range of learning rates. As predicted by theorem D.1 and D.2, IL-GN and BP-GN affect output layer values along minimum norm path toward the target. However, as predicted by theorem D.1 and D.2, we find BP-GN is significantly less stable than IL-GN (see F.4 right).
Propositions E.1 and E.2 show that under certain assumptions IL-GN weight updates should have a larger effect on h N when applied together than when applied alone, and the opposite is generally true of BP-GN. We test how well this result holds in practice in IL-SGD and BP-SGD networks trained on a regression tasks. We compute the effect each individual ∆W  When h b+1 N , we know ∆W n had a larger affect on the output layer when applied in conjunction with other weights than when applied alone. When ∆W n meets this condition we count it as compatible with the other weight updates. We then compute a compatibility score which stores a count of all compatible weight updates and uses it to compute the proportion of updates that were compatible.
We find that less than 50% of weight updates in IL-SGD and BP-SGD networks are compatible though for smaller learning rates a greater proportion of IL updates are compatible than BP updates (see figure 6). This supports the claim that IL updates interfere with one another less than do BP updates, both in their ability to minimize loss and in their ability to change the FF output layer values. Figure 6: (Left) In linear networks with small learning rates (α n < .1) IL weight updates have significantly larger compatibility scores (according to 2-sample T-test, 95% confidence intervals shown but are too small to see) than BP-SGD updates. (Right) The same trend exists in non-linear networks with ReLU at hidden layers for α n < .4.

F.5 Weight Norm Analysis
Here we report the mean and standard deviation of the parameter updates ∆θ from the training runs reported in 6. We can see that generally IL updates are on average smaller than BP-SGD updates and have smaller standard deviations. IL algorithms tended to improve accuracy faster than BP-SGD.
This suggests that IL did not improve accuracy faster because the overall changes in parameters were larger. One possibility is that IL improved accuracy faster because it updated parameters along a more direct path toward local minima. This suggestion is consistent with the findings above, that IL updates push FF output layer activities along a more direct path toward global target activities y.

G Methods
In this section, we provide a more detailed description of the algorithms used in our simulation, as well as a more detailed description of the experiments, so that they may be replicated in the future. We begin by describing the algorithms in more detail then discuss experiment details. 3 An Adam optimizer is assigned to each weight matrix. We used the standard/default hyper-parameters settings in the pytorch Adam optimizer: β 1 = .9, β 2 = .999, = 1 * 10 −8 .

IL-prox
This algorithm is our novel variant of IL that is made to more closely match the algorithm described in definition 4.1. Theorem 4.1 shows that minimizing the proximal loss w.r.t. hidden layer activities and output layer activities is equivalent to minimizing F under certain settings of the γ terms 4.1. The most important difference between these γ settings and the ones of IL-SGD is the way the learning rate is used. In particular, IL-SGD scales weight updates using the learning rate, whereas IL-prox scales weight updates by the norm of presynaptic activities using the NLMS rule (equation 15) instead of learning rate. Then, following theorem 4.1, IL-prox uses the learning rate α only to determine the degree to which h N is pushed toward global target y versus output layer prediction p N . In particular, following theorem 4.1, γ N = α N −1 α = 1 ĥ N −1 2 α , then the gradient update at the output layer is , and σ is the softmax in classification tasks and the sigmoid in the autoencoder task. In practice, we use a cross-entropy loss to update weights but find using gradients of the cross entropy loss to update output layer activities leads to some instability. Instead, we use the gradient of the MSE to update output layer activities. In particular, and we solve forĥ N to getĥ Thus, at each iteration of the inference phase, the output layer activityĥ N is updated to a weighted average between global target y and prediction p N . We see that as α → 0 thenĥ N → p N and as α → ∞ thenĥ N → y. Interestingly, we see the inverse relation betweenĥ N and the normalized learning rate. The normalized learning rate for W N −1 is ĥ N −1 −2 . We see that as ĥ N −1 −2 → 0 thenĥ N → y. This means that for very small normalized learning rates the output layer activity will be pushed toward the global target y, which will lead to larger output layer errors e N . As ĥ N −1 −2 → ∞ thenĥ N → p N . This means that for very large normalized learning rates the output layer activity will be pushed increasingly fast toward the output layer prediction p N , which will lead to smaller output layer errors. This illustrates how these activity updates manage weight update magnitude, as the magnitude of weight updates under the NLMS rule depend on the magnitude of post synaptic errors e N and the magnitude of the normalized learning rate. If the learning rate is small, e N can be larger and the weight update will still be small. If the learning rate is large, the activities are updated to keep e N small and thus the weight update small.
We find including the ĥ N −1 2 term in the output layer activity update is important for the stability and performance of IL-prox, but the weighting terms of theorem 4.1 are not particularly important at hidden layers for performance or stability of IL-prox. To reduce computation, we just use preset hyper-parameter γ values at hidden layers. IL-prox uses the NLMS rule (equation 5) to update weights. In practice, we find it sometimes improves performance to add a small constant to the NLMS rule such that it becomes The provides a smooth upper bound on the magnitude of the normalized learning rate and prevents division by zero errors. (When using biases, however, division by zero is generally not an issue, since one must treat the pre-synaptic vectorĥ n as having a 1 concatenated to the end of it which for most activation functions guarantees f (ĥ n ) −2 > 0.) In the case where mini-batches are used we first average the gradient across the batch, then multiply this averaged gradient by the average of the normalized learning rates across the mini-batch.
Algorithm 4: IL-prox begin // Feedforward Pasŝ h 0 ← x (b) for n = 0 to N − 1 do p n+1 ,ĥ n+1 ← W n f (h n ) end h N ← Eqn. 73 // Inference phase (Compute Local Targets) for i = 0 to 25 do for n = 1 to N dô h n ←ĥ n − β ∂F ∂ĥn p n+1 ← W n f (ĥ n ) end h N ← Eqn. 73 end // Update Weight Matrices Eqn. 74 end IL-prox Fast is the same as algorithm 4 except it only performs 12 updates over activities. IL-prox Adam uses IL-prox to compute (approximate implicit) gradients (i.e. uses the same inference phase as IL-prox then equation 74 to compute weight gradients), then inputs these gradients to an Adam optimizers, which updates the weights.
BP-prox was the algorithm we used as a control to compare against IL-prox in the stability tests. BP-prox softly clamps the output layer similarly to IL-prox and uses the NLMS rule similarly to IL-prox. The output layer activities are computed according tô Notice, unlike IL-prox, the FF activity h N is used rather than p N since G-BP does not compute local predictions. However, IL-prox initializes p N to h N so these values are the same upon initialization.
Here h N = σ(W N −1 f (h N −1 )) where σ is the softmax. The weights are updated as follows: where is a small scalar used to prevent division by zero and place upper bound on the normalized learning rate. (As we note below, in the stability tests the same small is used for both IL-prox and BP-prox for a fair comparison). Notice this is the same as the NLMS rule used by IL-prox (equation 74), except, in keeping with the G-BP framework, the pre-synaptic FF activities h n activities are used in the learning rate and gradient com-putation rather than target activitiesĥ n . A summary of the algorithm is presented below. Pytorch data-loaders were used to download and perform transformations on MNIST, Fashion-MNIST, and CIFAR-10 datasets. All data was transformed to pytorch tensors, but no other alterations were applied to the data (e.g., there was no normalization). All pixel values were in the range between 0 and 1. Pytorch auto-gradient library was used to compute gradients for BP models, and the local gradients for IL-SGD models. Gradients used to update weights and activities in IL-prox models were computed manually. Code is publicly available at ???.
Supervised Tasks For all classification tasks, ReLU activations were used at hidden layers and softmax at the output layer. For MNIST tasks, fully connected networks size 784-2x256-10 were used. For our main results with CIFAR-10 classification (table 2), we used fully connected networks size 3072-3x1024-10. For the convolutional classification we used a similar architecture as [3], who compared a wide range of local,biologically plausible learning algorithms. Convolutions were (5x5, 64, 2), (5x5, 128, 2), (3x3, 256, 2). Fully connected layer was 1024x10. (Note, these are the same dimensions as [3], though they used locally connected layers rather than convolutions.) A grid search was used to find the learning rates for CIFAR-10 and MNIST classification tasks shown in tables 2, 6, and 7. We first hand tuned the learning rates to find a range over which the algorithms performed the best, then chose 5-10 learning rates in the range and trained 2-5 seeds at each learning learning rate. We averaged the test accuracy achieved at the end of training over the seeds, then chose the learning rate with the highest average test accuracy at the end of training. We then reran the training run with the best learning rate across 5 seeds to get the data shown in tables 2, 6, and 7. Each seed was tested every 50 training iterations. Runs were averaged over the 5 seeds. The best test accuracy achieved from the averaged training run, is shown in the tables.
Learning rates used for each algorithm are shown below. The other main hyper-parameters are the γ terms in the energy equation (equation 3), the gradient step size used to update activities in the IL models, and the values in the IL-prox update (equation 74). The γ terms and the step sizes over activity updates in IL models, each scale the gradients of F. To simplify, we found it effective to just set γ decay = 0, then use the same two hyper-parameter terms at each hidden layer to differentially weight the two error gradient terms: γ bot and γ top . First, γ bot is multiplied by the 'bottom up' error gradient in equation 23, which comes to γ bot f (ĥ n )W T n e n+1 . Then we γ top is multiplied by the 'top-down' error gradient in equation 23, which comes to γ top e n . These same two hyper-parameters are used at each hidden layer of IL models and can be seen as incorporating both the actual gamma and step size hyper-parameters into a single hyper-parameter. For IL-SGD and IL-Adam, γ bot = .02 and γ top = .015. For IL-prox and IL-prox Adam, γ bot = .015 and γ top = .015. For IL-prox Fast, γ bot = .015 and γ top = 0. Finally, = .25 for IL-prox and IL-prox fast models.
Self-Supervised Tasks The same grid search procedure was used for the self-supervised task as for the supervised tasks to find learning rates for each algorithm. ReLUs were used at hidden layers and  Table 8: Learning rates for algorithms trained on supervised learning tasks with fully connected networks, with various data sets and mini-batch sizes (mb).
sigmoid at the output layer. The autoencoder used fully connected layers dimensions 3072-1024-500-100-500-1025-3072. Binary cross entropy (BCE) was used at the output layer. The BCE data shown in figure 2 was summed over pixels and averaged over the test set. The best test score from each of these averaged training runs is shown in the table. Figure 2 shows a run with mini-batch size 1, which were trained with 50,000 data point (1 epoch). Learning rates are shown in the table below. The γ settings are the same as those from the previous section. For the NLMS rule in IL-prox we used an = 20 and generally found it useful to have a larger (and thus a smaller upper bound on the normalized learning rate) than in the supervised task.
Importantly, the accuracy shown in table 4 and 5 are the test accuracy, averaged over the 5 seeds, at the end of training. We use the test accuracy at the end of training (rather than the best test accuracy achieved during training). The purpose of this was to test which algorithms not only achieved above chance accuracy, but also did so while converging in a stable manner. Also important, the values used in the IL-prox and BP-prox updates were set to zero to ensure any stability these algorithms showed in this test was not due simply to the down scaling the weight updates.
Output Layer Analysis For the output layer analyses presented in figures F.4, 6 we trained small toy models on a regression task. Models were FF networks with layer sizes, 10-5-5-5-5. Input data was generated by sampling vectors from a standard Gaussian. A teacher network of the same size with Tanh activations at hidden layers, was used to generate associated output target vectors associated with each input vector.
Linear toy networks and toy networks with ReLUs at hidden layers were trained with BP-SGD and IL-SGD across learning rates (.001, .005, .01, .05, .1, .25, .5, .75, 1, 1.5, 2, 3, 4, 5). Each algorithm trained 50 seeds at each learning rate for 200 training iterations. Each training iteration the compatibility score (equation 70) was computed. Additionally, each training iteration b, we computed the minimum norm path between the output layer activities and the output target y (b) as y (b) − h N after the update each training iteration. We then computed and stored the cosine similarity between the two. This data was then averaged over each iteration and seed for each learning rate. The result was the average compatibility score and minimum norm change similarity under each learning rate.


extfieldsofstudy:
Physics
Computer Science

text:
QuASK -- Quantum Advantage Seeker with Kernels

QuASK is a quantum machine learning software written in Python that supports researchers in designing, experimenting, and assessing different quantum and classical kernels performance. This software is package agnostic and can be integrated with all major quantum software packages (e.g. IBM Qiskit, Xanadu's Pennylane, Amazon Braket). QuASK guides the user through a simple preprocessing of input data, definition and calculation of quantum and classical kernels, either custom or pre-defined ones. From this evaluation the package provides an assessment about potential quantum advantage and prediction bounds on generalization error. Moreover, it allows for the generation of parametric quantum kernels that can be trained using gradient-descent-based optimization, grid search, or genetic algorithms. Projected quantum kernels, an effective solution to mitigate the curse of dimensionality induced by the exponential scaling dimension of large Hilbert spaces, are also calculated. QuASK can furthermore generate the observable values of a quantum model and use them to study the prediction capabilities of the quantum and classical kernels.

Introduction
Is there any clear benefit in using quantum kernels [1] to enhance machine learning applications? The definitive answer has still to be found [2]. The nofree lunch theorem guarantees that there is no best model in machine learning which outperforms all others, and efficiency is based on the type of training data distribution. Thus, there is no hope for a "universal" quantum kernel. However, quantum kernels have been proven to reach superior performance in learning from quantum experiments [3,4] and in solving artificially generated, classically intractable classification tasks [5]. Considering real-world data, they have been used to classify cosmological data with performance comparable to the classical Gaussian kernel [6]. Moreover, these can be effectively implemented on NISQ devices [6,7]. Quantum kernels must be carefully tailored to avoid problems related to the curse of dimensionality: the Hilbert space associated with a kernel scales exponentially in the number of qubits, leading to exponentially small inner products. This issue could be solved by projecting the data onto low-dimensional subspaces using projected quantum kernels [3].
The process of designing, experimenting, and evaluating different quantum kernels can be nontrivial, resulting in computationally intensive and errorprone calculations. The dataset must be appropriately preprocessed in order to reduce the number of features to accommodate the available resources of NISQ hardware and simulators. Furthermore, the features must be normalized and rebalanced according to the class distribution. The choice of the unitary transformation implementing the quantum feature map is arbitrary, as is the observable choice in projected quantum kernels. In addition, this unitary transformation may be parametric and optimization might be necessary according to some criteria [8]. Finally, the quantum kernel must be evaluated according to some metrics including the prediction accuracy, model complexity, and geometric difference [3].
The whole process described above can be automated with the proposed QuASK (Quantum Advantage Seeker with Kernel) library. This package, written in Python, implements a (quantum) machine learning pipeline that takes care of the most trivial but time-consuming part of working with quantum kernels. QuASK is available either as an interactive CLI (Command-Line Interface) or as a library whose components can be used to build new software. QuASK allows the complete reproducibility of the experiments, saving the information processed at any step of its pipeline. QuASK allows, also, the generation of data from quantum processes in order to study the performance of classical and quantum kernels using quantum data. Finally, in terms of results visualization, QuASK can generate a simple accuracy plot of all the trained kernel machines.

QuASK library
QuASK provides a library of elements that can be integrated with other projects. The library is composed of the following software modules.
quask.datasets and quask.metrics The quality of a quantum kernel can be empirically tested through the performance of a kernel machine with respect to a certain dataset, by evaluating some metrics. The quask.datasets module facilitates the researcher in choosing a suitable dataset, providing some of the most popular datasets from OpenML platform 1 and custom datasets generated also from quantum experiments (the latter allows us to reproduce results in [3]). The module quask.metrics contains the metrics to compare and evaluate the kernels, including the kernel polarity (i.e. Frobenius inner product between two Gram matrices), the targetkernel alignment [9], the training and testing accuracy of the Support Vector Machine with the precomputed kernel, the geometric difference [3] (which can be used to find a potential quantum advantage), and the model complexity [3].
quask.template pennylane and quask.template qiskit The quantum kernel requires the definition of a feature map φ which is implemented using a parameterized unitary transformation U φ (x). The Quantum Kernel Estimation algorithm [10] calculating the function k(x, x ) = φ(x)|φ(x ) is implemented through either the overlap test (Figure 1a) or the C-SWAP test (Figure 1b). The projected quantum kernel calculates classically the inner product between two feature vectors φ(x 1 ), φ(x 2 ), each one being the output of the quantum feature map φ(x) = 0|U † (x)HU (x)|0 (Figure 1c). The feature map crosses the quantum space first through U(x) and projects the data back in a classical representation when measuring with the Hermitian operator H. QuASK contains both some notable unitary transformations U from the literature and the code to use such unitary transformations as a kernel function through one of the three methods described above. The user can define its own unitary transformation and immediately get the kernel function. QuASK is agnostic with respect to the software framework used to define, simulate and execute the quantum circuits: we have implemented some unitary transformations in PennyLane (quask.template pennylane module) and some in Qiskit (quask.template qiskit). This allows also the use of the different functionalities offered by the different frameworks. For example, noiseless simulation with PennyLane can be speeded up using JAX 2 .

quask.kernels
The module quask.kernels keeps track of the quantum kernel available within the framework. The module quask.kernels contains the global object the kernel register, which is a list of kernels. This list is already filled with  The processing phase can analyze the system to find the best technique to analyze the system.
all the quantum kernels derived from the unitary transformations in modules quask.template xyz, while the user can add its own custom kernel too. The kernels do not need to be quantum. Such module can use trainable kernels (defined in quask.template pennylane) which are parametric quantum circuits that depend on both the input and some trainable parameters. These parameters are supposed to be trained with gradient-descent-based optimization algorithms. The user can take advantage of efficient optax optimization library. Nonetheless, the user can manually tune the parameters or even run a grid search. Furthermore, such an approach is extendable to optimization techniques that change the topology of the circuit itself such as genetic algorithms.

QuASK script
QuASK can be run as a standalone, interactive executable. The software is run with python3 -m quask <command>. The script implements the sequence of operation described in Figure 1d. The available commands are:

get-dataset command
The script is used to retrieve one of the datasets available in quask.datasets module. The output of the process is a pair of NumPy binary files representing the feature data and the corresponding labels.

preprocess-dataset command
The script can modify the dataset in several ways. Firstly, the user can vertically slice the dataset, keeping only a certain range of labels. Secondly, the user can apply dimensionality reduction techniques. These are especially important in the NISQ setting due to the lack of resources. The techniques available are PCA for numerical data and FAMD for mixed numerical and categorical data 3 . Thirdly, it is possible to fix the possible imbalanceness of the classes using random undersampling or random oversampling. When loading, the script already shows some statistics about the dataset, both for classification and regression tasks, which can guide the user through the preprocessing. The output of the process are the four files X train, y train, X test, y test which can now be fed to some kernel machine. > python3.9 -m quask preprocess-dataset Choose a random seed: 2323 Choose a random seed for splitting training and testing set: 3232 Where is the output folder: . Specify the path of classical feature data X (.npy array): iris_X. The user can use trainable quantum kernels. In this case, the user specifies the embedding and the variational form (and its number of repetitions), the optimizer (ADAM stochastic gradient descent or grid search), and the reward function to maximize.

plot-metric command
The script calculates the accuracy of a kernel machine using the training and testing Gram matrices given as input. The output is a plot comparing the different kernels. For each kernel matrix, the user specifies the label that appears at the x-axis of the plot. If multiple instances are specified with the same label these are interpreted as i.i.d. random experiments and will contribute to the error bars. The software also calculates, together with the accuracy, any metric present in quask.metrics.

Conclusions
We have introduced QuASK, a tool supporting researchers in creating powerful quantum kernels. The software takes care of the most time-consuming and errorprone aspects of the experimentation. It exploits theoretical metrics in QML, providing users with an environment to easily assess case for potential quantum advantage. This package offers the exciting perspective of testing these metrics on real-wolrd datasets. The QuASK project will be extended in future versions with a wider range of datasets and feature maps, both classical and quantum.