One Forward is Enough for Neural Network Training via Likelihood Ratio
  Method
While backpropagation (BP) is the mainstream approach for gradient
computation in neural network training, its heavy reliance on the chain rule of
differentiation constrains the designing flexibility of network architecture
and training pipelines. We avoid the recursive computation in BP and develop a
unified likelihood ratio (ULR) method for gradient estimation with just one
forward propagation. Not only can ULR be extended to train a wide variety of
neural network architectures, but the computation flow in BP can also be
rearranged by ULR for better device adaptation. Moreover, we propose several
variance reduction techniques to further accelerate the training process. Our
experiments offer numerical results across diverse aspects, including various
neural network training scenarios, computation flow rearrangement, and
fine-tuning of pre-trained models. All findings demonstrate that ULR
effectively enhances the flexibility of neural network training by permitting
localized module training without compromising the global objective and
significantly boosts the network robustness.
Jinyang Jiang, Zeliang Zhang, Chenliang Xu, Zhaofei Yu, Yijie Peng
# Training Neural Networks without Backpropagation:

###### Abstract

Backpropagation (BP) is the most important gradient estimation method for training neural networks in deep learning. However, the literature shows that neural networks trained by BP are vulnerable to adversarial attacks. We develop the likelihood ratio (LR) method, a new gradient estimation method, for training a broad range of neural network architectures, including convolutional neural networks, recurrent neural networks, graph neural networks, and spiking neural networks, without recursive gradient computation. We propose three methods to efficiently reduce the variance of the gradient estimation in the neural network training process. Our experiments yield numerical results for training different neural networks on several datasets. All results demonstrate that the LR method is effective for training various neural networks and significantly improves the robustness of the neural networks under adversarial attacks relative to the BP method.

## 1 Introduction

Neural networks have increasingly more successful applications, such as face recognition [1; 2; 3; 4], voice verification [5; 6; 7; 8] and automation vehicles [9; 10; 11]. Different model architectures are suitable for processing different tasks. For example, multi-layer perceptron (MLP) [12] has a good performance on structural data, convolutional neural network (CNN) [13] is typically adopted for image tasks, recurrent neural network (RNN) [14] is usually employed to process time series data, and graph neural network (GNN) [15] is often used to extract potential relationships between entities. Brain-like neural network structure includes spiking neural networks (SNNs) [16] where neurons communicate via discrete spike sequences, and they are ideal models for event-driven information processing.

Backpropagation (BP) [12] is the most important gradient estimation method used for training neural networks [17; 18]. It relies on the chain-rule of differentiation to backpropagate the residual error from the loss function to the previous layers for gradient computation. Despite its widespread use, recent research has highlighted the vulnerability of neural networks trained by BP, particularly under adversarial attacks [19; 20]. Adversarial attacks refer to the subtle manipulation of input data that is imperceptible to humans but can easily deceive neural networks [21]. Meanwhile, BP requires the activation and loss function of the neural networks to be continuous, which limits the freedom of network architecture for modeling complex environments [22].

In contrast to BP, some existing methods entirely abandon the backward process and only require the forward pass for gradient estimation in optimizing the neural networks [23; 24]. One family
of training methods without BP performs optimization following the gradient estimated from the correlation between the parameter perturbation and network outcomes. Unfortunately, they usually suffer from curse-of-dimensionality and require significant manual tuning to identify appropriate hyper-parameters [25; 26], so their applicability to deep neural network training is limited.

The likelihood ratio (LR) method has emerged as a promising alternative to BP for training neural networks [27]. It provides an unbiased gradient estimation without relying on the recursive gradient calculation. Previous study has developed an LR method to train MLPs on the MNIST dataset. However, there is still significant work to be done in extending the LR method to more widely used neural network architectures.

In our paper, we develop a new LR method to train a broader range of neural network architectures. We further propose several methods to reduce the variance of gradient estimation, thereby enhancing the stability of the neural network training process. We conduct extensive experiments on various neural network architectures, i.e., the CNN, RNN, GNN, and SNN, for corresponding classical tasks to verify the effectiveness of the proposed method, including classification on the Cifar-10 [28], Ag-News [29], Cora [30], MNIST [31], and Fashion-MNIST [32] datasets. Our experimental results demonstrate the efficacy of our proposed method in achieving comparable or even better performance with BP while also improving the robustness of the trained models.

## 2 Related Work

### Optimizations of neural networks

Optimization methods for neural network training can be roughly categorized into two types, with and without BP. BP-based optimization methods have been studied for a long time and have been effective in training various neural networks on different tasks. Among these methods, the vanilla BP method is utilized for computing the gradients of the neural network parameters with respect to the global loss function. BP relies on the perfect knowledge of computation details between the parameters and loss evaluation, which limits the flexibility of developing neural network architectures and training. Moreover, the neural networks trained by BP are vulnerable to adversarial attacks.

Some studies have explored alternative methods without BP for training neural networks. [23] proposes a simultaneous perturbation stochastic approximation (SPSA) method, which approximates the gradient by evaluating the objective function twice with perturbations of opposite signs on the parameters. SPSA is typically sensitive to the choice of hyperparameters. Another approach is the evolution strategy (ES) [24], which injects a population of random noises into the model and optimizes the parameters along the optimal direction of noises. Both SPSA and ES suffer from heavy performance drops as the number of parameters increases, which hinders the application to sophisticated models in deep learning.

[27] proposes a novel LR method for training neural networks, which is based on an unbiased gradient estimation with respect to the true gradient computed by the BP. However, previous research has primarily focused on the training of MLPs. In our work, we extend the LR to a wider range of neural network architectures.

### Robustness of neural networks

Neural networks are vulnerable to various types of input corruption, including natural noise, out-of-distribution inputs, and adversarial attacks. Natural noise can be caused by various factors such as sensor noise, compression, or transmission errors. The benchmark proposed by [33] evaluates the robustness of neural networks against various types of natural noises, including Gaussian noise, uniform noise, Poisson noise, etc. Out-of-distribution (OOD) robustness is another challenge for real-world applications of neural networks due to the distribution shift between the training and test data. Adversarial attacks are created by adding imperceptible perturbations to benign samples, specifically designed to fool deep neural networks. Although adversarial attacks were first identified in the image domain, subsequent research has revealed their existence across various modalities such as texts, structural graphs, etc. For images, the most widely used methods are the fast gradient sign method (FGSM), iterative fast gradient sign method (I-FGSM), and momentum-based iterative fast gradient sign method (MI-FGSM), which employs the sign of gradient to craft the adversarial perturbation efficiently. For texts, [34] proposes the genetic algorithm (GA) to craft the adversarial text examples
with words substituted. [35] proposes the probability-weighted word saliency (PWWS) algorithm, which generates adversarial text examples with semantic consistency by synonym substitution. For graphs, [36] proposes to fool the GNNs by randomly injecting fake edges into the graph (R.A.). [37] proposes to craft adversarial graph examples by manipulating the graph to disconnect internally and connect externally (DICE).

In this work, we evaluate the performance of different optimization methods on datasets corrupted by different types of noises as a benchmark. Our proposed LR method shows improved robustness while achieving comparable results to BP on clean data.

## 3 Likelihood Ratio Method for Training Neural Networks

### Gradient estimation using push-out LR

To avoid the chain rule of differentiation, [27] introduces additional noises to the forward feeding in each neuron and employs the push-out LR method to estimate the gradient of the MLP. Given an \(L\)-layer MLP, the input to the \(l\)-th layer is denoted as \(x^{l}\in\mathbb{R}^{d_{l}}\), and the output is given by

\[x^{l+1}=\varphi(v^{l}),\quad v^{l}=\theta^{l}x^{l}+b^{l}+z^{l}, \quad z^{l}=(\Sigma^{l})^{\frac{1}{2}}\varepsilon^{l},\quad\varepsilon^{l} \sim\mathcal{N}(0,I),\]

where \(\theta^{l}\in\mathbb{R}^{d_{l+1}\times d_{l}}\) and \(b^{l}\in\mathbb{R}^{d_{l+1}}\) are the weight and bias terms, \(\varphi\) is the activation function, and \(z^{l}\in\mathbb{R}^{d_{l+1}}\) is a zero-mean Gaussian noise with covariance matrix \(\Sigma^{l}\). Training an MLP is to solve the optimization problem: \(\min_{\theta,b,\Sigma}\frac{1}{|\mathcal{X}|}\sum_{x\in\mathcal{X}}\mathbb{E }\left[L(x^{L})\right]\), where \(L\) is the loss function, and a basic approach is to perform gradient descent following the gradient estimation of \(\mathbb{E}\left[L(x^{L})\right]\).

By noticing that conditional on \(x^{l}\), \(v^{l}\) follows \(\mathcal{N}\left(\tilde{v}^{l},\Sigma^{l}\right)\) with the density \(f^{l}(\xi)\), where \(\tilde{v}^{l}=\theta^{l}x^{l}+b^{l}\), we can push the weight \(\theta^{l}\) out of the performance function of \(L(x^{L})\) and into the conditional density as below:

\[\mathbb{E}\left[L(x^{L})\right]=\mathbb{E}\left[\int_{\mathbb{R}^ {d_{l}}}\mathbb{E}\left[L(x^{L})|\xi,x^{l}\right]f^{l}(\xi)d\xi\right].\]

Under mild integrability conditions, the expectation and gradient operators can be exchanged. Thus, the gradient estimation only requires the density function \(f^{l}\) to be differentiable. By appling the LR technique, it is straightforward to show that

\[\nabla_{\theta^{l}}\mathbb{E}\left[L(x^{L})\right]=\mathbb{E} \left[\int_{\mathbb{R}^{d_{l}}}\mathbb{E}\left[L(x^{L})|\xi,x^{l}\right]\nabla _{\theta^{l}}\ln f^{l}(\xi)f^{l}(\xi)d\xi\right].\]

With calculation, we can obtain

\[\nabla_{\theta^{l}}\mathbb{E}\left[L(x^{L})\right]=\mathbb{E} \left[L(x^{L})(\Sigma^{l})^{-\frac{1}{2}}\varepsilon^{l}(x^{l})^{\top}\right], \tag{1}\]

which implies that the gradient estimation in the whole model can be neuron-wise paralleled with the final evaluation of \(L(x^{L})\) and the local information \(x^{l}\) and \(\varepsilon^{l}\). The gradient estimation for \(b^{l}\) is a corollary of Eq. (1), and the LR method also allows us to optimize the noise magnitude as below:

\[\nabla_{\Sigma^{l}}\mathbb{E}\left[L(x^{L})\right]=\frac{1}{2} \mathbb{E}\left[L(x^{L})\left((\Sigma^{l})^{-\frac{1}{2}}\varepsilon^{l} \left(\varepsilon^{l}\right)^{\top}(\Sigma^{l})^{-\frac{1}{2}}-(\Sigma^{l})^{ -1}\right)\right],\]

which could benefit the capacity of representing the environment and robustness of the model [38].

### Extension to more architectures

**Convolutional neural networks**: Consider the convolutional module in a CNN. Denote the input as \(x=(x^{i}_{j,k})\in\mathbb{R}^{c_{\text{in}}\times h_{\text{in}}\times w_{\text {in}}}\), where the superscript represents the channel index. The output \(v=(v^{o}_{j,k})\in\mathbb{R}^{c_{\text{out}}\times h_{\text{out}}\times w_{ \text{in}}}\) is given by

\[v^{o}_{j,k}=(x*\theta)^{o}_{j,k}+b^{o}+z^{o}_{j,k}=\sum_{i=0}^{ c_{\text{in}}-1}\sum_{s=0}^{h_{\text{g}}-1}\sum_{t=0}^{w_{\text{g}}-1}x^{i}_{j+s,k+ }\theta^{o,i}_{s,t}+b^{o}+z^{o}_{j,k}, \tag{2}\]
where \(\theta=(\theta^{o,i}_{j,k})\in\mathbb{R}^{c_{\text{out}}\times c_{\text{in}}\times h _{\theta}\times w_{\theta}}\) and \(b=(b^{o})\in\mathbb{R}^{c_{\text{out}}}\) are the weight and bias terms, \(z^{o}_{j,k}=\sigma^{o}\varepsilon^{o}_{j,k}\) is the noise injected to perform the LR method, \(\sigma^{o}\) denotes the noise magnitude, and \(\varepsilon^{o}_{j,k}\) is a standard Gaussian random variable. Other modules in the CNN are abbreviated as a loss function \(L\) parameterized by \(\omega\). Next we study the gradient estimation of \(\mathbb{E}\left[L(v;\omega)\right]\) with respect to \(\theta\).

Due to the weight sharing, the minimum unit \(\theta^{o,i}_{j,k}\) is involved in the computation of the \(o\)-th output channel. Therefore, we push \(\theta^{o,i}_{j,k}\) into the conditional density of \(v^{o}\) given \(x^{i}\), i.e.,

\[\frac{\partial\mathbb{E}\left[L(v;\omega)\right]}{\partial\theta^{o,i}_{j,k} }=\frac{\partial}{\partial\theta^{o,i}_{j,k}}\mathbb{E}\left[\int_{\mathbb{R}^ {h_{\text{out}}\times\text{out}}}\mathbb{E}\left[L(v;\omega)|\xi,x^{i}\right] \ \prod_{s=0}^{h_{\text{out}}-1}\prod_{t=0}^{w_{\text{out}}-1}\frac{1}{\sigma^{o} }\phi\left(\frac{\xi_{s,t}-\tilde{v}^{o}_{s,t}}{\sigma^{o}}\right)d\xi\right],\]

where \(\tilde{v}^{o}_{j,k}=(x*\theta)^{o}_{j,k}+b^{o}\). With the LR technique, we can finally obtain the gradient estimation for each convolutional kernel channel as below:

\[\nabla_{\theta^{o,i}}\mathbb{E}\left[L(v;\omega)\right]=\mathbb{E}\left[\frac {1}{\sigma^{o}}L(v;\omega)x^{i}*\varepsilon^{o}\right]. \tag{3}\]

The formulation (3) is quite computationally friendly and can be implemented by the convolution operator in any standard toolkit.

When the kernel tensor is small, the noise generation in (2) may be costly with respect to the parameter size. An alternative way to perturb the neuron is injecting noises to the parameters:

\[v=x*\hat{\theta}+\hat{b},\quad\hat{\theta}^{o,i}_{j,k}=\theta^{o,i}_{j,k}+ \sigma^{o}\varepsilon^{o,i}_{j,k},\quad\hat{b}^{o}=b^{o}+\sigma^{o}\delta^{o},\]

where \(\varepsilon^{o,i}_{j,k}\) and \(\delta^{o}\) are independent standard Gaussian random variables. Then, we can derive the gradient estimation by pushing \(\theta^{o,i}_{j,k}\) into the conditional density of \(\hat{\theta}^{o,i}_{j,k}\) given \(\theta^{o,i}_{j,k}\) itself, i.e.,

\[\frac{\partial\mathbb{E}\left[L(v;\omega)\right]}{\partial\theta^{o,i}_{j,k} }=\frac{\partial}{\partial\theta^{o,i}_{j,k}}\mathbb{E}\left[\int_{\mathbb{R} }\mathbb{E}\left[L(v;\omega)|\xi,\theta^{o,i}_{j,k}\right]\ \frac{1}{\sigma^{o}}\phi\left(\frac{\xi-\theta^{o,i}_{j,k}}{\sigma^{o}}\right) d\xi\right]=\mathbb{E}\left[L(v;\omega)\frac{\varepsilon^{o,i}_{j,k}}{\sigma^{o}} \right]. \tag{4}\]

Since the variances of Eq. (3) and (4) both depend on the noise dimension, the selection between the two can be determined by the total dimensions of the output and convolutional kernels.

**Recurrent neural networks**: The most widely used RNN cells can be summarized as a common structure, where hidden states \(h_{t-1}\) and inputs \(x_{t}\) are delivered to nonlinear process \(\varphi\) after a linear transformation. Some RNN variants also have an extra data flow going through \(\varphi\) without any parameter [39]. Denote the \(t\)-th input and hidden states in a generic RNN as \(x_{t}\in\mathbb{R}^{d_{x}}\) and \(h_{t}\in\mathbb{R}^{d_{h}}\), respectively. The \(t\)-th cell state and hidden state are given by

\[(h_{t},c_{t})=\varphi(u_{t},v_{t}),\quad u_{t}=\theta^{hh}h_{t-1}+b^{hh}+z_{t} ^{hh},\quad v_{t}=\theta^{xh}x_{t}+b^{xh}+z_{t}^{xh},\]

where \(\theta^{hh}\in\mathbb{R}^{d_{h}\times d_{h}}\) and \(\theta^{xh}\in\mathbb{R}^{d_{h}\times d_{x}}\) are weight matrices, \(b^{hh},b^{xh}\in\mathbb{R}^{d_{h}}\) are bias terms, \(z_{t}^{hh}=(\Sigma^{hh})^{\frac{1}{2}}\varepsilon_{t}^{hh}\), \(z_{t}^{xh}=(\Sigma^{xh})^{\frac{1}{2}}\varepsilon_{t}^{xh}\), \(\varepsilon_{t}^{hh}\) and \(\varepsilon_{t}^{xh}\) are independent standard Gaussian random vectors. Other modules in the RNN are abbreviated as a loss function \(L\) parameterized by \(\omega\). The training objective is \(\frac{1}{|\mathcal{X}|}\sum_{x\in\mathcal{X}}\mathbb{E}\left[L(h,c;\omega)\right]\), where \(h=(h_{1},\cdots,h_{T})\), and \(c=(c_{1},\cdots,c_{T})\).

Though \(\theta^{hh}\) participates in the computation of each \(u_{t}\), it is hard to make use of their joint conditional density. Thus, we first treat \(\theta^{hh}\) as different in each calculation step, and then sum up the gradient estimates. Denote the conditional density of \(u_{t}\) given \(h_{t-1}\) as \(f^{l}(\xi)\). The gradient estimation can be unrolled as

\[\nabla_{\theta^{hh}}\mathbb{E}\left[L(h,c;\omega)\right]=\sum_{t=1}^{T}\nabla_{ \theta^{hh}}\mathbb{E}\left[\int_{\mathbb{R}^{d_{h}}}\mathbb{E}\left[L(h,c; \omega)|\xi,h_{t-1}\right]f_{t}(\xi)d\xi\right].\]

In the same manner as deriving Eq. (1), we can obtain

\[\nabla_{\theta^{hh}}\mathbb{E}\left[L(h,c;\omega)\right]=\mathbb{E}\left[L(h,c; \omega)(\Sigma^{hh})^{-\frac{1}{2}}\sum_{t=1}^{T}\varepsilon_{t}^{hh}h_{t-1}^{ \top}\right]. \tag{5}\]

The gradient estimates of other parameters are similar. The formulation (5) enables the gradients of the RNN family models to be computed in parallel over input sequences without relying on BP
Through Time, which is also one of the main advantages of Transformers [40] over RNNs. Eq. (5) is derived under a generic structure, and their specific forms could be simplified for different cells. Please refer to the appendix for details.

**Graph neural networks**: The computation in different GNNs utilizes the matrix representation of graph structures and can be transformed into a combination of matrix multiplication and non-linear mappings so that it is straightforward to train GNNs by LR. Consider the graph structure \(\mathcal{G}=(\mathcal{V},\mathcal{E})\) and a graph convolutional network (GCN) [41] layer, where the input is the matrix of node features \(h\in\mathbb{R}^{|\mathcal{V}|\times d_{\text{out}}}\). GCNs first aggregate local information with degree normalization and then extract output features \(h^{\prime}\) similar to MLPs. To derive an LR estimator, we can inject noises as follows:

\[h^{\prime}=\varphi(v),\quad v=\tilde{D}^{-\frac{1}{2}}\tilde{A}\tilde{D}^{- \frac{1}{2}}h\theta+z,\]

where \(\theta\in\mathbb{R}^{d_{\text{in}}\times d_{\text{out}}}\) is the weight matrix, \(\tilde{D}\) and \(\tilde{A}\) are the graph degree and adjacency matrices with additional self-loops, \(\varphi\) is the activation function, and \(z\in\mathbb{R}^{|\mathcal{V}|\times d_{\text{out}}}\) is the noise matrix. Graph attention networks (GATs) [42] introduce the attention mechanism in to GNNs and use trainable edge weights but not the adjacency matrix in information aggregation. Given a GAT layer with single-head attention, the extracted features of the \(i\)-th node and the attention coefficient of the \((i,j)\)-th edge are given by

\[v_{i}=h_{i}\omega+\xi_{i},\quad u_{i,j}=\theta(v_{i},v_{j})+\zeta_{i,j}, \tag{6}\]

where \(\omega\in\mathbb{R}^{d_{\text{in}}\times d_{\text{out}}}\) and \(\theta\in\mathbb{R}^{2d_{\text{in}}}\) are trainable weights, \(\xi_{i}\in\mathbb{R}^{d_{\text{out}}}\) and \(\zeta_{i,j}\in\mathbb{R}\) are injected noises. The normalized attention coefficients and output features are computed as follows:

\[\alpha_{i,j}=\frac{\exp(\varphi(u_{i,j}))}{\sum_{k\in\mathcal{N}_{i}}\exp( \varphi(u_{i,k}))},\quad h^{\prime}_{i}=\varphi\big{(}\sum_{j\in\mathcal{N}_{ i}}\alpha_{i,j}v_{i}\big{)}, \tag{7}\]

where \(\mathcal{N}_{i}\) is the neighborhood of the \(i\)-th node. Eq. (7) can be regard as a non-parameterized operator, so that we only need to focus on Eq. (6) and use the push-out LR technique to estimate the gradients.

**Spiking neural networks**: Consider an \(L\)-layer SNN with leaky integrate-and-fire neurons. The time series input of \(l\)-th layer at time step \(t\) is denoted as \(x^{t,l}\in\mathbb{R}^{d_{l}}\). The next membrane potential and spike out are given by

\[u^{t+1,l+1}=ku^{t,l+1}(\mathbf{1}-x^{t,l+1})+\theta^{l}x^{t+1,l}+z^{t+1,l}, \quad x^{t+1,l+1}=I(u^{t+1,l+1},V_{\text{th}}),\]

where \(\theta^{l}\in\mathbb{R}^{d_{l}\times d_{l-1}}\) is the synaptic weight, \(k\) is the delay factor decided by the membrane time constant, \(I\) is the Heaviside neuron activation function with a threshold \(V_{\text{th}}\), \(z^{t,l}=(\Sigma^{l})^{\frac{1}{2}}\varepsilon^{t,l}\) with \(\varepsilon^{t,l}\) being standard Gaussian random variable. For all time steps, the potentials are integrated into each neuron. When the current potential passes \(V_{\text{th}}\), the neuron releases a spike signal to the next layer and resets the membrane potential to zero. The spike signal of the last layer at the last time step, namely \(x^{T,L}\), is the final output of the SNN and is passed to the loss function \(L\). Similar to handling weight sharing over time in RNNs, we can obtain \(\nabla_{\theta^{l}}\mathbb{E}[L(x^{T,L})]=\mathbb{E}[L(x^{T,L})\sum_{t=1}^{T}( \Sigma^{l})^{-\frac{1}{2}}\varepsilon^{t,l}(x^{t,l})^{\top}]\).

### Variance reduction

**Layer-wise noise injection**: LR method relies on the correlation between noise and network output to estimate the gradient. When the dimension of noises is huge, the correlation is too weak to be identified. Since LR does not require BP, it is allowed to independently perturb each layer or even each neuron and estimate the corresponding gradient, which ensures that the correlation is sufficient to produce meaningful estimates. Intuitively, layer-wise noise injection and gradient estimation can be interpreted as a method for randomized block coordinate descent [43].

**Antithetic variable**: Note that Eq. (1) and other LR estimators calculate the gradient estimate through a Monte Carlo (MC) integration. A simple but effective variance reduction method we can use is the antithetic variates method [44]. Assume we repeatedly evaluate the integrand \(2n\) times to estimate the gradient. We can generate injected noises for the first \(n\) evaluations, and then multiply the noises by \(-1\) for the subsequent \(n\) evaluations.

**Quasi-Monte Carlo**: Quasi-Monte Carlo (QMC) [45] generates noise and evaluates the integrand using low-discrepancy sequences, which can fill the integration region better than the pseudo-random
sequence used in standard scientific computing toolkits. [46] points out that QMC has better performance than MC when the integrand is smooth and the dimension is relatively small. Given the smoothness of neural networks and the layer-wise perturbation mentioned earlier, it is desirable to try QMC in LR estimation.

### Relationships between LR and other methods

**LR and other perturbation-based methods**: Consider a neural module with noise injected into its parameter, i.e., \(v=\varphi(x;\theta+z)\), where \(x\) and \(v\) is the input and output, \(\varphi\) is a non-linear operator, and the noise \(z\) follows a density \(f\) defined on \(\Omega\). Denote the rest of the neural network as a loss function \(L\) with parameter \(\omega\). We can obtain the LR estimator by pushing \(\theta\) into the conditional density of \(z\) as below:

\[\nabla_{\theta}\mathbb{E}\left[L(v;\omega)\right]=\nabla_{\theta}\mathbb{E} \left[\int_{\Omega}\mathbb{E}\left[L(v;\omega)|\xi,\theta\right]f(\xi-\theta )d\xi\right]=\mathbb{E}\left[-L(v,\omega)\frac{\nabla_{z}f(z)}{f(z)}\right],\]

which incorporates several previous works [47; 48; 49; 24] into the framework of push-out LR method.

**LR and reinforcement learning**: We first show that the key theorem in reinforcement learning (RL) can be derived by the push-out LR method. The RL agent is modeled as a parameterized probability distribution \(\pi_{\theta}\), and \(p\) represents the transition kernel of the environment. Denote the action and state at the \(t\)-th time step as \(a_{t}\) and \(s_{t}\), where \(a_{t}\sim\pi_{\theta}(\cdot|s_{t})\), \(s_{t+1}\sim p(\cdot|s_{t},a_{t})\). The reward function is denoted as \(r\), and the cumulative reward is given by \(R(s,a)=\sum_{t=0}^{T-1}r(s_{t+1},a_{t},s_{t})\), which is a function of the state and action sequence. Policy-based RL aims at solving the optimization problem: \(\max_{\theta}\mathbb{E}\left[R(s,a)\right]\). In the same manner as deriving Eq. (5), we can treat \(\theta\) as different at each step and push it into the conditional density of \(a_{t}\) given \(s_{t}\) to perform the LR technique, i.e.,

\[\nabla_{\theta}\mathbb{E}\left[R(s,a)\right]=\sum_{t=0}^{T-1}\nabla_{\theta} \mathbb{E}\left[\int_{\mathcal{A}}\mathbb{E}\left[R(s,a)|\xi,s_{t}\right]\pi_{ \theta}(\xi|s_{t})d\xi\right]=\mathbb{E}\left[R(s,a)\sum_{t=0}^{T-1}\nabla_{ \theta}\ln\pi_{\theta}(a_{t}|s_{t})\right],\]

which results in the well-known policy gradient theorem [50].

Not only can we derive the RL gradient using the push-out LR method, but the reverse is also possible. Consider the LR-based RNN training discussed earlier, where the input \(x_{t}\) and hidden state \(h_{t-1}\) can be viewed as the RL state, and the results of the linear transformation \(u_{t}\) and \(v_{t}\) are viewed as the actions. The parameterized part and the rest in the RNN cell correspond to the RL agent and the simulation environment, respectively, and Eq. (5) can be modified directly from the policy gradient theorem.

Therefore, RL and deep learning can be unified from the perspective of the push-out LR framework. RL implicitly leverages the advantages of LR to overcome the problem of a black-box simulation environment that affects gradient backward propagation.

## 4 Evaluations

### Experiment settings

Our study focuses on classification tasks with different modalities using various neural network architectures. For image classification, we experiment with ResNet-5 and VGG-8 on the Cifar-10 dataset [28] and SNN on the MNIST [31] and Fashion-MNIST [32] datasets. We use RNNs (RNN, GRU [51], and LSTM [39]) to classify the articles on the Ag-News dataset [29] and use GNNs (GCN and GAT) to classify the graph nodes on the Cora dataset [30]. We compare our method to two baselines (BP and ES) using LR with noise injection on logits (LR-L) and weight/logits (LR-WL) in a hybrid manner. Supplementary materials provide network structures and training details.

### Numerical results on convolutional neural networks

**Evaluation criteria**: We evaluate all the methods for CNNs on the following performance criteria: 1) Task performance: the classification accuracy on the benign samples (Ori.); 2) Natural noise robustness: the classification accuracy on the natural noise corrupted dataset, where we adopt the Gaussian, uniform, and Poisson noises; 3) Adversarial robustness: the classification accuracy on
the adversarial examples, which are crafted by FGSM, I-FGSM, and MI-FGSM; 4) distribution shift robustness: the classification accuracy on the corrupted dataset with different distribution corresponding to the original dataset, where we transform the colored dataset into the grey images and apply the random mask to images respectively for the construction of OOD datasets. For the adversarial attacks, we set the maximum perturbation for each pixel as \(8/255\), and the maximum number of iterations as \(5\) for I-FGSM and MI-FGSM. For all the evaluations, we consider the corruption on the whole test set.

**Results**: As shown in Tab. 1, LR-WL can achieve comparable performance to BP with an improvement of \(0.4\%\) for ResNet-5 and a minor drop of \(2.6\%\) for VGG-8. Meanwhile, LR-WL has a significant improvement in the adversarial robustness of \(8.1\%\) on average. In ResNet-5, the classification accuracy relies more on the training of the later layers, where adding noise to logits incurs less cost, thus LR-L outperforms ES and achieves similar performance as LR-WL. In VGG-8, the situation is exactly the opposite. And the success of LR-WL in CNNs comes from combining the strengths of both noise injection modes.

### Numerical results on recurrent neural networks

**Evaluation criteria**: We evaluate all the methods for RNNs on task performance and robustness. For the robustness evaluation, we adopt two corruptions, 1) RanMask, in which we randomly mask a fixed ratio of \(90\%\) of the words in the whole article; 2) Shuffle, which random shuffles words in sentences; 3) adversarial attacks, including PWWS and GA. We set the maximum ratio of word substitutions as \(25\%\). Due to the low efficiency of adversarial text attacks, we consider \(200\) examples for evaluation on adversarial attacks and \(1,000\) for other perturbations.

**Results**: The evaluation results of the RNN family can be shown in Tab. 2. Due to the gradient vanishes problem in RNN for long sentences, the conventional BP method can not train the neural network well and it only achieves an accuracy of \(64.7\%\). Without relying on the chain-rule computation, the LR method can mitigate the gradient vanishes problem and improve the performance of RNN significantly by \(5.5\%\). For training GRU and LSTM, the LR method can achieve comparable performance on the clean dataset. On the other hand, neural network training using the LR method can achieve better robustness compared with the BP, namely \(2.3\%\) under the RanMask corruption, \(2.63\%\) under the Shuffle corruption, and \(4.0\%\) under the PWWS attack.

\begin{table}
\begin{tabular}{c c c c c c c c c c c} \hline \hline \multirow{2}{*}{Models} & \multirow{2}{*}{Alg.} & \multirow{2}{*}{Ori.} & \multicolumn{4}{c}{Natural Noise} & \multicolumn{4}{c}{Adversarial Noise} & \multicolumn{2}{c}{Distribution shift} \\ \cline{4-11}  & & & \multicolumn{2}{c}{Gaussian} & Uniform & Poisson & FGSM & I-FGSM & MI-FGSM & Grey & RanMask \\ \hline \multirow{4}{*}{ResNet-5} & BP & 64.4 & 59.4 & 60.7 & 44.2 & 12.1 & 1.1 & 1.7 & 53.5 & 54.1 \\  & ES & 35.1 & 35.1 & 35.0 & 29.5 & 20.8 & 4.1 & 4.9 & 25.3 & 34.3 \\  & LR-L & 62.3 & 61.4 & 61.9 & 50.9 & 26.0 & 6.4 & 9.4 & 55.9 & 55.5 \\  & LR-WL & 64.8 & 64.2 & 64.6 & 55.3 & 26.2 & 6.7 & 9.4 & 58.2 & 59.1 \\ \hline \multirow{4}{*}{VGG-8} & BP & 79.8 & 69.4 & 76.4 & 37.8 & 32.4 & 1.2 & 0.9 & 67.8 & 61.9 \\  & ES & 65.8 & 63.7 & 65.4 & 35.6 & 25.7 & 2.7 & 2.7 & 61.9 & 60.5 \\ \cline{1-1}  & LR-L & 58.7 & 50.6 & 56.2 & 30.7 & 27.9 & 11.5 & 10.5 & 48.9 & 47.6 \\ \cline{1-1}  & LR-WL & 77.3 & 72.7 & 76.7 & 42.9 & 34.6 & 2.8 & 3.5 & 71.4 & 64.9 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Results of ResNet-5 and VGG-8.

\begin{table}
\begin{tabular}{c c c c c c c c c c} \hline \hline \multirow{2}{*}{Attack} & \multicolumn{2}{c}{RNN} & \multicolumn{2}{c}{GRU} & \multicolumn{2}{c}{LSTM} \\ \cline{2-11}  & BP & ES & LR & BP & ES & LR & BP & ES & LR \\ \hline Ori. & 64.4 & 65.9 & 70.9 & 88.2 & 84.6 & 88.9 & 89.4 & 85.7 & 89.4 \\ RanMask & 53.4 & 53.4 & 56.7 & 63.9 & 59.2 & 65.5 & 60.1 & 56.2 & 62.1 \\ Shuffle & 62.7 & 63.2 & 69.0 & 87.1 & 80.2 & 88.3 & 88.9 & 82.6 & 89.3 \\ PWWS & 50.0 & 36.5 & 54.0 & 58.5 & 32.0 & 63.0 & 59.0 & 31.5 & 62.5 \\ GA & 37.0 & 23.5 & 40.5 & 56.5 & 20.0 & 60.5 & 59.5 & 21.0 & 62.0 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Results of the RNN family.
### Numerical results on graph neural networks

**Evaluation criteria**: We evaluate all the methods for GNNs on task performance and robustness. For the robustness evaluation, we adopt two adversarial attacks, including R.A. and DICE. For both the R.A. and DICE attacks, we perturb the edge ratio from \(10\%\) to \(100\%\) in \(10\%\) increments compared to the original number of edges and report the average results. For all the evaluations, we consider the corruption on the whole test set.

**Results**: As shown in Tab. 3, compared with ES, LR achieves a comparable performance relative to BP, with a performance gap of \(4.8\%\) for GCN and \(0.4\%\) for GAT on the clean dataset. LR further brings an average of \(4.3\%\) improvement for GCN and an average of \(5.1\%\) improvement for GAT on the model robustness. Moreover, the GAT training by LR does not require any additional copy of data, which has a similar computational cost compared to BP.

### Numerical results on spiking neural networks

**Evaluation criteria**: We evaluate the performance of different optimization methods by task performance and adversarial robustness. We use FGSM, I-FGSM, and MI-FGSM to craft the adversarial examples on MNIST and Fashion-MNSIT datasets, where the maximum perturbation for each pixel is \(0.1\) for all of the attack methods, the max number of iterations is \(15\), and the step size is \(0.01\) for I-FGSM and MI-FGSM. For all the evaluations, we consider the corruption on the whole test set.

**Results**: As shown in Tab. 4, ES fails in training SNN, and the LR method can achieve a higher classification accuracy on clean and adversarial datasets. Specifically, LR has an average improvement of \(1.55\%\) on two clean datasets and \(44.3\%\) on adversarial examples. The activation function in SNN is discontinuous, which hinders the application of the chain rule-based BP method, and the BP adopted here is based on an approximate computation. The unbiasedness of the LR gradient estimation leads to superior results.

### Ablation study

We study on the effect of the perturbation on different parts of neurons, the initialization for noise magnitudes \(\sigma\), and the variance reduction techniques. In our experiments, we use MC sampling to generate the noise for experiments except when we claim otherwise. We conduct experiments using ResNet-5 on the Cifar-\(10\) dataset for the ablation study. We use our LR method to estimate the gradient and compare it with that calculated by BP. We compute the cosine similarity to evaluate the performance. A cosine similarity closer to 1 means a more accurate gradient estimation.

**Impact of the perturbation on weights and logits**: We present the gradient estimation performance when perturbing the logits only in Fig. 1 (a), the weights only in Fig. 1 (b), weights and logits in a hybrid manner in Fig. 1 (c). Compared to Fig. 1 (b), Fig. 1 (a) shows that the estimated gradient of the last several layers, including the last two convolutional layers and the linear layer, is more accurate with the same number of copies when perturbing the weights rather than the logits. On the other hand, when perturbing the weights, though the estimated gradients for the linear layer is less accurate, it has a higher cosine similarity with that calculated by BP for the first several layers, including the first and second convolutional layers. Motivated by this phenomenon, we hybrid these two perturbing strategies, where we inject the noise on weights for the first two layers and logits for others, respectively. As shown in Fig. 1 (c), the performance of gradient estimation for all five layers has a significant improvement compared to purely adding perturbation on weights or logits.

**Impact of the initialization for \(\sigma\)**: Following the design of the aforementioned hybrid noise injection mode, we further study the effect of the initialization for \(\sigma\). From left to right in Fig. 1 (d)-(f), we

\begin{table}
\begin{tabular}{c c c c c c c} \hline \hline \multirow{2}{*}{Attack} & \multicolumn{3}{c}{GCN} & \multicolumn{3}{c}{GAT} \\ \cline{2-7}  & BP & ES & LR & BP & ES & LR \\ \hline Ori. & 80.4 & 57.3 & 75.6 & 82.2 & 34.4 & 81.8 \\ R.A. & 62.5 & 48.2 & 68.9 & 70.4 & 30.1 & 77.3 \\ Dice & 63.2 & 45.7 & 65.4 & 71.3 & 28.9 & 74.6 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Results of the GNN family.

\begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline \multirow{2}{*}{Attack} & \multicolumn{3}{c}{GCN} & \multicolumn{3}{c}{GAT} \\ \cline{2-7}  & BP & ES & LR & BP & ES & LR \\ \hline Ori. & 80.4 & 57.3 & 75.6 & 82.2 & 34.4 & 81.8 \\ R.A. & 62.5 & 48.2 & 68.9 & 70.4 & 30.1 & 77.3 \\ Dice & 63.2 & 45.7 & 65.4 & 71.3 & 28.9 & 74.6 \\ \hline \hline \end{tabular}
\end{table}
Table 4: Results of SNN.
perform LR with \(\sigma\) setting to \(1\times 10^{-1}\), \(1\times 10^{-2}\), and \(1\times 10^{-3}\), respectively. The optimal choice of \(\sigma\) varies for different noise injection modes. In this network architecture, a larger \(\sigma\) is more suitable for gradient estimation using noise injection on logits, whereas a smaller one is better when injecting noise on weights. Improper initialization of \(\sigma\) can result in low cosine similarities between the estimated gradient by LR and that calculated by BP, which would lead to reduced training efficiency.

**Impact of the variance-reduction techniques**: In Fig. 1 (g)-(i), we study the effect of three proposed variance-reduction techniques. Without the layer-wise noise injection, it is difficult to identify the correlation between the noise in each layer and the loss evaluation, which also makes it challenging to select appropriate initial \(\sigma\) for each layer. As shown in Fig. 1 (g), it requires a much larger number of copies to achieve comparable performance on the first two and last layers, and fails to efficiently estimate the gradients of the rest layers. From Fig. 1 (h), it can be seen that without the antithetic variable, there are low similarities between the estimated gradients by LR and that by BP for all convolutional layers even with a large number of copies. In Fig. 1 (i), we apply QMC to generate the noise for gradient estimation. Although the results are close to MC, the neuron-level parallelism is allowed by LR, resulting in a lower dimension of MC integration, which implies that QMC may lead to a more significant improvement.

## 5 Conclusion

In our work, we extend the LR method for gradient estimation to a broader range of neural network architectures, including CNNs, RNNs, GNNs, and SNNs. LR completely removes the need for BP during training, thus eliminating the dependency on recursive gradient computations. With variance reduction techniques, our method achieves comparable performance to the BP method on clean data and exhibits improved robustness, particularly in the presence of adversarial attacks. Moreover, we

Figure 1: Ablation study on the effect of the perturbation on different parts of neurons, the initialization for noise magnitude \(\sigma\), and the variance reduction technique.
also discuss the relationship between LR and other perturbation-based methods and unify RL and deep learning from a bidirectional perspective of the push-out LR framework. We believe that our method has the potential to be a valuable addition to the arsenal of methods available for training neural networks and look forward to further exploring its applications in future work.
