\subsubsection*{\bibname}}

\bibliographystyle{apalike}

\begin{document}

\twocolumn[

\aistatstitle{Adaptive Real Time Exploration and Optimization for Safety-Critical Systems}

\aistatsauthor{Buse Sibel Korkmaz \And Mehmet Mercangöz \And Marta Zagórowska}

\aistatsaddress{ Imperial College London \And Imperial College London \And ETH Zürich} ]

\begin{abstract}
We consider the problem of decision-making under uncertainty in an environment with safety constraints. Many business and industrial applications rely on real-time optimization with changing inputs to improve key performance indicators. In the case of unknown environmental characteristics, real-time optimization becomes challenging, particularly for the satisfaction of safety constraints. We propose the ARTEO algorithm, where we cast multi-armed bandits as a mathematical programming problem subject to safety constraints and learn the environmental characteristics through changes in optimization inputs and through exploration. We quantify the uncertainty in unknown characteristics by using Gaussian processes and incorporate it into the utility function as a contribution which drives exploration. We adaptively control the size of this contribution using a heuristic in accordance with the requirements of the environment. We guarantee the safety of our algorithm with a high probability through confidence bounds constructed under the regularity assumptions of Gaussian processes. Compared to existing safe-learning approaches, our algorithm does not require an exclusive exploration phase and follows the optimization goals even in the explored points, which makes it suitable for safety-critical systems. We demonstrate the safety and efficiency of our approach with two experiments: an industrial process and an online bid optimization benchmark problem.
\end{abstract}

\section{INTRODUCTION}
In most sequential decision-making problems under uncertainty, there exists an unknown function with noisy feedback that the decision-maker algorithm needs to sequentially estimate and optimize to reveal the decisions leading to the highest reward. Each decision incurs an instantaneous reward with an initially unknown distribution in a stochastic optimization setting. In this setting, initial decisions are made based on some heuristics, and incurred reward is memorized to exploit for the next decisions. Therefore, the uncertainty due to unknown characteristics decreases while new decisions are made based on previous reward observations. Even though more exploration can help to optimize the decisions by revealing more information in each iteration, it could be very expensive to evaluate the unknown function for many applications. Therefore, there is a need to balance exploring unknown decision points and exploiting previous experiences. The trade-off between exploration and exploitation is extensively studied in the literature, and the multi-armed bandit (MAB) approaches with confidence bounds are suggested to solve this problem by optimally balancing exploration and exploitation \citep{bubeck2012regret}.

\subsubsubsection\textbf{{Related work.}}\hspace{11pt}The idea of using confidence bounds to balance the exploration-exploitation trade-off through the optimism principle first appeared in \citet{1985Lai} with the utilization of the upper confidence bound (UCB). Since then, it led to the development of UCB algorithms for stochastic bandits with many arms \citep{banditbook}. Many efficient algorithms build for bandit problems having a cost or reward function under certain regularity conditions \citep{Dani2008, Bubeck2008}. In \citet{Srinavas2010}, the authors divided the stochastic optimization problem into two objectives: (1) unknown function estimation from noisy observations, and (2) optimization of the estimated function over the decision set. They used kernel methods and Gaussian processes (GPs) to model the reward function since these methods encode the regularity assumptions through kernels \citep{RasmussenGP}. Even though these methods are very successful to model and optimize unknown reward functions, they are not applicable to the specific setting that we are interested in this paper since they do not consider safety constraints.

In safety-critical systems, it is not possible to do exploration in some parts of the decision space due to safety concerns, which are modelled as safety constraints in optimization. Many industrial applications fall under the safety-critical systems due to their risk of  danger to human life, leading to substantial economic loss, or causing severe environmental damage. For example, we can consider a chemical process plant as a safety-critical operation since we need to satisfy the constraints of chemical reactions and surrounding processes to not cause a hazardous event for the environment and human operators. Therefore, if we want to approximate the unknown characteristics and optimize a process in this plant, we need to utilize safe exploration algorithms which allow only the exploration of feasible decision points by enforcing safety constraints. In this work, we define a feasible decision point in the given decision set as satisfying the safety and other constraints of the given problem.

The safe exploration problem has been studied by formalizing it as both bandit and Markov Decision Processes (MDPs) problems \citep{Schreiter2015, Turchetta2019, Wachi2018,Turchetta2016}. In \citet{Sui2015}, the authors have introduced the SafeOpt algorithm and showed it is possible that safely optimize a function with an unknown functional form by creating and expanding a safe decision set through safe exploration under certain assumptions. They also provided the safety guarantee of SafeOpt by using the confidence bounds construction method of \citet{Srinavas2010}. These safe exploration algorithms are applied to many control and reinforcement learning problems and proved their success in those domains \citep{app1,app2}. However, these algorithms either require an exploration phase or apply a trade-off strategy between optimal decisions and exploration. There are some real-world applications such as industrial processes which cannot afford an exclusive exploration phase or require following optimization goals of the system even in the explored points. For example, any deviation from target satisfaction in an industrial process may cause high costs to the plant or damage the reputation of the responsible party. Thus, there is a need to consider the exploration in an adaptive manner in accordance with the requirements of the environment and this aspect is not covered by existing safe exploration algorithms.

\subsubsubsection\textbf{{Our contributions.}}\hspace{11pt}In this paper, we propose a novel safe exploration algorithm, ARTEO, where we cast multi-armed bandits as a mathematical programming problem subject to safety constraints for the optimization of safety-critical systems. Our contribution is posing a further constraint to the safe exploration by instructing the algorithm to avoid decisions not satisfying the optimization goals of the safety application. This constraint discourages exploration which makes it harder to learn unknown rewards in the MAB setting. To compensate for that, we incorporate the exploration as a contribution to the utility function in an adaptive manner by encouraging explicitly to take decisions at points with high uncertainty which is quantified by the covariance of GP models. We adaptively control the size of this contribution using a heuristic in accordance with the requirements of the environment. ARTEO uses GPs to model how the unknown parameters of an application change based on external factors such as inputs and optimization goals, make new decisions when a change in the optimization goals or in the given input is occurred and explores new decision points with either a change in the optimization goals or based on uncertainty. We establish the safety of ARTEO by constructing the confidence bounds as \citet{Srinavas2010} and showing certain assumptions hold in our setting. We demonstrate the safety and efficiency of our algorithm using two real-world examples: the power management problem of the food processing plant refrigeration system and the high-dimensional online bid optimization problem. In the first experiment, the results show that we are able to successfully control three key performance indicators of the process by applying ARTEO. In the second experiment, our approach satisfies the constraints while yielding more profitable outcomes.

\section{PROBLEM STATEMENT AND BACKGROUND}

We want to find a sequence of decisions, $x_1,x_2, \ldots, x_T$, such that a certain utility function $f$ is minimised. At every iteration $t$, $t=1,\ldots,T$, the cost function depends on the decision $x_t \in D$ and on the system characteristics $p_t\in P$. The characteristics $p_t$ depends on the decision $x_t$, i.e. $p_t=h(x_t)$ where $p:D^\mathbb{N} \rightarrow  P^\mathbb{N}$ with $N$ arms. After making a decision $x_t$, we obtain a noisy measurement, $y_t=f(x_t,p(x_t))+\epsilon$ where $f:D^\mathbb{N} \times P^\mathbb{N} \rightarrow \mathbb{R}$ and $\epsilon$ is $R$-sub-Gaussian noise for a fixed constant $R \geq 0$ \citep{SubGaussian}. We assume that we know the functional form of $f(\cdot)$, but the functional form of $p(\cdot)$ is unknown. Furthermore, at every iteration, the decision $x_t$ must satisfy the  constraints $g_k(x_t,p(x_t))+ h_k\leq 0$, where $k=1,2,\ldots, K$ with $K$ denoting the number of constraints. The value of $h_k\geq 0$ is called a safety threshold. Thus, we can formalise our optimisation problem at time $t$ as:
\begin{equation}
 X_t=\argmin_{X} f(X,p(X)) \ s. \ t. \ g_k(X,p(X))+h_k\leq 0,  \forall k
\end{equation}

where $X = [x_{i}, ..., x_{N}]$ for arms ${i=1, \dots, N}$. If we know $p(x_t)$, we can solve the problem as an optimization problem with noise. For instance, we could use the concept of real-time optimization (RTO) from the process control domain \citet{RTOChemical} and use the approach proposed by \citet{SafeRTO}.
However, the characteristics $p(.)$ are unknown and need to be estimated. Thus, at every iteration $t$, we need to first solve an estimation problem to find $p$, then solve the optimization problem Equation (1). Solving an optimization problem by combining estimation then optimization is a common approach \citep{EstimationforRTO1, EstimationforRTO2}. However, few approaches quantify the uncertainty inherent in the estimation of $p$. In the current paper, we propose to estimate $p$ using Gaussian processes and use their properties to ensure safety.

\subsubsubsection{\textbf{Gaussian processes.}}\hspace{11pt}Gaussian processes are non-parametric models which can be used for regression. GPs are fully specified by a mean function $\mu(x)$ and a kernel  $k(x,x')$ which is a covariance function and decides the shape of prior and posterior in GP regression \citep{RasmussenGP}. The goal is to predict the value of the utility function $f$ including unknown components at the decision point $x^*$ by using GPs. Assuming having a zero mean prior, the posterior $p(f(x^*)|x^*, D)$ follows $N(\mu(x^*), \sigma^2(x^*))$ that satisfy,
\begin{equation}
\begin{aligned}
  \mu(x^*) = K_{N*}^T(K_{NN} + \sigma^2 I)^{-1}(y_1,...,y_N)^T\\
  \sigma^2(x^*) = k(x^*,x^*) - K_{N*}^T(K_{NN} + \sigma^2 I)^{-1}
\end{aligned}
\end{equation}

where $i,j \in \left\{1,.., N\right\}$ and denotes the index of observations, $[K_{N*}]_i = k(x_i,x^*)$, $K_{NN}$ is the $N\times N$ positive definite kernel matrix with  $[K_{NN}]_{i,j} = k(x_i, x_j)$ and $y_i$ is the noisy feedback of $i$th observation.

\subsubsubsection{\textbf{Regularity assumptions.}}\hspace{11pt}We do not have any prior knowledge of how unknown parameters in $f$ change based on external factors such as the optimization goals or inputs, and to provide safety with high probability at decision points and in exploration we need to make some assumptions, which are common in the literature \citep{Srinavas2010, Sui2015, Sui2018, app1}. We assume the decision set $D$ is compact (and finite in our case) as being a closed and bounded subset of Euclidean space \citep{Compactness}. Furthermore, we assume the utility function $f$ is a continuous function on the compact set $D$ such that it is possible to attain the minimum and maximum values for $f$ due to its smoothness \citep{Continuity}. The smoothness assumption holds for utility function $f$ since we sample it from a GP with a positive definite kernel \citep{RasmussenGP}.

Gaussian processes are related to reproducing kernel Hilbert space (RKHS) notion through its positive semidefinite kernel functions \citep{kernelcharacteristics}, which help us to construct confidence bounds in a safe manner later. The RKHS which is denoted by $\mathcal{H}_k(D)$ is formed by ``nice functions'' in a complete subspace of $L_2(D)$ and the inner product $\langle\cdot,\cdot\rangle_k$ of functions in RKHS follows the reproducing property: $\langle f,k(x,\cdot)\rangle_k = f(x)$ for all $f \in \mathcal{H}_k(D)$. The smoothness of a function in RKHS with respect to kernel function $k$ is measured by its RKHS norm $\|f\|_k = \sqrt{\langle f,f\rangle_k}$ and for all functions in $\mathcal{H}_k(D)$ $\|f\|_k < \infty$ \citep{Scholkopf2002}. Thus, we assume a known bound $B$ for the RKHS norm of the unknown function $f$: $\|f\|_k < B$. We use this bound $B$ to control the width of the confidence interval later in Equation (4).

\subsubsubsection{\textbf{Confidence bounds.}}\hspace{11pt}In ARTEO, we give safety constraints to the solver as hard constraints. The solver uses the bounds of the confidence interval which is constructed by using the standard deviation of conditioned GP on previous observations to decide the feasibility of a chosen point $x_t$. Hence, the correct classification of decision points in $D$ relies on the confidence-bound choice. Under the regularity assumptions stated before, Theorem 3 of \citet{Srinavas2010} and Theorem 2 of \citet{BoundProof} proved that it is possible to construct confidence bounds which include the true function with probability at least $1-\delta$ where $\delta \in (0,1)$ on a kernelized multi-armed bandit problem setting with no constraints. Moreover, as shown by \citet{Sui2015} in Theorem 1, this theorem is applicable to multi-armed bandit problems with safety constraints. Hence, we can state that the probability of the true value of safety function $g$ at the decision point $x$ is included by the confidence bounds in iteration $t$:
\begin{equation}
P\bigr[|g(x_t) - \mu_{t-1}(x_t)| \leq  \sqrt{\beta_t}\sigma_{t-1}(x_t)\bigr] \leq 1 - \delta, \ t \geq 1
\end{equation}

where $\mu_{t-1}(x)$ and $\sigma_{t-1}(x)$ denote the mean and the standard deviation of sampled functions at $x_t$ from a GP at iteration $t$, which is conditioned on previous $t-1$ observations to obtain the posterior. $\delta$ is a parameter that represents the failure probability in Equation (3) and $\beta_t$ controls the width of the confidence interval and satisfies Equation (3) when:

\begin{equation}
\sqrt{\beta_t} = B + R\sqrt{2(\gamma_{t-1} + 1 + \text{ln}(1/\delta))}
\end{equation}

where the noise in observations is $R$-sub-Gaussian and $\gamma_{t-1}$ represents the maximum information gain after $t-1$ iterations and it is formulated as:

\begin{equation}
\gamma_{t-1} = \max_{|A| \leq t-1} I(f; y_A)
\end{equation}

\section{ARTEO ALGORITHM}

We develop the ARTEO algorithm for safety-critical environments with high exploration costs. Therefore, the ARTEO algorithm is triggered by a change in the optimization goal or by a new event. At each iteration, the algorithm updates the posterior distributions of GPs with previous noisy observations as in Equation (2) and provides an optimized solution for the desired outcome based on how GPs model the unknown components. It does not require a separate training phase, instead, it learns during normal operation. The details of the safe learning and optimization are given next.

\subsection{Safe Learning}

In a setup with multi-arms, the decision set $D_i$ is defined for each arm $i$ as satisfying the introduced assumptions in Section 2. For each arm, a GP prior and initial ``safe seed'' set is introduced to the algorithm. The safe seed set $S_0$ includes at least one safe decision point with the true value of the safety function at that point. As in many published safe learning algorithms \citep{Sui2015, Sui2018,Turchetta2019}, without a safe seed set, an accurate assessment of the feasibility of any points is not possible. Each iteration of the algorithm could be triggered by time or an event. After receiving the trigger, the algorithm utilizes the past noisy observations to obtain the GP posterior of each arm to use in the optimization of the utility function, which includes the cost of decision and uncertainty. For the first iteration, past observations are given as safe seed sets.

The uncertainty in the utility function $f_t$ is quantified as:
\begin{equation}
    U_t(X_t, p(X_t)) = \sum_{i}^{n}{\sigma(x_{it})}
\end{equation}
where $\sigma(x_{it})$ is the standard deviation of $GP_i$ at the point $x_{it}$ for the arm $i$ in the iteration $t$. It is incorporated into the utility function $f$ by multiplying by an adjustable parameter $z$ as next:
\begin{equation}
    f_t(X_t, p(X_t)) = C_t(X_t, p(X_t)) + z U_t(X_t, p(X_t))
\end{equation}
where $C_t(X_t, p(X_t))$ represents the cost of decision at the evaluated points. In our experiments, the priority of the algorithm is optimizing the cost of decisions under given constraints. Hence, the uncertainty weight remains zero until the environment becomes available for exploration. Until that time, the algorithm follows optimization goals and learns through changes in the optimization goals such as operating in a different decision point to satisfy a new demand or a change in the given input such as starting a new advertising campaign in our second experiment.
\begin{figure}[t]
\vspace{-1em}
\begin{algorithm}[H]
\caption{ARTEO}\label{alg:cap}
\begin{algorithmic}
\Require Decision set $D_i$ for each arm $i \in \left\{1,..,n\right\}$
\\ \hspace{7mm} GP priors for each ${GP}_i$
\\ \hspace{7mm} Safe seed set for each GP as $ S_{i,0}$
\\ \hspace{7mm} Uncertainty weight as $\zeta$
\\ \hspace{7mm} Utility function $f$
\\ \hspace{7mm} Safety function $g$
\\ \hspace{7mm} Safety threshold $h$
\\ \hspace{7mm} Confidence interval width parameter $\beta_t$
\For{$t=1,...,T_0$}
\If{Steady-State = True}
\If{ the exploration heuristic rule holds}
\State $z \gets \zeta$
\Else{}
\State $z \gets 0$
\EndIf
\For{$i=1,..., n$}
\State Update ${GP}_i$ by conditioning on $S_{i,t-1}$
\EndFor
\State $x_{\left\{1, ..., n\right\},t}^* \gets \operatorname*{argmin}_{x_i \in D_i}{f\ \ s. \ t. \ \ g }$
\For{$i=1,..., n$}
\State $y_{i,t}^f \gets$ $f(x_{i,t}^*, p(x_{i,t}^*)) + \epsilon_t^f$
\State $y_{i,t}^g \gets$ $g(x_{i,t}^*, p(x_{i,t}^*)) + \epsilon_t^g$
\State $S_{i,t} \gets $ $S_{i,t-1} \cup {y_i,t} $
\EndFor
\EndIf
\EndFor
\end{algorithmic}
\end{algorithm}
\vspace{-3em}
\end{figure}

The exploration is controlled by the $z$ hyperparameter. It is possible to create different heuristic rules for changing the $z$ parameter during the execution to fine-tune it based on the needs of the simulated environment. This change in the utility function encourages the real-time optimizer to take decisions on unexplored points to decrease uncertainty as a part of the optimization. The safety constraints and the aim of satisfying the expectation with a minimum cost are still under consideration by the RTO during this phase, so, the algorithm does not violate safety constraints with a $1 - \delta$ probability as explained in the previous section.

Since we apply the ARTEO algorithm to a dynamical system simulation in our first experiment, we build a steady-state detector and evaluate costs and constraints in a steady-state. For the steady-state detection algorithm, we use a statistical test-based detector as common in many real-time optimization tools \citep{SteadyStateRTO}. The ARTEO algorithm is given in Algorithm 1.

\subsubsubsection{\textbf{Discussion.}} We want to discuss the significance and alternative implementations of parameter $z$ for the exploration. In a fixed environment with a fixed $z$, the contribution of the exploration could decrease as the environment is explored. As an extension to our current implementation, adding a forgetting factor can be considered to recover the excitement for exploration. For an environment with dynamically changing reward distributions of arms over time, the exploration can be driven by an adaptive version of $z$ where the distance between the prediction of GPs and the noisy feedback from the system is monitored as a metric for adapting (or adjusting) the exploration effort.

\subsection{Optimization}

The RTO incorporates the posterior of the Gaussian process into decision-making by modelling the unknown parameters for each arm by using the mean and standard deviation of GPs . The utility function is the objective function in the RTO formulation and the safety thresholds are constraints. In the utility function, the mean of GP posterior of each arm is used to evaluate the cost of decision and the standard deviation of GPs is used to measure uncertainty as in Equation (6). In the safety function, the standard deviation of the GP posterior of each arm is used to construct confidence bounds and then these bounds are used to assess the feasibility of evaluated points. The optimizer solves the minimization problem under safety constraints within the defined decision set of each arm. Any optimization algorithm that could solve the given problem can be used in this phase. We have constrained nonlinear problems in the experiments section, and we choose interior-point and sequential-least square programming (SLSQP) algorithms to solve our first and second problems, respectively.

\subsection{Complexity}

In each iteration, the ARTEO updates GP models by computing marginal likelihoods, and finds a feasible solution within constraints. The overall time complexity of each run of the algorithm is the number of iterations $t$ times the time complexity of each iteration. The first computationally demanding step in the ARTEO is fitting GPs on safe sets. The time complexity of training a full GP, i.e. exact inference, is $\mathcal{O}(n^3)$ due to the matrix inversion where $n$ is the number of observations \citep{GPcovariancecomplexity}. It is possible to reduce it further by using low-rank approximations which is not in the scope of our work \citep{GPlowrank1, GPlowrank2}. We introduce an individual GP for each arm, so the total complexity of GP calculations is $\mathcal{O}(kn^3)$, where $k$ is the number of arms.

The next demanding step is optimization. The computational demand of RTO depends on chosen optimization algorithm and the required $m$ number of steps to converge. The most computationally expensive step for all iterations in both mentioned optimization algorithms is the LDL factorization of a matrix with a $\mathcal{O}(n^3)$ complexity \citep{interiorpointmethod, SLSQP}. Hence, the complexity of RTO becomes $\mathcal{O}(mn^3)$. In our implementation, $m \gg k$, so, the time complexity of each iteration in ARTEO scales with the RTO complexity. Therefore, the overall time complexity of one iteration of ARTEO is $\mathcal{O}(mn^3)$. The memory complexity of the algorithm is $\mathcal{O}(kn^2)$, which is dominated by matrix storage in GPs and optimization.

\section{EXPERIMENTS}

In this section, we evaluate our approach on two real-life applications: a fish-freezing process plant refrigeration system and online bid optimization. The former case study is introduced first in \citet{fishfreezingpaper} and the latter one in \citet{ipinyoudataset}. We conduct all experiments on an M1 Pro chip with 16 GB memory. The first case study is developed with MATLAB/Simulink and the latter with Python. The code is also available \footnote{https://github.com/buseskorkmaz/ARTEO}.

\subsection{The Industrial Refrigeration System}

The industrial refrigeration system is formed by two subsystems with different pressure levels, which consists of 11 compressors such as 3 reciprocating, and 8 screw compressors to freeze stored fish within the refrigeration system. The screw compressors have varying ranges of cooling capacities and different power consumption characteristics. The first subsystem has 5 screw compressors and the second subsystem has 3 screw compressors. The system with 5 screw compressors is chosen to apply our approach. Figure \ref{Figure1} illustrates a simplified schematic of this subsystem.

\begin{figure}[tb]
  \centering
  \framebox{
  \includegraphics[scale=0.135]{Drawing.png}}
  \caption{A simplified schematic of the refrigeration system}
  \label{Figure1}
\end{figure}

\begin{figure}[t]
  \centering
  \framebox{
  \includegraphics[scale=0.45]{load_scenario.png}}
  \caption{Cooling load scenario}
  \label{Figure2}
\vspace{-1em}
\end{figure}

\begin{figure*}[t]
  \centering
  \framebox{
  \includegraphics[scale=0.3]{start-compressors.png}}
  \caption{The first iteration of the algorithm. The safe seed set includes three points for each compressor within their respective operating intervals. The blue-shaded area represents the uncertainty which is high due to unknown regions.}
  \label{Figure3}
\end{figure*}

\begin{figure*}[h]
  \centering
  \framebox{
  \includegraphics[scale=0.3]{end-compressors.png}}
  \caption{The last iteration in the simulation. GP predicted performance lines are converged to the actual performance curves of each compressor. Purple-coloured sample points are chosen by exploration when the environment becomes available. The algorithm encourages exploration specifically in the bounds of operating intervals since these points have high uncertainty as seen in Figure 3.}
  \label{Figure4}
\end{figure*}

The cooling demand is subject to change based on the stored amount of fish and environmental factors. Figure \ref{Figure2} demonstrates the production load during fishing season when the plant is needed to be cooled at full capacity for most days. This period is chosen to optimize since large production requires high energy consumption, which is a constraint in most operations due to high energy prices or due to a shortage in energy supply in case of full production.
More specifically, we demonstrate our approach on Saturday, Sunday and Monday because those dates have significant fluctuations in the cooling demand, which lead the plant to work at different operating points. Operating at unknown points may cause the violation of constraints due to unknown performance characteristics at those points.

For the dynamics of the screw compressors, 200 seconds is considered to be an average time to move to a new steady-state based on the dynamic experiments from \citet{dynamiccompressor}. Therefore, ARTEO is triggered by 200 seconds time intervals. For the performance of screw compressors, the data is obtained from \citet{perfmulcompressor}. The system consists of one small-sized, one medium-sized and three large-sized compressors. The minimum and maximum cooling capacities for each compressor type are given in Table \ref{Table1}. The coefficient of performance for all compressors is assumed to be 1.6 \citep{fishfreezingpaper}.

\begin{table}[h]
\caption{Operating Intervals of Screw Compressors}
\label{Table1}
\begin{center}
\begin{tabular}{|p{0.28\linewidth}||p{0.28\linewidth}|p{0.25\linewidth}|}
\hline
Compressor Size & Min. Capacity (kW) & Max. Cooling Capacity (kW)\\
\hline
Small & 56 & 220\\
\hline
Medium & 237 & 537\\
\hline
Large & 194 & 795\\
\hline
\end{tabular}
\end{center}
\end{table}

We proceed as next to put this problem into our framework. We first model the plant in MATLAB/Simulink with the plant's dynamic characteristics. We choose the Squared Exponential kernel, which is a positive definite kernel function, as the GP prior for each compressor based on the suggestion of the GP hyperparameter optimization tool in Matlab \citep{RasmussenGP}. The minimum and maximum cooling capacities in Table 1 are introduced as decision-set boundaries to the optimization as in Equation (8).
\begin{equation}
    x_{i\min}\leq x_{it} \leq x_{i\max} \quad \forall i,t
\end{equation}
Then, we create the safe seed set including three safe decision points for each compressor to obtain the posterior of GPs before making any decision. We define the utility function $f$ as follows:
\begin{equation}
f_t(x) = \bigl[\sum_{i=1}^{5}{\mu_{P_i(x_{it})}}\bigr]^2 + [E_t - \sum_{i=1}^{5}{x_{it}}]^2 + z\sum_{i=1}^{5}{\sigma_{P_i(x_{it})}}
\end{equation}
where $\mu_{P_i(x_{it})}$ represents the mean prediction of power consumption at production load $x$ and $\sigma_{P_i(x_{it})}$ represents the standard deviation of prediction to represent uncertainty in predictions at load $x_{it}$. $E_t$ denotes the desired production load at time $t$ based on the introduced cooling scenario. We develop the following two conditions and assign a positive value $\zeta$ = 1000 to $z$ when these conditions hold: (1) the cooling demand is the same with two immediate previous states, (2) there is at least one feasible decision point that satisfies the cooling demand with a tolerance $\alpha$ = 10 that is explored in two immediate previous states. The safety threshold for power consumption is set as 80\% of the maximum power consumption for the sum of all compressors' power consumptions. The safety function $g$ is defined next as:
\begin{equation}
g_t(x) = \sum_{i=1}^{5}{\bigl[\mu_{P_i(x_{it})} + \sqrt{\beta_t}\sigma_{P_i(x_{it})}\bigr]}
\end{equation}
where $\beta_t$ is chosen as explained in Equation (4).

\begin{figure}[t]
  \centering
  \framebox{
  \includegraphics[scale=0.35]{end-metrics-target-load.png}}
  \caption{Target and achieved cooling demands during simulation. The ARTEO algorithm follows and satisfies the expected demand and safety constraints.}
  \label{Figure5}
\end{figure}

The experiment is simulated with the given cooling load trajectory in Figure \ref{Figure2}. Figure \ref{Figure3} demonstrates the first iteration of the algorithm. In the first iteration, the uncertainty is high and it can be recognized from the wide confidence bounds. At each iteration, GPs are updated with new observations, and the algorithm takes new decisions by the updated GPs. Figure \ref{Figure4} shows the final results for performance curve estimations at the end of the simulated cooling load scenario. The ARTEO algorithm is able to successfully approximate the power consumption curves as close to their real shape.

Figure \ref{Figure5} illustrates the expected and achieved production loads. The expected production load is satisfied unless: (1) it is higher than what compressors can achieve without crossing the safety threshold, (2) it requires compressors to operate lower than their minimum operating points. Figure \ref{Figure6} demonstrates that when the total power consumption is on the limit of the maximum power threshold, expected load could not be satisfied. At those production loads, our approach finds safe operating points which minimize the stated objective function.

\begin{figure}[t]
  \centering
  \framebox{
  \includegraphics[scale=0.35]{end-metrics-power-consumption.png}}
  \caption{The total power consumption of all compressors over time. The algorithm does not violate the maximum power constraint.}
  \label{Figure6}
\end{figure}

\subsection{Online Bid Optimization}

In the second experiment, we investigate the implementation of ARTEO in a multi-dimensional problem of online bid optimization from the advertiser perspective. In bid optimization, the advertiser sets bid values with the aim of achieving high volumes by maximizing the number of shown advertisements and high profitability by maximizing the return-on-investment (ROI) ratio. In most agreements, unsatisfied ROI causes financial losses for advertisers. It becomes more challenging to sustain high ROIs when the number of advertisements increases. Constraining the ROI to remain above a certain threshold is a common approach, however, this method does not guarantee to satisfy the ROI constraint with zero violation \citep{SafeRTB}. ROI is measured by the revenues and costs, which are unknown to the bidding algorithms and brings uncertainty to the online bid optimization problem. Therefore, safe optimization algorithms could be useful to set bid values under the uncertainty of the revenues.

We apply the ARTEO algorithm to the iPinYou dataset \citep{ipinyoudataset}. This dataset has been released by a leading DSP (Demand-Side Platform) in China and consists of relevant information for personalized ads such as creative metadata, interests of users, and advertisement slot properties which include slot width, height, price and visibility. We simulate our approach by creating different campaign subsets from the original data, and for each campaign $t$, we minimize the following utility function:
\begin{equation}
\begin{split}
f_{t}{(x)} = \sum_{j=1}^{m}{cx_{tj}\mu_{C_{tj}}}  + \sum_{j=1}^{m}{\sqrt{(x_{tj} - \mu_{BP_{tj}})^2}} - \\ z\sum_{j=1}^{m}{\bigl(\sigma_{C_{tj}}+\sigma_{BP_{tj}}\bigr)}
\end{split}
\end{equation}
where $j$ denotes the ad number in the campaign, $\mu_{C_{tj}}$ is the mean prediction of GP for the $j$th advertisement to get a click with set bid values $x$, and $\mu_{BP_{tj}}$ is the mean prediction of GP for the bid price of $j$th ad in campaign $t$ based on previous similar observations. The fixed budget constraint for $m$ number of ads in campaign $t$ is formulated as:
\begin{equation}
\begin{split}
&\sum_{j=1}^{m}{x_j} \leq 180 m \\
\end{split}
\end{equation}
The safe ROI constraint is constructed for the threshold $h_t$ for the campaign $t$ as follows:
\begin{equation}
\begin{split}
&{\frac{\sum_{j=1}^{m}\mu_{C_{tj}} - \sqrt{\beta_t} \sigma_{C_{tj}}}{\sum_{j=1}^{m}{x_{tj}}}} \geq h_t \\
\end{split}
\end{equation}
Lastly, the bid values $x_{tj}$ are bounded with non-negativity for all campaigns $t$ and for all advertisements $j$.
\begin{figure}[t]
  \centering
  \framebox{
  \includegraphics[scale=0.37]{ROI-fig.png}}
  \caption{The minimum and achieved ROIs for campaigns. The ARTEO algorithm is able to remain above the threshold for all campaigns.}
  \label{Figure7}
\end{figure}

As opposed to the first experiment, where the changes in optimization goals were driven by changes in cooling demand, an increase in $t$ in the online bid optimization example is driven by the new input as starting a new campaign. We construct two GPs in this experiment, the first GP learns the bid prices from past observations, and the second one models the impressions, which are represented in binary for clicks. The impressions are traditionally predicted by classifiers due to their binary representation. However, it is possible to cast it as a regression problem where we decide the binary representations after thresholding. Since we use covariance functions of GPs to model uncertainty, we cast it as a regression and guide our RTO with continuous values. The optimization algorithm bids an ad comparatively high when its value is higher than others. Different feature sets in the dataset are used to compute posteriors based on relevance to the predictions. The GP of the bid price is initialized with the Matern kernel with $\nu$ = 1.5 and is trained over 143 features whereas the impression GP has the Squared Exponential kernel and 69 features. The safe seeds start with 30 samples, which is higher than our first experiment since this is a higher dimensional problem. As a minimum ROI threshold, 90\% of the given benchmark data ROI is set due to having a strict budget and ROI requirements in our setup. We partition the selected subset of the dataset into 25 campaigns. Each campaign has its ROI threshold and budget, which are calculated as Equation (12) and Equation (13).

The algorithm starts with a safe seed set to compute the posteriors of GPs, and then for each campaign, it utilizes the mean and standard deviation of GP posteriors to measure ROI and click probability. During the RTO phase, the higher bid prices for higher estimated click values are encouraged within a fixed budget, and the difference between the predicted bid price by GP and the proposed bid price by RTO for each advertisement is accumulated and introduced as a penalty in the objective function. Thus, the algorithm does not put the entire budget into the highest-valued ad within the campaign. At the end of each campaign, true bid prices and clicks with additive Gaussian noise are used to update the posteriors of GPs. The feedback is given only for ads in the campaign with a non-negative optimized bid price which leads to high standard deviations for non-bid similar ads. The environment becomes available for exploration after spending less than the sum of predicted bid prices and satisfying minimum thresholds in two consecutive campaigns. Hence, $z$ is set to $\zeta$ = 100 to excite the RTO to take decisions at points that could reduce uncertainty in predictions. The results of the simulation are given in Figure \ref{Figure7} and Figure \ref{Figure8}. Our approach remains above the safety threshold while proposing lower bid prices compared to the benchmark dataset.

\begin{figure}[t]
  \centering
  \framebox{
  \includegraphics[scale=0.38]{Comparison-bidding-prices.png}}
  \caption{The mean bid prices for the given benchmark and ARTEO for each campaign. ARTEO achieves higher ROI with lower costs.}
  \label{Figure8}
\end{figure}

\section{CONCLUSIONS}

In this work, we study real-time optimization of sequential-decision making in a safety-critical system under uncertainty due to unknown parameters. We propose a safe exploration and optimization algorithm which follows the optimization goals and does exploration without violating the constraints. We adaptively control the contribution of exploration via the uncertainty term in the utility function. We demonstrate the applicability of our approach to real-world datasets with industrial process optimization and online bid optimization. Our algorithm is able to model the unknown parameters in both applications and satisfy the optimization goals and safety constraints. We believe our study provides important examples to show how uncertainty in a safety-critical system can be modelled and guide the exploration.


\section{Introduction}
Artificial neural networks (ANNs) have been adopted for a wide range of classification tasks, such as image or sound classification. Both training and inference of the model mostly take place on powerful platforms that ideally contain GPUs for acceleration. However, with more neural network applications moving closer to the edge for reasons of data privacy and immediacy, machine learning faces new challenges brought on by the power constraints and limited computing resources found there~\cite{li2020survey}.

Spiking neural networks (SNNs)~\cite{neil2016learning, farabet2012comparison, o2013real} are a promising technology for deployment of AI at the edge~\cite{lu2021autonomous}. Experiments, which directly compared them with conventional ANNs, have revealed significant energy and latency advantages~\cite{lee2021accurate, blouw2019benchmarking, stromatias2015scalable}.  A reason for that is their closer resemblance to biological brains, which inspires new computing paradigms and benefits from custom hardware architectures~\cite{yang2022coreset, luo2021nc, aung2021deepfire, gerlinghoff2021e3ne, lin2020scalable, zhang2020low}. SNNs encode information in spike trains, which are series of binary values to indicate the presence (one) or absence (zero) of a spike event. Various methods exist to encode integer or floating point values into spike trains~\cite{kempter1999spike, borst1999information, gutig2014spike, tang2020rank, sboev2020solving}. In this work, we use the firing-rate encoding as detected in the brain by Zhang~\& Linden~\cite{zhang2003other}, where the number of spikes within the fixed-length spike train is proportional to the real value.

Training the weights of SNNs to yield the correct output spike train, however, poses its own difficulties. Neuron models, such as leaky integrate-and-fire~\cite{koch1998methods}, are not differentiable. This means that the classical backpropagation used for ANNs cannot be used directly for SNN training. Instead, there exist three methods for the training of SNNs:
\begin{enumerate}
    \item {\em Transfer learning.} An equivalent ANN model is trained using backpropagation while considering the constraints of the SNN. The weights are then transferred onto the SNN~\cite{rueckauer2017conversion, sengupta2019going}.
    \item {\em Approximation.} Backpropagation of errors is done using approximations or substitution of the non-differentiable aspects of the neuron model~\cite{bohte2000spikeprop, thiele2019spikegrad, shrestha2018slayer, qiao2021direct}.
    \item {\em Direct training methods.} Weights are adjusted based directly on spike events of the neurons. There exist numerous variants of the unsupervised spike-timing-dependent plasticity (STDP) and the supervised Widrow-Hoff learning rule~\cite{gerstner1996neuronal, caporale2008spike, ponulak2006supervised, mohemmed2012span}.
\end{enumerate}

High classification accuracy has been achieved with the first approach of transfer learning. These algorithms, however, are computationally expensive and may not be suitable for certain use cases, especially on the edge, where the model is required to adapt to changes in the input data through online learning. This work presents a learning algorithm that is based on STDP. It assigns desire values to all neurons, which enables supervision for weight updates to reduce the output error. The learning rule is applied to both output and hidden layers and thus allows the training of multi-layer SNN.

This paper first provides the fundamentals of SNN operation and learning in Section~\ref{sec: background}. This is followed by a description of our proposed learning rule in Section~\ref{sec: proposed learning procedure}.
The results of our experiments are presented in Section~\ref{sec: experiment}. An overview over existing spike-based learning algorithms is given in Section~\ref{sec: related works}. Lastly, we draw a conclusion.

\section{Background} \label{sec: background}
Like traditional artificial neural networks, SNNs have a layered structure with each layer consisting of an array of neurons. Adjacent neurons are connected via synapses, which can be understood as neuronal junctions to transfer spikes. The two neurons connected by a synapse are referred to as pre-synaptic and post-synaptic neuron, respectively. In this section, we consider a generic fully-connected layer. Post-synaptic neurons of this layer are indexed with $o$. They receive spikes from pre-synaptic neurons $i$.

\subsection{Neuron model} \label{sec: neuron model}
Neuron models are a mathematical description of the neuron behavior observed in biological experiments. Both the Hodgkin-Huxley model~\cite{hodgkin1952quantitative} and the Izhikevich model~\cite{izhikevich2003simple} are accurate representations, but are computationally expensive when used for large networks consisting of thousands of neurons. The Spike Response Model (SRM)~\cite{gerstner1993spikes} simulates the neuron's response to a spike using filter kernels. The SRM is a generalization of the Leaky Integrate-and-Fire (LIF) model~\cite{koch1998methods}. Both are biologically plausible, albeit simpler than the first two models. For that reason, they are commonly used for efficient SNN implementations~\cite{fang2020encoding, cassidy2013cognitive, girau2007massively}. Our experiments are also based on LIF neurons.

The information transfer between neurons comprises binary spike events, which can occur at any time step $t \in [0, T-1]$, with $T$ denoting the length of the spike trains. A LIF neuron $o$ retains its internal state, i.e. its membrane potential $p_o$, between the time steps. Additionally, a leak term $\beta_p \in [0, 1]$ causes a decay of the membrane potential over time. An input spike $s_i$ causes $p_o$ to be increased by weight $w_{oi}$, where $w_{oi}$ can be a signed value. If the membrane potential surpasses a threshold $\theta_p$, an output spike $s_o$ is fired and the membrane potential is reset. Due to the non-linearity inherent in the neuron model, no additional activation function is needed at the neuron output. The behavior of the LIF neuron can be illustrated as in Figure~\ref{fig: lif neuron}, and expressed by Equation~\ref{eqn: lif neuron}.

\begin{equation} \label{eqn: lif neuron}
    p_o[t+1] = \beta_p p_o[t] + \sum_i w_{oi} s_i[t] - \theta_p s_o[t]
\end{equation}

\subsection{Spike-timing-dependent plasticity} \label{sec: spike-timing-dependent plasticity}
According to Hebb~\cite{hebb1949organization}, a weight between two neurons is influenced by the spiking activity of both these neurons. This led to the development of the spike-timing-dependent plasticity (STDP) rule~\cite{gerstner1996neuronal, caporale2008spike}, which is a purely spike-based method of training spiking neural networks. It follows the idea of \textit{neurons that fire together wire together}~\cite[p.~64]{shatz1992developing}. Synaptic weight updates are dependent on relative differences between the spike timing of pre- and post-synaptic neurons. A pre-synaptic spike occurring before a post-synaptic spike is an indication that the pre-synaptic neuron $i$ contributes to the firing of the post-synaptic neuron $o$. In this case, the connection between $i$ and $o$ is strengthened by increasing weight $w_{oi}$. And vice versa, if pre-synaptic neuron $i$ fires before post-synaptic neuron $o$, the weight connecting them is decreased. The update of weight $w_{oi}$ at a time step $t$ can be expressed as:

\begin{equation} \label{eqn: stdp}
    \Delta w_{oi}[t] = \eta^+ s_o[t] r_i[t] - \eta^- s_i[t] r_o[t].
\end{equation}

\noindent
where $s_i$ and $s_o$ are binary spike trains of input and output neuron, respectively. $\eta^+$ and $\eta^-$ are learning rates. Decaying spike traces $r_i$ and $r_o$ account for the time delay between pre- and post-synaptic spikes. The magnitude of the weight update increases as the time period between pre- and post-synaptic spikes becomes shorter. Spike traces are generated by first decaying any existing trace value with an exponential kernel. After that, the trace is incremented in case the neuron fires a spike (see Equation~\ref{eqn: traces}). The decay rate $\beta_r < 1$ is a hyper-parameter, which can be trimmed to adjust how the delay between spikes influences the weight update. Figure~\ref{fig: stdp} shows the dynamics of spike traces and weight update expressed by Equation~\ref{eqn: stdp}.

\begin{equation} \label{eqn: traces}
    r_{\{o|i\}}[t] = \beta_r * r_{\{o|i\}}[t-1]+ s_{\{o|i\}}[t]
\end{equation}

\begin{figure}[t]
    \centering
    \begin{minipage}[t]{0.48\textwidth}
        \centering
        \includegraphics[width=0.84\textwidth]{figures/neuron_model.png}
        \caption{Input spikes $s_i$ and membrane potential $p_o$ of LIF neuron $o$ for every time step $t \in [0, 6]$. If threshold $\theta_p$ is exceeded, output spikes $s_o$ are generated.}
        \label{fig: lif neuron}
    \end{minipage}
    \hfill
    \begin{minipage}[t]{0.48\textwidth}
        \centering
        \includegraphics[width=0.8\textwidth]{figures/spike_traces.png}
        \caption{Input spikes $s_i$ and output spikes $s_o$ with respective traces $r$ for six time steps $t$.}
        \label{fig: stdp}
    \end{minipage}
\end{figure}

\section{Proposed learning rule} \label{sec: proposed learning procedure}
Spike-timing-dependent plasticity is an effective way of adjusting the relationship between neurons in a way that increases their sensitivity to important features present in the dataset~\cite{fu2021ensemble}. The weight updates are, however, required to be guided, if one wishes to perform supervised learning. In that case, besides extracting features, neuron weights are tuned for reduction of the network's output error. Each neuron in the network is assigned a {\em desire} value, which is derived from the training labels and determines the direction of the weight update. Through local loss computation, desire values can be assigned to the neurons in the hidden layers.

Since we adopt the backpropagation algorithm, the training routine follows that of traditional deep neural networks. During the forward pass, samples of the training dataset, which are encoded as spike trains, generate spikes and spike traces in the hidden and output layers. The backward pass employs the desire backpropagation algorithm to derive desire values from training labels and output error. Lastly, the weights are updated following the STDP rule, utilizing spike trains/traces and being guided by desire values.

\subsection{Desire-based weight update} \label{sec: weight update}
The original STDP rule updates the weight $w_{oi}$ whenever there is a spike in either neuron $i$ or neuron $o$. For supervised learning, however, weight updates should be tied to the activity of the post-synaptic neuron~$o$ only, as those are used to compute the subsequent layer and ultimately lead to the output result. Our weight update therefore only takes the post-synaptic spike train $s_o$ and pre-synaptic trace $r_i$ into account. Weights between neurons with higher activity experience larger adjustments, following the idea of Hebbian learning. This is expressed in Equation~\ref{eqn: desire weight update}, which applies equally to all weights in the network.

\begin{equation} \label{eqn: desire weight update}
    \Delta w_{oi}[t] = d_o * \eta s_o[t] r_i[t].
\end{equation}

The STDP term of the update equation is not sufficient for supervised learning. That is, instead, accomplished by the desire term $d_o$, which serves two purposes. Firstly, it decides whether the weight is increased or decreased. Secondly, the weight update is masked if the neuron's spikes have an insignificant influence on the successive layer. Therefore, we define a ternary value $d_o \in \{-1, 0, +1\}$ with the three states describing that:
\begin{itemize}
    \item $+1$: the neuron $o$ ought to spike and its weights $w_{oi}$ should increase for all $i$,
    \item $-1$: the neuron $o$ should not spike and its weights should decrease, or
    \item $0$: the activity of neuron $o$ is indifferent and the weights should not change
\end{itemize}

Besides disabling the weight update, a desire of $0$ will cause the neuron to not be considered during backpropagation, as will be seen in the next section. The proposed learning rule incorporates STDP, while the desire value depends on the output error. Therefore, it is considered a spike-based supervised learning rule.

\subsection{Desire backpropagation} \label{sec: desire backpropagation}
So far, we have used indices $i$ and $o$ to index pre- and post-synaptic neurons of a generic linear layer. In the following, however, there is a differentiation between the output layer and the hidden layers of the network. Figure~\ref{fig: network architecture} sketches a network with two subsequent hidden layers $G$ and $H$, and the last two layers $J$ and $K$, each with a minimal number of neurons for the sake of readability. We replace indices $o$ and $i$ according to the respective layer descriptors. The figure gives an overview of the steps involved in the computation of the desire values $d_k$ and $d_g$, based on the output error and local errors, respectively.

\begin{figure}[t]
    \centerline{\includegraphics[width=0.75\textwidth]{figures/network_architecture.png}}
    \caption{Minimal neural network, which only shows two arbitrary hidden layers $G$ \& $H$ and the last two layers $J$ \& $K$. The steps of computing output desire $d_k$ and hidden desire $d_g$ are visualized by the green and blue arrows, respectively.}
    \label{fig: network architecture}
\end{figure}

\subsubsection{Output layer}
Supervised learning algorithms aim to reduce the output error of the last network layer. In rate-encoded SNNs, each of the network's output neurons $k$ generates a spike sequence $s_k$ with $T$ binary elements. Accumulation over that sequence yields the spike count. The output error $e_k$ of a neuron is the difference between the normalized spike count and a target $\hat{s_k}$, given in Equation~\ref{eqn: desire error output}. The target is derived from the training labels. If neuron $k$ represents the target class, it should fire $T$ spikes, and remain silent otherwise.

\begin{equation} \label{eqn: desire error output}
    e_k = \frac{1}{T} \sum_{t=0}^{T-1} s_k[t] - \hat{s_k}
    \quad \text{with} \quad
    \hat{s_k} =
    \begin{cases}
        1, & \text{if } k \text{ is target class} \\
        0, & \text{otherwise}
    \end{cases}
\end{equation}

To quantify the performance of the network we apply the squared error loss function in Equation~\ref{eqn: desire loss} to the output error $e_k$. The derivative of the loss function with respect to a single output spike train $s_k$ is required for backpropagation (see Equation~\ref{eqn: desire loss derivative}).

\begin{equation} \label{eqn: desire loss}
    L = \frac{1}{2} \sum_k e_k^2
\end{equation}

\begin{equation} \label{eqn: desire loss derivative}
    \frac{\partial L}{\partial s_k} = e_k
\end{equation}

The update of weight $w_{kj}$ is guided by a desire value $d_k$ according to Equation~\ref{eqn: desire weight update}. Due to the structural similarity between feed-forward ANN and SNN, we derive the desire by applying gradient decent. The gradient of the loss function with respect to the weight is given in Equation~\ref{eqn: desire backpropagation 1}, where variables are denoted using the SNN symbols from the previous sections. Due to the non-differential LIF neuron model, we follow~\cite{shrestha2019approximating} and approximate the second term with a straight-through estimator, which has a gradient of one. The first term describes the impact of the weight on the membrane potential $p_k$. The binary nature of the spikes causes the gradient to be 1, resulting in Equation~\ref{eqn: desire backpropagation 2}.

\begin{equation} \label{eqn: desire backpropagation 1}
    \frac{\partial L}{\partial w_{kj}} =
    \frac{\partial p_k}{\partial w_{kj}}
    \frac{\partial s_k}{\partial p_k}
    \frac{\partial L}{\partial s_k}
\end{equation}

\begin{equation} \label{eqn: desire backpropagation 2}
    \frac{\partial L}{\partial w_{kj}} =
    \frac{\partial L}{\partial s_k}
\end{equation}

By following the descending gradient of the loss function, we obtain an expression for the desired change in spiking activity, which leads us to our desire value $d_k$. Since the magnitude of the weight update is determined by the spike timings using STDP, we apply the ternarization function $\mathscr{T}_\theta$ (see Equation~\ref{eqn: desire ternarization}). $\sgn(\cdot)$ is the sign function, which returns $-1$ for negative and $+1$ for positive values. $\theta$ is a threshold value to stabilize the training for small errors, which we call $\theta_{\text{out}}$ for the output layer in particular. That yields the final expression in Equation~\ref{eqn: desire value output}. In case of the output layer, only weights connecting to the target neuron can get increased, while all others might get decreased.

\begin{equation} \label{eqn: desire ternarization}
    \mathscr{T}_\theta(x) =
    \begin{cases}
        \sgn x, & \text{if } |x| > \theta \\
        0,      & \text{otherwise}
    \end{cases}
\end{equation}

\begin{equation} \label{eqn: desire value output}
    d_k =
    \mathscr{T}_{\theta_{\text{out}}} \left(-\frac{\partial L}{\partial w_{kj}} \right) =
    \mathscr{T}_{\theta_{\text{out}}} \left(-\frac{\partial L}{\partial s_k} \right) =
    \begin{cases}
        -\sgn e_k, & \text{if } |e_k| > \theta_{\text{out}} \\
        0, & \text{otherwise}
    \end{cases}
\end{equation}

\subsubsection{Hidden layers}
Desire backpropagation enables the training of multi-layer neural networks. In contrast to previous methods, such as~\cite{tavanaei2019bp}, where the output error of the last layer is backpropagated, we compute the error locally. That allows us to train networks with a higher number of layers. As shown in Figure~\ref{fig: network architecture}, desire values $d_g$ of hidden layer $G$ are determined based on the desire values and spike trains of layer $H$.

In the previous section, we derived the dependency of a desire value on the loss gradient with regards to the spike train. To determine this gradient for a hidden neuron, we use the classical backpropagation as a starting point (see Equation~\ref{eqn: desire backpropagation 3}). Just as before, the straight through estimator and the squared error loss is applied to obtain Equation~\ref{eqn: desire backpropagation 4}.

\begin{equation} \label{eqn: desire backpropagation 3}
    \frac{\partial L}{\partial s_{g}} = \sum_h
    \frac{\partial p_h}{\partial s_g}
    \frac{\partial s_h}{\partial p_h}
    \frac{\partial L}{\partial s_h}
\end{equation}

\begin{equation} \label{eqn: desire backpropagation 4}
    \frac{\partial L}{\partial s_{g}} = \sum_h
    w_{hg} * e_h
\end{equation}

To calculate the loss locally for layer $H$, the error is only dependent on the spike output and the desire values of this layer. Equation~\ref{eqn: desire error hidden} transforms the ternary desire values $d_h$ into targets of $\{0, 1\}$ for the normalized spike count. If the desire value is zero, however, an error cannot be calculated due to the indifference in spiking desire of the neuron. The loss is obtained by squaring the error as done before in Equation~\ref{eqn: desire loss}.

\begin{equation} \label{eqn: desire error hidden}
    e_h =
    \begin{cases}
        \frac{1}{T} \sum_{t=0}^{T-1} s_h[t] - (d_h + 1) / 2, & \text{if } d_h \neq 0 \\
        0, & \text{otherwise}
    \end{cases}
\end{equation}

The desire values of the preceding hidden layer $G$ are obtained by supplying the negative gradient to the ternarization function $\mathscr{T}_\theta$ with a threshold $\theta_{\text{hid}}$ specific for hidden layers. The final expression of a desire value $d_g$ is given in Equation~\ref{eqn: desire value hidden}. In summary, the computation only relies on desire values, spike trains and weights of layer $H$.

\begin{equation} \label{eqn: desire value hidden}
    d_g =
    \mathscr{T}_{\theta_{\text{hid}}} \left(-\frac{\partial L}{\partial s_g} \right) =
    \begin{cases}
        -\sgn \left( \sum_h w_{hg} e_h \right), & \text{if } \left| \sum_h w_{hg} e_h \right| > \theta_{\text{hid}} \\
        0, & \text{otherwise}
    \end{cases}
\end{equation}

Using Equation~\ref{eqn: desire value hidden}, the desire values of layer $G$ are set such that the desire in layer $H$ is satisfied. Although a pre-synaptic neuron $g$ cannot fulfill the desire of all post-synaptic neurons $d_h$, the sum over all $w_{hg} e_h$ gives an indication about the preference in the spiking activity of $g$. If neuron $g$ does not significantly contribute to the average fulfillment of the desire in layer $H$, no adjustment is made to its weights. This is reflected by a desire $d_g$ of zero.

If the network has been trained and has reached an equilibrium, the local errors are small enough to draw desire values to zero. In that case, weights theoretically stabilize and are not adjusted. The network's ability to generalize depends, however, on the number of neurons and layers. In practise, the output error will not be zero for all samples, which leads to continuing weight adjustments even towards the end of the training.

\subsection{Dropout} \label{sec: dropout}
After many training epochs, the accuracy of classifying the training dataset usually increases, while the accuracy on the validation dataset might remain unchanged or even decrease. This is the phenomenon of overfitting.  Dropout~\cite{srivastava2014dropout} is a technique to increase the generality of a network model, and has been shown to be effective for SNNs~\cite{lee2020enabling}.

Dropout simulates the training of multiple neural networks by randomly excluding neurons of each layer during the training phase. During the evaluation phase, dropout is not applied, which is equivalent to combining the results of all those neural networks. The parameter $P_{\text{drop}}$ gives the probability of a neuron to be dropped out. Since the neuron configuration needs to remain constant for all time steps, a binary dropout mask $m$ is generated for each layer at the beginning of each sample. An input neuron $i$, which is silenced by $m_i = 0$, does not contribute to the desire backpropagation and does not get its weights updated.

Dropout leads to a reduction of the number of spikes arriving at the next layer. To compensate for that, the influence of the active neurons on the membrane potential is increased by a factor $1 / (1 - P_{\text{drop}})$. Equation \ref{eqn: lif neuron} for the training phase is modified to include dropout as follows:

\begin{equation} \label{eqn: dropout}
    p_o[t+1] = \beta_p p_o[t] + w_{oi} s_i[t] * \frac{m_i} {1 - P_{\text{drop}}} - s_o[t].
\end{equation}

\section{Experiments and results} \label{sec: experiment}
\subsection{Spike generation}
For the sake of comparability, we used popular image datasets for classification. Since those datasets are only available as two-dimensional pixel data, spike trains with additional timing information need to be generated. For rate encoding, the number of spikes in the input spike train is proportional to the brightness value of the respective pixel. We do not carry out any data preprocessing or augmentation.

Our input layer contains LIF neurons (see Section~\ref{sec: neuron model}), which increment their membrane potential at each time step by a constant, which is equal to the the brightness value of the input $\in [0, 1]$. The membrane potential is not decayed, i.e. $\beta_p = 1$. A brighter pixel leads to larger increments of the membrane potential $p_o$ at every time step. Hence, the threshold $\theta_p$ is reached faster and more often and more spikes are generated. The dropout rate $P_{drop}$ can be adjusted separately for the input layer.

\subsection{MNIST classification} \label{sec: MNIST classification}
\begin{figure}[t]
    \centering
    \begin{minipage}[t]{0.49\textwidth}
        \centering
        \includegraphics[width=\textwidth]{figures/plot_architecture.png}
        \caption{Classification accuracy (in \%) plotted against the number of hidden neurons $N$ in a 3-layer MLP architected as 784 -- $2N$ -- $N$ -- 10 and $L=20$.}
        \label{fig: plot architecture}
    \end{minipage}
    \hfill
    \begin{minipage}[t]{0.49\textwidth}
        \centering
        \includegraphics[width=\textwidth]{figures/plot_timesteps.png}
        \caption{Accuracy (in \%) versus spike train length / number of time steps $T$ with $N$ set to 800.}
        \label{fig: plot timesteps}
    \end{minipage}
\end{figure}

\begin{figure}[t]
    \centerline{\includegraphics[width=0.50\textwidth]{figures/plot_local_errors.png}}
    \caption{Average local test losses for each neuron in layer 1, 2 and 3 throughout the training with 150 epochs. Accuracy on test dataset for each epoch (gray line).}
    \label{fig: local errors}
\end{figure}

\begin{figure}[t]
    \centerline{\includegraphics[width=1.0\textwidth]{figures/plot_spike_activity.png}}
    \caption{Color-coded spike activity of neurons in the three layers over five training epochs. Plots are separated into columns for each of the ten handwritten digits 0 to 9. Four regions of interest are marked by green circles.}
    \label{fig: spike activity}
\end{figure}

\begin{figure}[t]
    \centerline{\includegraphics[width=1.0\textwidth]{figures/plot_contribution.png}}
    \caption{Contribution of each of the input layer neurons to satisfy on average the desire of the neurons in the next layer, measured for the first 1000 samples in the first training epoch. White color represents no contribution, while the shades of red and blue stand for negative and positive contribution, respectively.}
    \label{fig: contribution}
\end{figure}

We trained a multilayer perceptron (MLP) on 60,000 training images of the handwritten digits dataset MNIST~\cite{lecun1998mnist}. The accuracy values in this section were obtained from the 10,000 test images. The network consists of three fully-connected layers with a variable neuron count of 784 -- $2 N$ -- $N$ -- 10. We use a variable number of time steps $T$ for the spike trains. Our PyTorch implementation of the spiking neural network allows an exponential decay of the learning rate. We start off with a learning rate $\eta = 1*10^{-5}$, which is reduced by 4\% at the end of every training epoch. Training is carried out for a total of 150 epochs. Other hyper-parameters include the firing threshold $\theta_p = 1.0$, the desire thresholds $\theta_{\text{hid}} = 0.05$ and $\theta_{\text{out}} = 0.30$, and the dropout rate $P_{\text{drop}} = 30\%$ for hidden layers and no dropout at the input layer.

A classification is deemed correct, if the output neuron, which corresponds to the target class, fires the strict maximum number of spikes. That is, all other output neurons fire fewer spikes. We plotted the test accuracy against the number of neurons $N$ and time steps $T$ as seen in Figures \ref{fig: plot architecture} and \ref{fig: plot timesteps}, respectively. A higher number of neurons increases the ability of the network to recognize details in the input samples. As expected, the classification accuracy of our network tends to increase with the neuron count. Saturation occurs when additional neurons do not lead to a considerable increase in accuracy. The most efficient network structure contains $2*800 + 800 = 2400$ hidden neurons, with $N=800$. It reaches a test accuracy of 98.41\%. The highest accuracy of 98.48\% can be reached at the cost of increase in neurons. With regards to Figure~\ref{fig: plot timesteps}, the hyperparameter tuning was performed for $T = 20$, explaining the peak accuracy for this configuration. With spike-based learning, every time step influences the update of weights. Although longer spike trains can represent more information, it also requires additional adjustments of learning rate and thresholds. Using only five time steps causes significant information loss in the activations, which is reflected in the drop of classification accuracy.

For a training with parameters $N = 800$ and $T = 20$, we analyze the learning rule's impact on the neurons of the network. After each epoch, we evaluate the network with the test dataset and record the local loss of each neuron in the three layers according to Equation~\ref{eqn: desire loss}. The result is plotted in Figure~\ref{fig: local errors}, which shows a clear decrease in local losses throughout the training process. Concurrently, the classification accuracy tends to increase. A deeper insight in the spiking activity of individual neurons for the first five training epochs is given in Figure~\ref{fig: spike activity}, with special regions of interest marked by green circles. We randomly selected 20 neurons from the two hidden layers and, together with the ten output neurons, plot their color-coded relative spike count. The activity is measured separately for each of the ten classes of digits 0 to 9. The last row of plots represents the output layer and, as expected, the spike activity is concentrated at the correct output neuron for each class, i.e. neuron 0 spikes for class 0, etc. One can observe that the activity of target neurons is initially lower and increases throughout the training (region~1). Accordingly, non-target output neurons tends transition towards lower activity (region~2). In the hidden layers, represented by the upper two rows of plots, neurons are shown to be sensitive to certain input classes. For example, neuron 11 in the first layer spikes for class 6, but not for any other class (region~3). This sensitivity is established mainly within the first training epoch. Bright regions at the left edges of some plots show that, initially, many neurons were active and later silenced as result of the training (region~4). Lastly, in Figure~\ref{fig: contribution}, colors encode how a neuron contributes on average to satisfy the desire of the neurons in the next layer. Contribution is a function of layer weights, input/output spikes, and desire values. Thereby, red represents a bad influence, meaning that spike activity in this neuron will cause neurons in the next layer to fire in an undesired manner. Blue color reflects positive contribution to the overall desire satisfaction. The contribution was measured for each of the 784 neurons of the first layer over the first 1000 training samples of the first epoch. Because neurons at the edges of the input image remain silent for all samples, their contribution is zero throughout. The others display negative contribution at the start of the training, gradually improving as the network is exposed to more training samples. The red patch at around 500 samples indicate a temporary drop in contribution, likely caused by out-of-distribution training samples.

Table~\ref{tbl: results mnist} compares our desire backpropagation with state-of-the-art works for SNN training. For a fair comparison, we include learning rules, which perform supervised spike-timing-driven weight updates, and employ fully-connected layers. Our method outperforms most of those learning algorithms in terms of accuracy, such as~\cite{shrestha2017stable, neftci2017event, kheradpisheh2018stdp, tavanaei2019bp, shrestha2019approximating, hao2020biologically, comsa2020temporal, liu2021sstdp, mirsadeghi2021stidi, tang2021biograd}. Those achieving a higher accuracy~\cite{zhang2018plasticity, falez2019multi} require a significant increase in the number of neurons and/or time steps. CNN architectures~\cite{kheradpisheh2018stdp, falez2019multi} especially require considerably more computational effort, with only a small impact on accuracy. While~\cite{neftci2017event, shrestha2019approximating, tang2021biograd} show advantages in terms of neuron count, they require feedback weights in addition to the feed-forward weights, which might not be suitable for memory-constraint devices. While all listed learning rules are able to train at least the last network layer in a supervised, some learning rules cannot be generalized for deeper networks~\cite{tavanaei2019bp, falez2019multi, hao2020biologically}. The table shows that desire backpropagation is a competitive learning rule, which achieves high accuracy while keeping the number of neurons low.

\begin{table}[ht]
    \caption{Comparison of state-of-the-art learning rules for the MNIST dataset with regards to accuracy (in \%), number of neurons and time steps.}
    \normalsize
    \renewcommand{\arraystretch}{1.1}
    \begin{center}
    \begin{tabular}{l|l l l|r r r}
        \hline
        \textbf{Model} & \textbf{Layers} & \textbf{Learning Rule} & \textbf{Acc.} & \textbf{Neurons} & \textbf{T.steps} \\
        \hline
        Shrestha (2017) \cite{shrestha2017stable}        & 2 FC & Stable STDP      & 89.7 &  1600 &  200 \\
        Neftci (2017) \cite{neftci2017event}             & 3 FC & eRBP             & 98.0 &  1000 & ~125 \\
        Kheradpisheh (2018) \cite{kheradpisheh2018stdp}  & CNN  & STDP+SVM         & 98.4 & 12220 &   30 \\
        Zhang (2018) \cite{zhang2018plasticity}          & 2 FC & Equ. Learn.+STDP & 98.5 &  4500 &  100 \\
        Tavanaei (2019) \cite{tavanaei2019bp}            & 3 FC & BP-STDP          & 97.2 &   650 &  ~10 \\
        Falez (2019) \cite{falez2019multi}               & CNN  & STDP+SVM         & 98.6 & 35328 &   -- \\
        Shrestha (2019) \cite{shrestha2019approximating} & 3 FC & EMSTDP           & 97.3 &  1000 &  200 \\
        Hao (2020) \cite{hao2020biologically}            & 2 FC & DA-STDP          & 96.7 & 10000 &  700 \\
        Comsa (2020) \cite{comsa2020temporal}            & 2 FC & Error BP         & 98.0 &   340 &   10 \\
        Mirsadeghi (2021) \cite{mirsadeghi2021stidi}     & 2 FC & STiDi-BP         & 97.4 &   500 & ~300 \\
        Liu (2021) \cite{liu2021sstdp}                   & 2 FC & SSTDP            & 98.1 &   300 &   16 \\
        Tang (2021) \cite{tang2021biograd}               & 3 FC & BioGrad          & 98.1 &   600 &   20 \\
        \textbf{This work}                               & 3 FC & STDP+Desire BP   & 98.4 &  2400 &   20 \\
        \hline
    \end{tabular}
    \label{tbl: results mnist}
    \end{center}
\end{table}

\subsection{Fashion-MNIST classification}
Just like MNIST, Fashion-MNIST~\cite{xiao2017fashion} consists of 60,000 training and 10,000 validation images that are 28x28 pixels. They belong to ten kinds of clothing, such as shirts, trousers, and shoes. The photographic nature of the images contains more detail than digits of the MNIST dataset and are hence more challenging to classify correctly. However, for the same reason, Fashion-MNIST is considered a more realistic dataset for computer vision applications~\cite{hao2020biologically}. We train a network with three fully-connected layers and 784 -- 1000 -- 100 -- 10 neurons. Spike trains have a length of $T = 20$. A total of 600 epochs had to be computed due to higher dropout rates of 40\% and 5\% for hidden and input layers, respectively. All other hyperparameters are equal to the ones used for MNIST classification in Section \ref{sec: MNIST classification}.

With those settings, we achieved a classification accuracy of 87.56\%. The comparison with Shrestha~et~al. \cite{shrestha2019approximating} and Hao et al.~\cite{hao2020biologically} reflects the result of the MNIST classification (see Table \ref{tbl: results fashion mnist}). Desire backpropagation achieved the highest accuracy, while using either a similar number of neurons or less.

\begin{table}[ht]
    \caption{Results and comparison of Fashion-MNIST classification with indication of accuracy (in \%), number of neurons and time steps.}
    \normalsize
    \renewcommand{\arraystretch}{1.1}
    \begin{center}
    \begin{tabular}{l|l l|r r r}
        \hline
        \textbf{Model} & \textbf{Layers} & \textbf{Learning Rule} & \textbf{Acc.} & \textbf{Neurons} & \textbf{T.steps} \\
        \hline
        Shrestha (2019) \cite{shrestha2019approximating} & 3 FC & EMSTDP         & 86.1 & 1000 & 200 \\
        Hao (2020) \cite{hao2020biologically}            & 2 FC & DA-STDP        & 85.3 & 6400 &  -- \\
        \textbf{This work}                               & 3 FC & STDP+Desire BP & 87.6 & 1100 &  20 \\
        \hline
    \end{tabular}
    \label{tbl: results fashion mnist}
    \end{center}
\end{table}

\subsection{Computational complexity}
\begin{figure}[t]
    \centering
    \begin{minipage}{0.45\textwidth}
        \centering
        \captionof{table}{Comparison of operations needed for three sub-processes using desire backpropagation and classical backpropagation. $I$/$O$ are the numbers of inputs/outputs to a neuron, $T$ is the number of time steps.}
        \normalsize
        \renewcommand{\arraystretch}{1.1}
        \setlength{\tabcolsep}{4pt}
        \begin{tabular}{l|l l l|l l}
            \hline
            \textbf{Process} & \multicolumn{3}{c|}{\textbf{SNN: Desire BP}} & \multicolumn{2}{c}{\textbf{ANN: BP}} \\
            & Mult. & Add. & Cond. & Mult. & Add. \\
            \hline
            Forward  & $1$ & $T I$ & $T I$ & $I$  & $I$ \\
            Backward & $O$ & $2O$  & $1$   & $2O$ & $O$ \\
            Update   & $1$ & --    & $1$   & $2$  & --  \\
            \hline
        \end{tabular}
        \label{tbl: complexity}
    \end{minipage}
    \hfill
    \begin{minipage}{0.52\textwidth}
        \centering
        \includegraphics[width=\textwidth]{figures/plot_complexity.png}
        \caption{Normalized execution time of forward pass, backward pass, and weight update, and memory usage during training. Numbers on plot lines are the scaling factors to indicate relationship between execution times.}
        \label{fig: complexity}
    \end{minipage}
\end{figure}

Spiking neural networks have been touted as an efficient alternative to conventional artificial neural networks. In this section, we take a deeper look into the computational complexity of desire backpropagation. For forward pass, backward pass and weight update of a single neuron, the number of operations is determined and compared with classical backpropagation in Table~\ref{tbl: complexity}. It shows the dependency of the SNN forward pass on the number of time steps $T$ and the number of inputs $I$. Binary spikes allow the use of conditional statements instead of multiplications. Only one multiplication is needed for the decay of the membrane potential, whereas ANNs require multiplications and additions for all $I$ inputs. In the backward pass, Desire BP replaces one multiplication with an addition due to the local computation of losses. One conditional statement is used in the ternarization function. The ternary nature of the desire values benefits the weight update, as one multiplication operation can be saved. Multiplications are computationally more expensive than additions, which are in turn more expensive than conditional statements. Spike-based processing in combination with our learning rule can therefore be expected to be more computationally efficient overall.

From table~\ref{tbl: complexity} we can derive linear time complexity with respect to inputs $I$ and outputs $O$. To confirm this hypothesis, an SNN for the classification of MNIST digits was constructed, consisting of input, hidden, and output layer. The number of neurons in the hidden layer is varied in logarithmic steps from 100 to 30000. We then executed the training on an NVIDIA GeForce RTX 2080 Ti GPU and measured the execution time of forward pass, backward pass and weight update. In addition, we recorded the memory utilization of the process running the Python script.

Before displaying the results in Figure~\ref{fig: complexity}, the data were individually normalized in terms of offset and range, allowing a direct comparison of the dependency on the hidden neuron number. The plot uses logarithmic scale for both horizontal and vertical axes. As expected, a roughly linear trend can be detected for time complexity. Because each neuron in the network comprises an equal number of internal variables, the memory shows linear behaviour as well. Due to parallel execution within the GPU, execution time does not always increase with larger neuron count. Instead, it maintains a constant level before increasing by a bigger step. The absolute relationship between the execution time for different processes is indicated by the scaling factors on the plot lines. The forward pass takes around $20\times$ more time to run. With the number of time steps being 20, the forward pass has to execute 20 times, while backward pass and weight update is only performed once.

\section{Related research} \label{sec: related works}
The learning rules of spiking neural networks are inspired by the inner workings of biological brains. In order for a neuron model and learning rule to be deemed biologically plausible, training and inference should be spike-based, i.e., operate on the temporal information of the spike train and interface using spike information~\cite{hao2020biologically}. Spike-timing-dependent plasticity is such a rule, that only relies on locally available spike information, and was used for unsupervised learning~\cite{masquelier2007unsupervised}. The same team of researchers enhanced this method by learning neurons in competition through a winner-takes-all mechanism that inhibits the spiking of laterally adjacent neurons to prevent them from learning the same patterns~\cite{masquelier2010learning}. Kheradpisheh et al.~\cite{kheradpisheh2018stdp} coupled STDP-based network layers with a linear SVM classifier to enable supervised learning. Similarly, Thiele et al.~\cite{thiele2018timescale} uses a custom event-based classifier to allow the learning to depend on the output error.

To use the STDP rule directly for supervised learning, it was separated into the STDP and anti-STDP processes, which refer to weight strengthening and weakening, respectively~\cite{ponulak2006supervised, sporea2013supervised, wang2014online}. In his doctoral dissertation, Ponulak~\cite{ponulak2006supervised} showed how to train a network for a timing-accurate spike sequence. The STDP process is applied whenever a spike in the desired sequence occurs, while anti-STDP is applied for every actual output spike. Once the actual spike train matches the desired one, STDP and anti-STDP process counteract each other and weights are not modified. A \textit{learning window} is used for the spike timing to scale the magnitude of the weight update. Ponulak's learning rule was adopted by Wang et al.~\cite{wang2014online}. Weights of their output layer are updated at the time of actual output spikes only. The sign of the weight update is determined by whether or not the output neuron is desired to spike. The idea of spiking desire was used earlier by the Tempotron~\cite{gutig2006tempotron} to determine the sign of the weight update, while the magnitude is proportional to the output error. In a similar manner, PBSNLR~\cite{xu2013new} tries to match the output spike timings with a desired spike train by adjusting weights, whenever there is a spike at an undesired time or a failure to spike at a desired time step. Those rules are only applied to a single-layer MLP or the output layer of the SNN as hidden neurons lack the information about spiking desire.

Other spike-based supervised learning algorithms are based on the Widrow-Hoff rule~\cite{mohemmed2012span, mohemmed2013training}. It adjusts the neurons of the output layer according to the contribution to their respective output errors. This rule is, again, applied to the last layer due to the ability to directly compute the output error. Tavanaei \& Maida~\cite{tavanaei2019bp} improved this rule to train the last two layers of a network.

Seeing the benefits of spike-based learning, recent works have adopted traditional backpropagation to SNN neuron models. Shrestha et al.~\cite{shrestha2019approximating} uses error neurons to record local errors during the forward pass. They exhibit the same neuron behaviour as forward neurons. Mirsadeghi et al.~\cite{mirsadeghi2021stidi} approximate error backpropagation in SNN and execute the weight adjustment based on STDP principles. Liu at al.~\cite{liu2021sstdp} approximate partial derivatives of the backpropagation with an STDP-like expression, enabling the error to be propagated back through the network.

Learning of spiking neural networks has also been adopted from conventional artificial neural networks. However, learning rules do often not adhere to the characteristic of biological plausibility mentioned earlier. Wu et al.~\cite{wu2019deep} apply error backpropagation to SNNs by updating the synaptic weights based on the total number of output spikes rather than spike timing. Similarly, GLSNN~\cite{zhao2020glsnn} does not consider the temporal relation between spikes when updating their weight through feedback alignment.

\section{Conclusion}
In this paper, we propose desire backpropagation, a novel lightweight learning algorithm for spiking neural networks. With STDP as its foundation, our learning rule utilizes the benefits of spike-based based learning such as energy efficiency and biological plausibility. To enable multi-layer learning, we introduced the ternary spiking desire of a neuron. It determines, whether the weight needs to be increased or decreased during the STDP update to minimize the output error. The desire value of the output layer can directly be derived from the labels of the training samples. For hidden layers, we utilized local losses, which can be calculated based on already propagated desire values. The backpropagation sets desire values in a way that helps satisfying the desire of the successive layer. The ternary desire values are coupled with a binary dropout mask to avoid overfitting of the model.

We tested our learning rule on the MNIST and Fashion-MNIST datasets and achieved remarkable classification accuracy of 98.41\% and 87.56\%, respectively. It not only shows superior accuracy compared to many other learning rules, but also uses less neurons and time steps than most other algorithms. We further compared the required operations with classical backpropagation and demonstrated its advantages in terms of computational complexity. Its performance and efficiency makes desire backpropagation a potential candidate for deployment and online learning on edge devices.

\section{Acknowledgement}
This work was supported by the Singapore Government’s Research, Innovation and Enterprise 2020 Plan (Advanced Manufacturing and Engineering domain) under Grant A1687b0033.