text:
\section{Introduction}\label{sec:introduction}

In recent years, Spiking Neural Networks (SNNs) have shown promise towards enabling low-power machine intelligence with event-driven neuromorphic hardware. Founded on bio-plausibility, the neurons in an SNN compute and communicate information through discrete binary events (or 'spikes') a significant shift from the standard artificial neural networks (ANNs), which process data in a real-valued (or analog) manner. The binary all-or-nothing spike-based communication combined with sparse temporal processing precisely make SNNs a low-power alternative to conventional ANNs. With all its appeal for power efficiency, training SNNs still remains a challenge. The discontinuous and non-differentiable nature of a spiking neuron (generally, modeled as leaky-integrate-and-fire (LIF), or integrate-and-fire (IF)) poses difficulty to conduct gradient descent based backpropagation. Practically, SNNs still lag behind ANNs, in terms of performance or accuracy, in traditional learning tasks. Consequently, there has been several works over the past few years that propose different learning algorithms or learning rules for implementing deep convolutional SNNs for complex visual recognition tasks \citep{wu2019direct, hunsberger2015spiking, cao2015spiking}. Of all the techniques, conversion from ANN-to-SNN \citep{diehl2016conversion, diehl2015fast, sengupta2019going, hunsberger2015spiking} has yielded state-of-the-art accuracies matching deep ANN performance for Imagenet dataset on complex architectures (such as, VGG~\citep{simonyan2014very} and ResNet~\citep{he2016deep} ). In conversion, we train an ANN with ReLU neurons using gradient descent and then convert the ANN to an SNN with IF neurons by using suitable threshold balancing \citep{sengupta2019going}. But, SNNs obtained through conversion incur large latency of $2000{-}2500$ time steps (measured as total number of time steps required to process a given input image\footnote{SNNs process Poisson rate-coded input spike trains, wherein, each pixel in an image is converted to a Poisson-distribution based spike train with the spiking frequency proportional to the pixel value}). The term `time step' defines an unit of time required to process a single input spike across all layers and represents the network latency. The large latency translates to higher energy consumption during inference, thereby, diminishing the efficiency improvements of SNNs over ANNs. To reduce the latency, spike-based backpropagation rules have been proposed that perform end-to-end gradient descent training on spike data. In spike-based backpropagation methods, the non-differentiability of the spiking neuron is handled by either approximating the spiking neuron model as continuous and differentiable \citep{huh2018gradient} or by defining a surrogate gradient as a continuous approximation of the real gradient \citep{wu2018spatio,bellec2018long,neftci2019surrogate}. Spike-based SNN training reduces the overall latency by ${\sim}10{\times}$ (for instance, $200-250$ time steps required to process an input \citep{lee2019enabling}) but requires more training effort (in terms of total training iterations) than conversion approaches. A single feed-forward pass in ANN corresponds to multiple forward passes in SNN which is proportional to the number of time steps. In spike-based backpropagation, the backward pass requires the gradients to be integrated over the total number of time steps that increases the computation and memory complexity. The multiple-iteration training effort with exploding memory requirement (for backward pass computations) has limited the applicability of spike-based backpropagation methods to small datasets (like CIFAR10) on simple few-layered convolutional architectures.

In this work, we propose a hybrid training technique which combines ANN-SNN conversion and spike-based backpropagation that reduces the overall latency as well as decreases the training effort for convergence. We use ANN-SNN conversion as an initialization step followed by spike-based backpropagation incremental training (that converges to optimal accuracy with few epochs due to the precursory initialization).  Essentially, our hybrid approach of taking a converted SNN and incrementally training it using backpropagation yields improved energy-efficiency as well as higher accuracy than a model trained from scratch with only conversion or only spike-based backpropagation.

In summary, this paper makes the following contributions:
\begin{itemize}
    \item We introduce a hybrid computationally-efficient training methodology for deep SNNs. We use the weights and firing thresholds of an SNN converted from an ANN as the initialization step for spike-based backpropagation. We then train this initialized network with spike-based backpropagation for few epochs to perform inference at a reduced latency or time steps.

    \item We propose a novel spike time-dependent backpropagation (STDB, a variant of standard spike-based backpropagation) that computes surrogate gradient using neuron's spike time. The parameter update is triggered by the occurrence of spike and the gradient is computed based on the time difference between the current time step and the most recent time step the neuron generated an output spike. This is motivated from the Hebb's principle which states that the plasticity of a synapse is dependent on the spiking activity of the neurons connected to the synapse.

    \item Our hybrid approach with the novel surrogate gradient descent allows training of large-scale SNNs without exploding memory required during spike-based backpropagation. We evaluate our hybrid approach on large SNNs (VGG, ResNet-like architectures) on Imagenet, CIFAR datasets and show near iso-accuracy compared to similar ANNs and converted SNNs at lower compute cost and energy.
\end{itemize}

\begin{wrapfigure}{r}{0.5\textwidth}
    \centering
    \includegraphics[width=0.5\textwidth]{gradient}
    \caption{Surrogate gradient of the spiking neuron activation function (Eq. \ref{eq:gradient}). $\alpha = 0.3, \beta = 0.01$. The gradient is computed for each neuron and $\Delta t$ defines the time difference between current simulation time and the last spike time of the neuron. For example, if a neuron spikes at $t_s = 12$ its gradient will be maximum at $t=12 (\Delta t = 0)$ and gradually decrease for later time steps. If the same neuron spikes later at $t_s = 24$ its previous spike history will be overwritten and the gradient computation for $t=24$ onward will only consider the most recent spike. This avoids the overhead of storing all the spike history in memory.}
    \vspace{-40pt}
    \label{fig:gradient}
\end{wrapfigure}

\section{Spike Timing Dependent Backpropagation (STDB)}\label{sec:stdb}

In this section, we describe the spiking neuron model, derive the equations for the proposed surrogate gradient based learning, present the weight initialization method for SNN, discuss the constraints applied for ANN-SNN conversion, and summarize the overall training methodology.

\subsection{Leaky Integrate and Fire (LIF) Neuron Model}\label{subsec:lif_model}

The neuron model defines the dynamics of the neuron's internal state and the trigger for it to generate a spike. The differential equation
\begin{equation}\label{eq:neuron_differential}
    \tau\frac{dU}{dt} = -(U - U_{rest}) + RI
\end{equation}
is widely used to characterize the leaky-integrate-and-fire (LIF) neuron model where, \(U \) is the internal state of the neuron referred as the membrane potential, \( U_{rest}\) is the resting potential, \( R\) and \( I\) are the input resistance and the current, respectively. The above equation is valid when the membrane potential is below the threshold value (\( V\)). The neuron geneartes an output spike when $U{\geqslant}V$ and \( U\) is reduced to the reset potential. This representation is described in continuous domain and more suitable for biological simulations. We modify the equation to be evaluated in a discrete manner in the Pytorch framework \citep{wu2018spatio}. The iterative model for a single post-neuron is described by
\begin{equation}\label{eq:neuron_iterative}
    u^{t}_i = \lambda u^{t-1}_i + \sum_j w_{ij}o_j^t - vo_i^{t-1}
\end{equation}
\begin{equation}\label{eq:thresholding}
    o_i^{t-1} =
    \begin{cases}
    1, & \text{if}\ u_i^{t-1} > v \\
    0, & \text{otherwise}
    \end{cases}
\end{equation}

where \( u\) is the membrane potential, subscript \( i\) and \( j\) represent the post- and pre-neuron, respectively, superscript \( t\) is the time step, \( \lambda\) is a constant (\( <1\)) responsible for the leak in membrane potential, \( w\) is the weight connecting the pre- and post-neuron, \( o\) is the binary output spike, and \( v\) is the firing threshold potential. The right hand side of Equation \ref{eq:neuron_iterative} has three terms: the first term calculates the leak in the membrane potential from the previous time step, the second term integrates the input from the previous layer and adds it to the membrane potential, and the third term which is outside the summation reduces the membrane potential by the threshold value if a spike is generated. This is known as soft reset as the membrane potential is lowered by $v$ compared to hard reset where the membrane potential is reduced to the reset value. Soft reset enables the spiking neuron to carry forward the excess potential above the firing threshold to the following time step, thereby minimizing information loss.

\begin{minipage}[t]{0.48\textwidth}
\vspace{-20pt}
\begin{algorithm}[H]
\caption{ANN-SNN conversion: initialization of weights and threshold voltages}\label{algo:ann-snn_conversion}
\KwIn{Trained ANN model ($A$), SNN model ($N$), Input ($X$)}
// Copy ann weights to snn \\
\For{l=1 \KwTo L}{
$N_l.W \leftarrow A_l.W$
}
// Initialize threshold voltage to $0$ \\
$V \leftarrow [0, \cdots , 0]_{L-1}$ \\
\For{l=1 \KwTo L-1}{
    $v \leftarrow 0$ \\
    \For{t=1 to T}{
        $O^t_0 \leftarrow PoissonGenerator(X)$ \\
        \For{k=1 \KwTo l}{
            \If{$k<l$}{
                // Forward (Algorithm \ref{algo:snn_training}) \\
                }
            \Else{
                // Pre-nonlinearity ($A$) \\
                $A \leftarrow N_l(O_{k-1}^t)$ \\
                \If{$max(A) > v$}{
                    $v \leftarrow max(A)$
                }
            }
        }
    }
    $V[l] \leftarrow v$
}
\end{algorithm}
\end{minipage}
\hfill
\begin{minipage}[t]{0.48\textwidth}
\vspace{-20pt}
\begin{algorithm}[H]
\caption{Initialize the neuron parameters. Membrane potential ($U$), last spike time ($S$), dropout mask ($M$). The initialization is performed once for every mini-batch.}
\KwIn{Input($X$), network model($N$)}\label{algo:neuron_init}

$b\_size = X.b\_size$ \\
$h = X.height$ \\
$w = X.width$ \\
\For{l=1 \KwTo L}{
    \If{$isintance (N_l, Conv)$}{
    $U_l = zeros(b\_size, N_l.out, h, w)$ \\
    $S_l = ones(b\_size,N_l.out, h, w) * (-1000)$
    }
     \ElseIf{$isintance (N_l, Linear$)}{
    $U_l = zeros(b\_size, N_l.out)$ \\
    $S_l = ones(b\_size, N_l.out) * (-1000)$
    }
    \ElseIf{$isintance (N_l, Dropout$)}{
    // Generate the dropout map that will be fixed for all time steps \\
    $M_l = N_l(ones(U_{l-1}.shape))$
    }
    \ElseIf{$isintance (N_l, AvgPool)$}{
    // Reduce the width and height after average pooling layer\\
    $h = h//kernel\_size$ \\
    $w = w//kernel\_size$

    }
}
\end{algorithm}
\end{minipage}

\subsection{Spike Timing Dependent Backpropagation (STDB) Learning Rule}\label{subsec:spike_timing_dependent_backpropagation_learning_rule}
The neuron dynamics (Equation~\ref{eq:neuron_iterative}) show that the neuron's state at a particular time step recurrently depends on its state in previous time steps. This introduces implicit recurrent connections in the network \citep{neftci2019surrogate}. Therefore, the learning rule has to perform the temporal credit assignment along with the spatial credit assignment. Credit assignment refers to the process of assigning credit or blame to the network parameters according to their contribution to the loss function. Spatial credit assignment identifies structural network parameters (like weights), whereas temporal credit assignment determines which past network activities contributed to the loss function.  Gradient-descent learning solves both credit assignment problem:
spatial credit assignment is performed by distributing error spatially across all layers using the chain rule of derivatives, and temporal credit assignment is done by unrolling the network in time and performing backpropagation through time (BPTT) using the same chain rule of derivatives \citep{werbos1990backpropagation}. In BPTT, the network is unrolled for all time steps and the final output is computed as the sum of outputs from each time step. The loss function is defined on the summed output.

The dynamics of the neuron in the output layer is described by Equation~(\ref{eq:neuron_output_integrate}), where the leak part is removed (\( \lambda = 1\)) and the neuron only integrates the input without firing. This eliminates the difficulty of defining the loss function on spike count \citep{lee2019enabling}.
\begin{equation}\label{eq:neuron_output_integrate}
    u^{t}_i = u^{t-1}_i + \sum_j w_{ij}o_j
\end{equation}

The number of neurons in the output layer is the same as the number of categories in the classification task. The output of the network is passed through a softmax layer that outputs a probability distribution. The loss function is defined as the cross-entropy between the true output and the network's predicted distribution.
\begin{equation}
    L = - \sum_i y_i log(p_i)
\end{equation}
\begin{equation}
    p_i = \frac{e^{u_i^T}}{\sum_{k=1}^{N} e^{u_k^T}}
\end{equation}

\( L\) is the loss function, \( y\) the true output, \( p\) the prediction, $T$ the total number of time steps, \( u^T \) the accumulated membrane potential of the neuron in the output layer from all time steps, and \( N\) the number of categories in the task. For deeper networks and large number of time steps the truncated version of the BPTT algorithm is used to avoid memory issues. In the truncated version the loss is computed at some time step $t'$ before T based on the potential accumulated till $t'$. The loss is backpropagated to all layers and the loss gradients are computed and stored. At this point, the history of the computational graph is cleaned to save memory. The subsequent computation of loss gradients at later time steps ($2t',3t',...T$) are summed together with the gradient at $t'$ to get the final gradient. The optimizer updates the parameters at $T$ based on the sum of the gradients. Gradient descent learning has the objective of minimizing the loss function. This is achieved by backpropagating the error and updating the parameters opposite to the direction of the derivative.  The derivative of the loss function w.r.t. to the membrane potential of the neuron in the final layer is described by,
\begin{equation}
    \frac{\partial L}{\partial u_i^T} = p_i - y_i
\end{equation}

\begin{algorithm}[t]
\caption{Training an SNN with surrogate gradient computed with spike timing. The network is composed of $L$ layers. The training proceeds with mini-batch size ($batch\_size$)}\label{algo:snn_training}
\KwIn{Mini-batch of input ($X$) - target ($Y$) pairs, network model ($N$), initial weights ($W$), threshold voltage ($V$)}

$U, S, M = InitializeNeuronParameters(X)$ [Algorithm \ref{algo:neuron_init}] \\
// Forward propagation \\
\For{t=1 \KwTo T}{
   $O_0^t = PoissonGenerator(X)$ \\
   \For{l=1 \KwTo L-1}{
        \If{$isintance (N_l, [Conv, Linear])$}{
            // accumulate the output of previous layer in $U$, soft reset when spike occurs \\
            $U_l^t = \lambda U_l^{t-1} + W_lO_{l-1}^t - V_l*O^{t-1}_l$ \\
            // generate the output ($+1$) if $U$ exceeds $V$ \\
            $O^t_l = STDB (U_l^t, V_l, t)$ \\
            // store the latest spike times for each neuron\\
            $S_l^t[O_l^t==1] = t $
            }
        \ElseIf{$isintance (N_l, AvgPool)$}{
            $O^t_l = N_l(O_{l-1}^t)$
            }
        \ElseIf{$isintance (N_l, Dropout)$}{
        $O^t_l = O_{l-1}^t * M_l$
        }
    }
    $U_L^t = \lambda U_L^{t-1} + W_LO_{L-1}^t$
}
// Backward Propagation \\
Compute $\frac{\partial L}{\partial U_L}$ from the cross-entropy loss function using BPTT \\
\For{t=T \KwTo 1}{
    \For{l=L-1 \KwTo 1}{
       Compute $\frac{\partial L}{\partial O^t_l}$ based on if $N_l$ is linear, conv, pooling, etc. \\
       $\frac{\partial L}{\partial U_l^t} = \frac{\partial L}{\partial O^t_l} \frac{\partial O^t_l}{\partial U_l^t} = \frac{\partial L}{\partial O^t_l} * \alpha e^{-\beta S_l^t}$
    }
}
\end{algorithm}

To compute the gradient at current time step, the membrane potential at last time step ($u^{t-1}_i$ in Equation~\ref{eq:neuron_output_integrate}) is considered as an input quantity. Therefore, gradient descent updates the network parameters \( W_{ij}\) of the output layer as,
\begin{equation}
    W_{ij} = W_{ij} - \eta\Delta W_{ij}
\end{equation}
\begin{equation}
    \Delta W_{ij} = \mathlarger{\sum_t} \frac{\partial L}{\partial W_{ij}^t} =  \mathlarger{\sum_t} \frac{\partial L}{\partial u_i^T} \frac{\partial u_i^T}{\partial W_{ij}^t} = \frac{\partial L}{\partial u_i^T}  \mathlarger{\sum_t} \frac{\partial u_i^T}{\partial W_{ij}^t}
\end{equation}

where \( \eta \) is the learning rate, and $W_{ij}^t$ represents the copy of the weight used for computation at time step $t$.  In the output layer the neurons do not generate a spike, and hence, the issue of non-differentiability is not encountered. The update of the hidden layer parameters is described by,
\begin{equation}
    \Delta W_{ij} = \mathlarger{\sum_t}\frac{\partial L}{\partial W_{ij}^t} = \mathlarger{\sum_t}\frac{\partial L}{\partial o_i^t} \frac{\partial o_i^t}{\partial u_i^t} \frac{\partial u_i^t}{\partial W_{ij}^t}
\end{equation}
where \( o_i^t \) is the thresholding function (Equation~\ref{eq:thresholding}) whose derivative w.r.t to \( u_i^t\) is zero everywhere and not defined at the time of spike.
The challenge of discontinuous spiking nonlinearity is resolved by introducing a surrogate gradient which is the continuous approximation of the real gradient.
\begin{equation}\label{eq:gradient}
    \frac{\partial o_i^t}{\partial u_i^t} = \alpha e^{- \beta \Delta t}
\end{equation}
where \( \alpha\) and \( \beta\) are constants, \( \Delta t\) is the time difference between the current time step ($t$) and the last time step the post-neuron generated a spike ($t_s$). It is an integer value whose range is from zero to the total number of time steps (\( T\)).
\begin{equation}\label{eq:delta_t_bound}
    \Delta t = (t - t_s) ,
    0<\ \Delta t<T \ , \ \Delta t \ \epsilon \ \mathbb{Z}
\end{equation}

The values of $\alpha$ and $\beta$ are selected depending on the value of $T$. If $T$ is large $\beta$ is lowered to reduce the exponential decay so a spike can contribute towards gradients for later time steps. The value of $\alpha$ is also reduced for large $T$ because the gradient can propagate through many time steps. The gradient is summed at each time step and thus a large $\alpha$ may lead to exploding gradient. The surrogate gradient can be pre-computed for all values of \( \Delta t\) and stored in a look-up table for faster computation. The parameter updates are triggered by the spiking activity but the error gradients are still non-zero for time steps following the spike time. This enables the algorithm to avoid the `dead neuron' problem, where no learning happens when there is no spike. Fig. \ref{fig:gradient} shows the activation gradient for different values of $\Delta t$, the gradient decreases exponentially for neurons that have not been active for a long time. In Hebbian models of biological learning, the parameter update is activity dependent. This is experimentally observed in spike-timing-dependent plasticity (STDP) learning rule which modulates the weights for pair of neurons that spike within a time window~\citep{song2000competitive}.
\begin{wrapfigure}{R}{0.35\textwidth}
\centering
\includegraphics[width=0.35\textwidth]{Figures/resnet.pdf}
\caption{Residual architecture for SNN}
\vspace{-40pt}
\label{fig:resnet}
\end{wrapfigure}
\section{SNN Weight Initialization}\label{sec:snn_weight_initialization}

A prevalent method of constructing SNNs for inference is ANN-SNN conversion \citep{diehl2015fast, sengupta2019going}. Since the network is trained with analog activations it does not suffer from the non-differentiablity issue and can leverage the training techniques of ANNs. The conversion process has a major drawback: it suffers from long inference latency (${\sim} 2500$ time steps) as mentioned in Section~\ref{sec:introduction}. As there is no provision to optimize the parameters after conversion based on spiking activity, the network can not leverage the temporal information of the spikes. In this work, we propose to use the conversion process as an initialization technique for STDB. The converted weights and thresholds serve as a good initialization for the optimizer and the STDB learning rule is applied for temporal and spatial credit assignment.

Algorithm \ref{algo:ann-snn_conversion} explains the ANN-SNN conversion process. The threshold voltages in SNN needs to be adjusted based on the ANN weights. \citet{sengupta2019going} showed two ways to achieve this: weight-normalization and threshold-balancing. In weight-normalization the weights are scaled by a normalization factor and threshold is set to 1, whereas in threshold-balancing the weights are unchanged and the threshold is set to the normalization factor. Both have a similar effect and either can be used to set the threshold. We employ the threshold-balancing method and the normalization factor is calculated as the maximum output of the corresponding convolution/linear layer in SNN. The maximum is calculated over a mini-batch of input for all time steps.

There are several constraints imposed on training the ANN for the conversion process \citep{sengupta2019going, diehl2015fast}. The neurons are trained without the bias term because the bias term in SNN has an indirect effect on the threshold voltage which increases the difficulty of threshold balancing and the process becomes more prone to conversion loss. The absence of bias term eliminates the use of Batch Normalization \citep{ioffe2015batch} as a regularizer in ANN since it biases the input of each layer to have zero mean. As an alternative, Dropout \citep{srivastava2014dropout} is used as a regularizer for both ANN and SNN training. The implementation of Dropout in SNN is further discussed in Section \ref{sec:overall_proposed_algorithm}. The pooling operation is widely used in ANN to reduce the convolution map size. There are two popular variants: max pooling and average pooling \citep{boureau2010theoretical}. Max (Average) pooling outputs the maximum (average) value in the kernel space of the neuron's activations. In SNN, the activations are binary and performing max pooling will result in significant information loss for the next layer, so we adopt the average pooling for both ANN and SNN \citep{diehl2015fast}.

\section{Network Architectures} \label{sec:network_architectures}
In this section, we describe the changes made to the VGG \citep{simonyan2014very} and residual architecture \citep{he2016deep} for hybrid learning and discuss the process of threshold computation for both the architectures.

\subsection{VGG Architecture} \label{subsec:vgg_architecture}
The threshold balancing is performed for all layers except the input and output layer in a VGG architecture. For every hidden convolution/linear layer the maximum input\footnote{input to a neuron is the weighted sum of spkies from pre-neurons $\sum_j {w_{ij}o_j}$} to the neuron is computed over all time steps and set as threshold for that layer. The threshold assignment is done sequentially as described in Algorithm \ref{algo:ann-snn_conversion}. The threshold computation for all layers can not be performed in parallel (in one forward pass) because in the forward method (Algorithm \ref{algo:snn_training}) we need the threshold at each time step to decide if the neuron should spike or not.

\subsection{Residual Architecture}\label{subsec:residual_architecture}
Residual architectures introduce shortcut connections between layers that are not next to each other. In order to minimize the ANN-SNN conversion loss various considerations were made by \citet{sengupta2019going}. The original residual architecture proposed by \citet{he2016deep} uses an initial convolution layer with wide kernel ($7{\times}7$, stride $2$). For conversion, this is replaced by a pre-processing block consisting of a series of three convolution layer ($3{\times}3$, stride $1$) with dropout layer in between (Fig.~\ref{fig:resnet}). The threshold balancing mechanism is applied to only these three layers and the layers in the basic block have unity threshold.

\begin{table}[t]
\caption{Classification results (Top-1) for CIFAR10, CIFAR100 and ImageNet data sets. Column-1 shows the network architecture. Column-2 shows the ANN accuracy when trained under the constraints as described in Section \ref{sec:snn_weight_initialization}. Column-3 shows the SNN accuracy for $T= 2500$ when converted from a ANN with threshold balancing. Column-4 shows the performance of the same converted SNN with lower time steps and adjusted thresholds. Column-5 shows the performance after training the Column-4 network with STDB for less than 20 epochs.}\label{table:results}
\renewcommand{\arraystretch}{1.4}
\begin{center}
\begin{tabular}{c|>{\centering\arraybackslash}p{1.5cm}|>{\centering\arraybackslash}p{2cm}|>{\centering\arraybackslash}p{3cm}|>{\centering\arraybackslash}p{3.5cm}}
\hline
Architecture & ANN   & ANN-SNN Conversion ($T=2500$) & ANN-SNN Conversion (reduced time steps) & Hybrid Training (ANN-SNN Conversion + STDB)   \\ \hline\hline
\multicolumn{5}{c}{CIFAR10}                                      \\ \hline\hline
VGG5       & $87.88\%$ & $87.64\%$          & $84.56\%$ ($T=75$)     & \textbf{$86.91\%$ ($T=75$)} \\ \hline
VGG9     & $91.45\%$ & $90.98\%$           & $87.31\%$ ($T=100$)      & \textbf{$90.54\%$ ($T=100$)} \\ \hline
VGG16     & $92.81\%$ & $92.48\%$             & $90.2\%$ ($T=100$)      & \textbf{$91.13\%$ ($T=100$)} \\ \hline
ResNet8   & $91.35\%$ & $91.12\%$            & $89.5\%$ ($T=200$)      & \textbf{$91.35\%$ ($T=200$)} \\ \hline
ResNet20  & $93.15\%$ & $92.94\%$            & $91.12\%$ ($T=250$)     & \textbf{$92.22\%$ ($T=250$)} \\ \hline\hline
\multicolumn{5}{c}{CIFAR100}                                     \\ \hline\hline
VGG11     & $71.21\%$ & $70.94\%$             & $65.52\%$ ($T=125$)     & \textbf{$67.87\%$ ($T=125$)} \\ \hline\hline
\multicolumn{5}{c}{ImageNet}                                     \\
\hline\hline
ResNet34      & $70.2\%$ & $65.1\%$             & $56.87\%$ ($T=250$)     & \textbf{$61.48\%$ ($T=250$)} \\ \hline
VGG16     & $69.35\%$ & $68.12\%$           & $62.73\%$ ($T=250$)       & \textbf{$65.19\%$ ($T=250$)} \\
\hline
\end{tabular}
\end{center}
\vspace{-20pt}
\end{table}

\section{Overall Training Algorithm}\label{sec:overall_proposed_algorithm}
Algorithm \ref{algo:ann-snn_conversion} defines the process to initialize the parameters (weights, thresholds) of SNN based on ANN-SNN conversion. Algorithm \ref{algo:neuron_init} and \ref{algo:snn_training} show the mechanism of training the SNN with STDB. Algorithm~\ref{algo:neuron_init} initializes the neuron parameters for every mini-batch, whereas Algorithm \ref{algo:snn_training} performs the forward and backward propagation and computes the credit assignment.
 The threshold voltage for all neurons in a layer is same and is not altered in the training process. For each dropout layer we initialize a mask ($M$) for every mini-batch of inputs. The function of dropout is to randomly drop a certain number of inputs in order to avoid overfitting. In case of SNN, inputs are represented as a spike train and we want to keep the dropout units same for the entire duration of the input. Thus, a random mask ($M$) is initialized (Algorithm \ref{algo:neuron_init}) for every mini-batch and the input is element-wise multiplied with the mask to generate the output of the dropout layer \citep{lee2019enabling}. The Poisson generator function outputs a Poisson spike train with rate proportional to the pixel value in the input. A random number is generated at every time step for each pixel in the input image. The random number is compared with the normalized pixel value and if the random number is less than the pixel value an output spike is generated. This results in a Poisson spike train with rate equivalent to the pixel value if averaged over a long time. The weighted sum of the input is accumulated in the membrane potential of the first convolution layer. The STDB function compares the membrane potential and the threshold of that layer to generate an output spike. The neurons that output a spike their corresponding entry in $S$ is updated with current time step ($t$). The last spike time is initialized with a large negative number (Algorithm~\ref{algo:neuron_init}) to denote that at the beginning the last spike happened at negative infinity time. This is repeated for all layers until the last layer. For last layer the inputs are accumulated over all time steps and passed through a softmax layer to compute the multi-class probability. The cross-entropy loss function is defined on the output of the softmax and the weights are updated by performing the temporal and spatial credit assignment according to the STDB rule.

\section{Experiments} \label{sec:experiments}
We tested the proposed training mechanism on image classification tasks from CIFAR~\citep{krizhevsky2009learning} and ImageNet~\citep{imagenet_cvpr09} datasets. The results are summarized in Table \ref{table:results}.
CIFAR10: The dataset consists of labeled $60,000$ images of $10$ categories divided into training ($50,000$) and testing ($10,000$) set. The images are of size $32{\times}32$ with RGB channels.

CIFAR100: The dataset is similar to CIFAR10 except that it has 100 categories.

ImageNet: The dataset comprises of labeled high-resolution $1.2$ million training images and $50,000$ validation images with $1000$ categories.

\section{Energy-Delay Product Analysis of SNNs}
A single spike in an SNN consumes a constant amount of energy \citep{cao2015spiking}. The first order analysis of energy-delay product of an SNN is dependent on the number of spikes and the total number of time steps. Fig.~\ref{fig:average_spike} shows the average number of spikes in each layer when evaluated for $1500$ samples from CIFAR10 testset for VGG16 architecture. The average is computed by summing all the spikes in a layer over $100$ time steps and dividing by the number of neurons in that layer. For example, the average number of spikes in the $10^{th}$ layer is $5.8$ for both the networks, which implies that over a $100$ time step period each neuron in that layer spikes $5.8$ times on average over all input samples. Higher spiking activity corresponds to lower energy-efficiency. The average number of spikes is compared for a converted SNN and SNN trained with conversion-and-STDB. The SNN trained with conversion-and-STDB has $1.5{\times}$ less number of average spikes over all layers under iso conditions (time steps, threshold voltages, inputs, etc.) and achieves higher accuracy compared to the converted SNN. The converted SNNs when simulated for larger time steps further degrade the energy-delay product with minimal increase in accuracy~\citep{sengupta2019going}.

\begin{figure}[t]
\centering
\includegraphics[width=0.7\textwidth]{Figures/average_spike.pdf}
\vspace{-10pt}
\caption{Average number of spikes for each layer in a VGG16 architecture for purely converted SNN and SNN trained with hybrid technique. The converted SNN and SNN trained with hybrid technique achieve an accuracy of 89.20\% and 91.87\%, respectively, for the randomly selected 1500 samples from the test set. Both the networks were inferred for 100 time steps and `v' represents the threshold voltage for each layer obtained during the conversion process (Algorithm \ref{algo:ann-snn_conversion}).}
\vspace{-10pt}
\label{fig:average_spike}
\end{figure}

\section{Related Work}\label{sec:realted_work}
\iffalse
ANN-SNN conversion - \citet{cao2015spiking, hunsberger2015spiking, sengupta2019going, diehl2015fast} \\
\fi
\cite{bohte2000spikeprop} proposed a method to directly train on SNN by keeping track of the membrane potential of spiking neurons only at spike times and backpropagating the error at spike times based on only the membrane potential. This method is not suitable for networks with sparse activity due to the `dead neuron' problem: no learning happens when the neurons do not spike. In our work, we need one spike for the learning to start but gradient contribution continues in later time steps as shown in Fig.~\ref{fig:gradient}. \citet{zenke2018superspike} derived a surrogate gradient based method on the membrane potential of a spiking neuron at a single time step only. The error was backpropagated at only one time step and only the input at that time step contributed to the gradient. This method neglects the effect of earlier spike inputs. In our approach, the error is backpropagated for every time step and the weight update is performed on the gradients summed over all time steps. \citet{shrestha2018slayer} proposed a gradient function similar to the one proposed in this work. They used the difference between the membrane potential and the threshold to compute the gradient compared to the difference in spike timing used in this work. The membrane potential is a continuous value whereas the spike time is an integer value bounded by the number of time steps. Therefore, gradients that depend on spike time can be pre-computed and stored in a look-up table for faster computation. They evaluated their approach on shallow architectures with two convolution layer for MNIST dataset. In this work, we trained deep SNNs with multiple stacked layers for complex calssification tasks. \citet{wu2018spatio} performed backpropagation through time on SNN with a surrogate gradient defined on the membrane potential. The surrogate gradient was defined as piece-wise linear or exponential function of the membrane potential. The other surrogate gradients proposed in the literature are all computed on the membrane potential \citep{neftci2019surrogate}. \citet{lee2019enabling} approximated the
neuron output as continuous low-pass filtered spike train. They used this approximated continuous value to perform backpropagation. Most of the works in the literature on direct training of SNN or conversion based methods have been evaluated on shallow architectures for simple classification problems. In Table \ref{table:result_comparison} we compare our model with the models that reported accuracy on CIFAR10 and ImageNet dataset. \citet{wu2019direct} achieved convergence in $12$ time steps by using a dedicated encoding layer to capture the input precision. It is beyond the scope of this work to compute the hardware and energy implications of such encoding layer. Our model performs better than all other models at far fewer number of time steps.

\iffalse
Linear surrogate gradient: \citet{courbariaux2016binarized, essera2016convolutional, bellec2018long}

SuperSpike: \citet{zenke2018superspike}

Exponential surrogate gradient SLayer: \citet{shrestha2018slayer}

Smooth spiking activation function: \citet{huh2018gradient,lee2019enabling}
\fi

\begin{table}[t]
\centering
\renewcommand{\arraystretch}{1.4}
\caption{Comparion of our work with other SNN models on CIFAR10 and ImageNet datasets}\label{table:result_comparison}
\begin{tabular}{p{2cm}|l|p{2cm}|p{1.7cm}|p{1.7cm}|p{1.4cm}}
\hline
Model     & Dataset  & Training Method    & Architecture  & Accuracy & Time-steps \\ \hline
\citet{hunsberger2015spiking}         & CIFAR10  & ANN-SNN Conversion  & 2Conv, 2Linear & $82.95\%$  & $6000$       \\ \hline
\citet{cao2015spiking}         & CIFAR10  & ANN-SNN Conversion  & 3Conv, 2Linear & $77.43\%$  & $400$        \\ \hline
\citet{sengupta2019going}         & CIFAR10  & ANN-SNN Conversion  & VGG16         & $91.55\%$  & $2500$       \\ \hline
\citet{lee2019enabling}         & CIFAR10  & Spiking BP         & VGG9          & $90.45\%$  & $100$        \\ \hline
\citet{wu2019direct}        & CIFAR10  & Surrogate Gradient & 5Conv, 2Linear & $90.53\%$  & $12$         \\ \hline
\textbf{This work} & \textbf{CIFAR10}  & \textbf{Hybrid Training} & \textbf{VGG16}         & \boldmath{$91.13\%$ $92.02\%$}\unboldmath  & \boldmath{$100$ \hspace{10pt} $200$}\unboldmath       \\ \hline\hline
\citet{sengupta2019going}         & ImageNet & ANN-SNN Conversion  & VGG16         & $69.96\%$  & $2500$       \\ \hline
\textbf{This work} & \textbf{ImageNet} & \textbf{Hybrid Training} & \textbf{VGG16}         & \boldmath{$65.19\%$}\unboldmath  & \boldmath{$250$}\unboldmath        \\ \hline
\end{tabular}
\vspace{-15pt}
\end{table}

\section{Conclusions} \label{sec:conclusion}

The direct training of SNN with backpropagation is computationally expensive and slow, whereas ANN-SNN conversion suffers from high latency. To address this issue we proposed a hybrid training technique for deep SNNs. We took an SNN converted from ANN and used its weights and thresholds as initialization for spike-based backpropagation of SNN. We then performed spike-based backpropagation on this initialized network to obtain an SNN that can perform with fewer number of time steps. The number of epochs required to train SNN was also reduced by having a good initial starting point. The resultant trained SNN had higher accuracy and lower number of spikes/inference compared to purely converted SNNs at reduced number of time steps. The backpropagation through time was performed with surrogate gradient defined using neuron's spike time that captured the temporal information and helped in reducing the number of time steps. We tested our algorithm on CIFAR and ImageNet datasets and achieved state-of-the-art performance with fewer number of time steps.

\subsubsection*{Acknowledgments}
This work was supported in part by the National Science Foundation, in part by Vannevar Bush Faculty Fellowship, and in part by C-BRIC, one of six centers in JUMP, a Semiconductor Research Corporation (SRC) program sponsored by DARPA.


text:
\section{Introduction}
\label{sec1}

Discovering dynamical system models from data is critically important across science and engineering.
Traditionally, models are derived from first principles, although this approach may be prohibitively challenging in many fields, such as climate science, finance, and biology.
Fortunately, data-driven model discovery (i.e., system identification) is a rapidly developing field~\cite{Brunton2019book}, with a range of techniques including classic linear approaches~\cite{Nelles2013book,ljung2010arc}, dynamic mode decomposition (DMD)~\cite{Schmid2010jfm,Kutz2016book} and Koopman theory~\cite{Budivsic2012chaos,Mezic2013arfm,Williams2015jnls,klus2017data}, nonlinear autoregressive models~\cite{Akaike1969annals,Billings2013book}, neural networks~\cite{yang2018physics,Wehmeyer2018jcp,Mardt2018natcomm,vlachas2018data,pathak2018model,lu2019deepxde,Raissi2019jcp,Champion2019pnas,raissi2020science}, Gaussian process regression~\cite{Raissi2017arxiva,Raissi2017arxiv}, nonlinear Laplacian spectral analysis~\cite{Giannakis2012pnas}, diffusion maps~\cite{Yair2017pnas}, genetic programming~\cite{Bongard2007pnas,Schmidt2009science,Daniels2015naturecomm}, and sparse regression~\cite{Brunton2016SINDy,Rudy2017PDE-find,Schaeffer2017prsa}, to highlight some of the recent developments.
Of particular note is a recent push towards \emph{parsimonious} modeling~\cite{Bongard2007pnas,Schmidt2009science,Brunton2016SINDy}, which favors Pareto-optimal models with the lowest complexity required to describe the observed data.  These models benefit from being interpretable, and they tend to generalize and prevent overfitting.
The sparse identification of nonlinear dynamics (SINDy) algorithm~\cite{Brunton2016SINDy} discovers parsimonious models through a sparsity-promoting optimization to select only a few model terms from a library of candidate functions.
SINDy has been widely adopted in the community~\cite{Sorokina2016oe,Dam2017pf,Schaeffer2017prsa,narasingam2018data,JC2018ConGalerkin,Quade2018chaos,boninsegna2018sparse,Loiseau2018jfm,Zhang2018arxiv,Hoffmann2018arxiv,mangan2019model,Thaler2019jcp,lai2019sparse,wu2018numerical,de2019discovery,Gelss2019mindy,Guang2018}, but it relies on the dynamics having a sparse representation in a pre-defined library, making it difficult to discover implicit dynamics and rational functions. The implicit-SINDy extension~\cite{Mangan2016ImSINDY} makes it possible to identify these implicit functions, although this algorithm is extremely sensitive to noise.
In this work, we develop a robust, parallel algorithm for the sparse identification of implicit dynamics, making it possible to explore entirely new classes of systems that were previously inaccessible.

Parsimonious modeling has a rich history, with many scientific advances being argued on the basis of Occam's razor, that the simplest model is likely the correct one.
SINDy exemplifies this principle, identifying a potentially nonlinear model with the fewest terms required to describe how the measurement data changes in time.
The basic idea behind SINDy may be illustrated on a one-dimensional system $\dot{x}=f(x)$; the general formulation for multidimensional dynamics will be described in the following sections.  An interpretable form of the nonlinear dynamics may be learned by writing the rate of change of the state of the system ${x}$ as a sparse linear combination of a few terms in a library of candidate functions, $\boldsymbol{\Theta}(x)=\begin{bmatrix}\theta_1(x) & \theta_2(x) & \dotsc & \theta_p(x)\end{bmatrix}$:
\begin{align}\label{Eq:BasicSINDy}
    \dot{x}(t) = {f}({x}(t)) \approx \boldsymbol{\Theta}({x}(t))\boldsymbol{\xi}.
\end{align}
where each $\theta_j(x)$ is prescribed candidate term (e.g. $x, x^2, \sin(x), \cdots$).
The derivative of the state and the library of candidate functions may both be computed from measured trajectory data.   It then remains to solve for a sparse vector $\boldsymbol{\xi}$ with nonzero entries $\xi_j$ indicating which functions $\theta_j(x)$ are active in characterizing the dynamics.
The resulting models strike a balance between accuracy and efficiency, and they are highly interpretable by construction.
In a short time, the SINDy algorithm has been extended to include inputs and control~\cite{Kaiser2018prsa}, to identify partial differential equations~\cite{Rudy2017PDE-find,Schaeffer2017prsa}, to incorporate physically relevant constraints~\cite{JC2018ConGalerkin}, to include tensor bases~\cite{Gelss2019mindy}, and to incorporate integral terms for denoising~\cite{Schaeffer2017IntegralSINDy,Reinbold2020pre}.
These extensions and its simple formulation in terms of a generalized linear model in \eqref{Eq:BasicSINDy} have resulted in SINDy being adopted in the fields of fluid mechanics~\cite{JC2018ConGalerkin,Loiseau2018jfm}, nonlinear optics~\cite{Sorokina2016oe}, plasma physics~\cite{Dam2017pf}, chemical reactions~\cite{Hoffmann2018arxiv,boninsegna2018sparse,narasingam2018data}, numerical methods~\cite{Thaler2019jcp}, and structural modeling~\cite{lai2019sparse}.

The generalized linear model in \eqref{Eq:BasicSINDy} does not readily lend itself to representing implicit dynamics and rational functions, which are not naturally expressible as sum of a few basis functions.
Instead, the implicit-SINDy algorithm~\cite{Mangan2016ImSINDY} reformulates the SINDy problem in an implicit form:
\begin{align}\label{Eq:ImplicitSimple}
    \boldsymbol{\Theta}(x,\dot{x})\boldsymbol{\xi} = 0.
\end{align}
This formulation is flexible enough to handle a much broader class of dynamics with rational function nonlinearities, such as $\dot{x} = {N(x)}/{D(x)}$ which may be rewritten as $\dot{x}D(x) + N(x) = 0$.
However, the sparsest vector $\boldsymbol{\xi}$ that satisfies \eqref{Eq:ImplicitSimple} is the trivial solution $\boldsymbol{\xi}=\mathbf{0}$.
Thus, the implicit-SINDy algorithm leverages a recent non-convex optimization procedure~\cite{Wright2009ieeetpami,Qu2014ADM} to find the sparsest vector $\boldsymbol{\xi}$ in the null space of $\boldsymbol{\Theta}(x,\dot{x})$, which differs from other approaches~\cite{zhu2015review,zhu2018control} that identify the rational dynamics. For even small amounts of noise, the dimension of the null space will become prohibitively large, making this approach extremely sensitive to noise and compromising the model discovery process.

This work develops an optimization and model selection framework that recasts implicit-SINDy as a convex problem, making it as noise robust as the original non-implicit SINDy algorithm and enabling the identification of implicit ODEs and PDEs that were previously inaccessible.
The key to making the implicit-SINDy algorithm robust is the realization that if we know even a single term in the dynamics, corresponding to a non-zero entry $\xi_j$, then we can rewrite \eqref{Eq:ImplicitSimple} in a non-implicit form
\begin{align}\label{Eq:SimplePI}
    \theta_j(x,\dot{x}) = \boldsymbol{\Theta}'(x,\dot{x})\boldsymbol{\xi}'
\end{align}
where $\boldsymbol{\Theta}'$ and $\boldsymbol{\xi}'$ have the $j$-th element removed.
Because none of these terms are known \emph{a priori}, we sweep through the  library, term by term, testing \eqref{Eq:SimplePI} for a sparse model that fits the data.
This procedure is highly parallelizable and provides critical information for model selection.
Our approach is related to the recent work of Zhang et al.~\cite{Guang2018}, which also makes the implicit problem more robust by testing candidate functions individually.
However, there are a number of key differences in the present approach.
Our work explicitly considers \emph{rational} nonlinearities to discover exceedingly complex implicit PDEs, such as a simplified model of the Belousov-Zhabotinsky (BZ) reaction.
Our framework also provides several new greedy algorithms, including parallel and constrained formulations.
We further extend this method to include the effect of control inputs, making it applicable to robotic systems~\cite{Koryakovskiy2018}, and we use this procedure to discover Hamiltonians.
Finally, our approach provides guidance on model selection, a comprehensive comparison with previous methods, and a careful analysis of noise robustness.

\section{Background}
\label{sec2}
We briefly introduce the full multidimensional SINDy and implicit-SINDy algorithms, which will provide a foundation for our robust implicit identification algorithm in Sec.~\ref{sec3}.

\subsection{Sparse Identification of Nonlinear Dynamics}
The goal of SINDy~\cite{Brunton2016SINDy} is to discover a dynamical system
\begin{align}\label{Eq:Dynamics}
    \frac{d}{dt}{\boldsymbol{x}}(t)=\boldsymbol{f}(\boldsymbol{x}(t)),
\end{align}
from time-series data of the state $\boldsymbol{x}(t)= [{x_{1}(t)},{\ldots}, {x_{n}(t)}]^{T} \in \R^{n}$.
We assume that the dynamics, encoded by the function $\boldsymbol{f}$, admit a sparse representation in a library of candidate functions:
\begin{equation}
    \boldsymbol{\Theta}(\boldsymbol{x}) = \begin{bmatrix} \theta_1(\boldsymbol{x}) & \theta_2(\boldsymbol{x}) & \cdots & \theta_p(\boldsymbol{x})\end{bmatrix}.
\end{equation}
Thus, each row equation in \eqref{Eq:Dynamics} may be written as
\begin{equation}
    \frac{d}{dt}x_k(t) = f_k(\boldsymbol{x}(t)) \approx \boldsymbol{\Theta}(\boldsymbol{x})\boldsymbol{\xi}_k
\end{equation}
where $\boldsymbol{\xi}_k$ is a sparse vector, indicating which terms are active in the dynamics.

We determine the nonzero entries of $\boldsymbol{\xi}_k$ through sparse regression based on trajectory data.
The time-series data is arranged into a matrix ${\mathbf{X} = \begin{bmatrix} \boldsymbol{x}(t_1) & \boldsymbol{x}(t_2) & \cdots & \boldsymbol{x}(t_m)\end{bmatrix}^T}$, and the associated time derivative matrix ${\dot{\mathbf{X}}=\begin{bmatrix}\dot{\boldsymbol{x}}(t_1)& \dot{\boldsymbol{x}}(t_2)& \cdots & \dot{\boldsymbol{x}}(t_m)\end{bmatrix}^T}$ is computed using an appropriate numerical differentiation scheme~\cite{Chartrand2011TVReg,Brunton2016SINDy,Brunton2019book}.
It is then possible to evaluate the library $\boldsymbol{\Theta}$ on trajectory data in $\mathbf{X}$ so that each column of $\boldsymbol{\Theta}(\mathbf{X})$ is a function $\theta_j$ evaluated on the $m$ snapshots in $\mathbf{X}$.

It is now possible to write the dynamical system in terms of a generalized linear model, evaluated on trajectory data:
\begin{equation}\label{Eq:SINDyMatrix}
    \dot{\mathbf{X}}=\boldsymbol{\Theta}(\mathbf{X}) \boldsymbol{\Xi}.
\end{equation}
There are several approaches to identify the sparse matrix of coefficients $\boldsymbol{\Xi}$, including sequentially thresholded least squares~(STLSQ)~\cite{Brunton2016SINDy,zhang2019convergence}, LASSO~\cite{Tibshirani1994}, sparse relaxed regularized regression~(SR3)~\cite{ZhengSR3,Champion2019SINDySR3}, stepwise sparse regression~(SSR)~\cite{boninsegna2018sparse}, and Bayesian approaches~\cite{Guang2018,Pan2016BayesianSINDy}.
It is possible to augment the library to include partial derivatives for the identification of partial differential equations (PDEs)~\cite{Rudy2017PDE-find,Schaeffer2017prsa}.
Similarly, it is possible to include external forcing terms in the library $\boldsymbol{\Theta}$, enabling the identification of forced and actively controlled systems~\cite{Kaiser2018prsa}. To alleviate the effect of noise, it is possible to reframe the SINDy problem in terms of an integral formulation~\cite{Schaeffer2017IntegralSINDy,Reinbold2020pre}. There are a number of factors that affect the robustness of SINDy, some of which are discussed in App.~\ref{App:Robustness}.

\subsection{Implicit Sparse Identification of Nonlinear Dynamics}
The implicit-SINDy algorithm~\cite{Mangan2016ImSINDY} extends SINDy to identify implicit differential equations
\begin{equation}\label{Eq:ImplicitDynamics}
    \boldsymbol{f}(\boldsymbol{x},\dot{\boldsymbol{x}})=0,
\end{equation}
and in particular, systems that include rational functions in the dynamics, such as chemical reactions and metabolic networks that have a separation of timescales.

The implicit-SINDy generalizes the library $\boldsymbol{\Theta}(\mathbf{X})$ in \eqref{Eq:SINDyMatrix} to include functions of $\boldsymbol{x}$ and $\dot{\boldsymbol{x}}$:
\begin{equation}\label{Eq:implicitSINDyMatrix}
    \boldsymbol{\Theta}(\mathbf{X},\dot{\mathbf{X}}) \boldsymbol{\Xi}=\mathbf{0}.
\end{equation}
However, this approach requires solving for a matrix $\boldsymbol{\Xi}$ whose columns $\boldsymbol{\xi}_k$ are sparse vectors in the null space of $\boldsymbol{\Theta}(\mathbf{X},\dot{\mathbf{X}})$.
This approach is non-convex, relying on the alternating directions method (ADM)~\cite{Mangan2016ImSINDY,Qu2014ADM}, and null space computations are highly ill-conditioned for noisy data~\cite{Gavish2014ieeetit,Mangan2016ImSINDY,Brunton2019book}, thus inspiring the current work and mathematical innovations.

\begin{figure}[t]
\vspace{-.2in}
    \centering
    \includegraphics[width=1\textwidth]{Figures/DL_SINDy.pdf}
    \caption{The illustration of the SINDy-PI algorithm on Michaelis-Menten dynamics. (a) The Michaelis-Menten system is simulated, and measurement data is provided to SINDy-PI. (b) Multiple possible left-hand side functions are tested at the same time. (c) The candidate model prediction error is calculated, and the best model is selected.}
    \label{Fig_DL_SINDy}
\end{figure}
\section{SINDy-PI: Robust Parallel Identification of Implicit Dynamics}\label{sec3}
We have developed the SINDy-PI (parallel, implicit) framework for the robust identification of implicit dynamics, bypassing the null space approach discussed in Sec.~2.2.
The idea is that if even a single term $\theta_j(\boldsymbol{x},\dot{\boldsymbol{x}})\in\boldsymbol{\Theta}(\boldsymbol{x},\dot{\boldsymbol{x}})$ in the dynamics \eqref{Eq:ImplicitDynamics} is known, it is possible to rewrite \eqref{Eq:implicitSINDyMatrix} as
\begin{equation}\label{Eq:SINDyPIMatrix}
    \theta_j(\mathbf{X},\dot{\mathbf{X}}) = \boldsymbol{\Theta}(\mathbf{X},\dot{\mathbf{X}}|\theta_j(\mathbf{X},\dot{\mathbf{X}}) \boldsymbol{\xi}_j,
\end{equation}
where $\boldsymbol{\Theta}(\mathbf{X},\dot{\mathbf{X}}|\theta_j(\mathbf{X},\dot{\mathbf{X}}))$ is the library $\boldsymbol{\Theta}(\mathbf{X},\dot{\mathbf{X}})$ with the $\theta_j$ column removed.
Equation~\eqref{Eq:SINDyPIMatrix} is no longer in implicit form, and the sparse coefficient matrix corresponding to the remaining terms may be solved for using previously developed SINDy techniques~\cite{Brunton2016SINDy,ZhengSR3,Champion2019SINDySR3,boninsegna2018sparse,Guang2018,Pan2016BayesianSINDy,Rudy2017PDE-find,Schaeffer2017prsa,Schaeffer2017IntegralSINDy,Reinbold2020pre}. In particular, we solve for a sparse coefficient vector $\boldsymbol{\xi}_j$ that minimizes the following loss function:
\begin{equation}\label{Eq:SINDyPILoss}
     \| \theta_j(\mathbf{X},\dot{\mathbf{X}}) - \boldsymbol{\Theta}(\mathbf{X},\dot{\mathbf{X}}|\theta_j(\mathbf{X},\dot{\mathbf{X}}) \boldsymbol{\xi}_j\|_2 +\beta\norm{\boldsymbol{\xi}_j}_0,
\end{equation}
where $\beta$ is the sparsity promoting parameter.
There are numerous relaxations of the non-convex optimization problem in~\eqref{Eq:SINDyPILoss}, for example the sequentially thresholded least-squares algorithm~\cite{Brunton2016SINDy}. Because there is no null space calculation, the resulting algorithm is considerably more robust to noise than the implicit-SINDy algorithm~\cite{Mangan2016ImSINDY}, i.e. we longer have to deal with an ill-conditioned null space problem.

In general, the entire point of SINDy is that the dynamics are not known ahead of time, and so it is necessary to test each candidate function $\theta_j$ until one of the models in \eqref{Eq:SINDyPIMatrix} admits a sparse and accurate solution.
When an incorrect candidate term is used, then the algorithm results in a dense (non-sparse) model $\boldsymbol{\xi}_j$ and an inaccurate model fit, and when a correct term is included, the algorithm identifies a sparse model $\boldsymbol{\xi}_j$ and an accurate model fit.
In this way, it is clear when the algorithm has identified the correct model.
Moreover, there is a wealth of redundant information, since each term in the correct model may be used as the candidate function on the left hand side, and the resulting models may be cross-referenced.
This approach is highly parallelizable, and each candidate term may be tested simultaneously in parallel.
The non-parallel formulation in \eqref{Eq:SINDyPIMatrix} was recently introduced by Zhang et al.~\cite{Guang2018} in the context of Bayesian regression, where they also make the implicit problem more robust by testing candidate functions individually; however, they do not consider dynamics with rational function nonlinearities or control inputs.
In this work, we extend the robust implicit formulation to identify several challenging implicit ODE and PDE systems with rational function nonlinearities, which are ubiquitous in engineering and natural systems, and systems with external forcing and control inputs.
We also introduce the parallel formulation and model selection frameworks.
Further, we will introduce a constrained optimization framework to simultaneously test all candidate functions.

\subsection{Model Selection}
For each candidate function in \eqref{Eq:SINDyPIMatrix}, we obtain one candidate model.
When the candidate function $\theta_j$ is not in the true dynamics, then the resulting coefficient vector $\boldsymbol{\xi}_j$ will not be sparse and there will be large prediction error.
In contrast, when a correct candidate function is selected, then we obtain a sparse coefficient vector $\boldsymbol{\xi}_j$ and small prediction error.
For an implicit dynamical system, there may be several different implicit equations that must be identified, resulting in several candidate functions that admit sparse models.
The sequentially thresholded least squares (STLSQ) algorithm that we use here, and whose convergence properties are considered by Zhang and Schaeffer~\cite{zhang2019convergence}, iteratively computes a least-squares solution to minimize $\|\theta_j(\mathbf{X},\dot{\mathbf{X}})-\boldsymbol{\Theta}(\mathbf{X},\dot{\mathbf{X}}|\theta_j(\mathbf{X},\dot{\mathbf{X}}))\boldsymbol{\xi}_j\|_2$ and then zeros out small entries in $\boldsymbol{\xi}_j$ that are below a set threshold $\lambda$.
This threshold $\lambda$ is a hyperparameter that must be tuned to select the model that most accurately balances accuracy and efficiency.
Thus, we must employ model selection techniques to identify the implicit models that best supports the data, while remaining as simple as possible.

There are several valid approaches to model selection.  To select a parsimonious yet accurate model we can also employ the Akaike information criterion~(AIC)~\cite{Akaike1998InformationTheory, Akaike1974StatisticalModelIdentification} and Bayesian information criterion~(BIC)~\cite{Schwarz1978EstimatingDomensionOfAModel}, as in~\cite{Mangan2017ModelSelection}.
It is also possible to sweep through the parameter $\lambda$ and candidate functions $\theta_j$, and then choose the Pareto optimal model from a \emph{family} of models on the Pareto front balancing accuracy and efficiency; this is the approach in the original SINDy work~\cite{Brunton2016SINDy} and in earlier work leveraging genetic programming to discover dynamics~\cite{Bongard2007pnas,Schmidt2009science}.
In this work, we take a different approach, selecting models based on performance on a test data set $\mathbf{X}_t$ that has been withheld for model validation to automate the model selection process.  For each threshold $\lambda$, the resulting model is validated on the test set $\mathbf{X}_t$, and the model with the lowest test error is selected.
One error function is the model fit:
\begin{equation}
    \label{eq27}
    \text{Error}={\frac{\norm{\theta_j(\mathbf{X}_t,\dot{\mathbf{X}}_t)-\boldsymbol{\Theta}(\mathbf{X}_t,\dot{\mathbf{X}}_t|\theta_j(\mathbf{X}_t,\dot{\mathbf{X}}_t)) \boldsymbol{\Xi}}_2}{\norm{\theta_j(\mathbf{X}_t,\dot{\mathbf{X}}_t)}_2}}.
\end{equation}
In practice, for rational dynamics, we select based upon the predicted derivative $\dot{\mathbf{X}}_t$:
\begin{equation}
    \label{eq_n1}
    \text{Error}={\frac{\norm{\dot{\mathbf{X}}_t-\dot{\mathbf{X}}^{\text{model}}_t}_2}{\norm{\dot{\mathbf{X}}_t}_2}}.
\end{equation}
For implicit dynamics where each state derivative may be written as a rational function
\begin{equation}
    \dot{x}_k = f_k(\boldsymbol{x}) = \frac{N_k(\boldsymbol{x})}{D_k(\boldsymbol{x})},
\end{equation}
then we restrict the candidate functions to $\theta_j(\boldsymbol{x},\dot{\boldsymbol{x}})=\dot{x}_k\theta_j(\boldsymbol{x})$ for some $\theta_j(\boldsymbol{x})\in\boldsymbol{\Theta}(\boldsymbol{x})$ to identify a separate sparse model for each $\dot{x}_k$.  Several candidate functions may provide accurate and sparse models.  These different models may further be cross-references to check that the same terms are being selected in each model, providing additional information for model selection and validation.

\subsection{Constrained Optimization Formulation}
\label{sec3.2}

\begin{figure}[t]
\vspace{-.2in}
    \centering
    \includegraphics[width=1\textwidth]{Figures/DL_SINDy_Constrained.pdf}
    \vspace{-.2in}
    \caption{Schematic illustrating the constrained formulation of the SINDy-PI algorithm. }
    \label{Fig_DL_SINDy_Con}
\end{figure}

In \eqref{Eq:SINDyPIMatrix} each candidate function was tested individually in a parallel optimization.
However, each of these individual equations may be combined into a single constrained system of equations
\begin{equation}
    \label{eq43}
    \boldsymbol{\Theta}(\mathbf{X},\dot{\mathbf{X}})=\boldsymbol{\Theta}(\mathbf{X},\dot{\mathbf{X}})\boldsymbol{\Xi} \quad\text{such that}\quad \boldsymbol{\Xi}_{jj}=0.
\end{equation}
We constrain $\boldsymbol{\Xi}$ to have zero entries on the diagonal, as shown in Fig.~\ref{Fig_DL_SINDy_Con}, which is the same as removing the candidate function from the library in the separate optimization problems in \eqref{Eq:SINDyPIMatrix}.
Without this constraint, the trivial solution $\boldsymbol{\Xi}=\mathbb{I}_{p\times p}$ will provide the sparsest $\boldsymbol{\Xi}$ and the most accurate model.
This may be written as a formal constrained optimization problem:
\begin{equation}
    \label{eq45}
    \begin{aligned}
        \min_{\boldsymbol{\Xi}} \quad & ||\boldsymbol{\Theta}(\mathbf{X},\dot{\mathbf{X}})-\boldsymbol{\Theta}(\mathbf{X},\dot{\mathbf{X}})\boldsymbol{\Xi}||_2 + \beta\|\boldsymbol{\Xi}\|_0,\\
        \textrm{s.t.} \quad & \textrm{diag}(\boldsymbol{\Xi})=\mathbf{0}.
    \end{aligned}
\end{equation}
This optimization is non-convex, although there are many relaxations that result in accurate and efficient proxy solutions~\cite{Tibshirani1994,Brunton2016SINDy,ZhengSR3}.
In this work, we will use sequentially thresholded least squares, so that any entry $\boldsymbol{\Xi}_{ij}<\lambda$ will be set to zero; the sparsity parameter $\lambda$ is a hyperparameter, and each column equation may require a different parameter $\lambda_j$.
The constrained formulation in \eqref{eq45} can be solved efficiently in modern optimization packages, and we use CVX~\cite{CVX_Software,Grant_Boyd2008_Nonsmooth_Convex_Programs}.
After solving \eqref{eq45} we have numerous candidate models, one for each column $\boldsymbol{\xi}_k$ of $\boldsymbol{\Xi}$, given by
\begin{equation}
    \label{eq46}
    \boldsymbol{\Theta}(\mathbf{X},\dot{\mathbf{X}})\boldsymbol{\xi}_j = 0.
\end{equation}
The sparse models that result in an accurate fit are candidate implicit models, and they may be assessed using the model selection approaches outlined above.
These various models may be cross-referenced for consistency, as the same models will have the same sparsity pattern.
This information can then be used to refine the library $\boldsymbol{\Theta}$, for example to only include the nonzero entries in the sparse columns of $\boldsymbol{\Xi}$.

\subsection{Noise Robustness}
\label{sec3.3}
We now compare the noise sensitivity of SINDy-PI and implicit-SINDy on the one-dimensional Michaelis–Menten model for enzyme  kinetics~\cite{Johnson2011MMK_Translate,Mangan2016ImSINDY,Cornish_Bowden2015MMK}, given by
\begin{equation}
    \label{s4eq1}
    \dot{x}=j_{x}-\frac{V_{\max } x}{K_{m}+x},
\end{equation}
where $x$ denotes the concentration of the substrate, $j_x$ denotes the influx of the substrate, $V_{max}$ denotes the maximum reaction time, and $K_m$ represents the concentration of half-maximal reaction.
We use the same parameters as in~\cite{Mangan2016ImSINDY}, with $j_x=0.6$, $V_{max}=1.5$, and $K_m=0.3$.
Figure~\ref{TVReg_k0} shows the result of the noise robustness of SINDy-PI and implicit-SINDy.
In this example, SINDy-PI is able to handle over $10^5$ more measurement noise than implicit-SINDy, while still accurately recovering the correct model.  Details are provided in App.~\ref{A1}, and key factors that limit robustness are discussed inn App.~\ref{App:Robustness}.

\begin{figure}[t]
    \vspace{-.2in}
    \centering
    \includegraphics[width=1\textwidth]{Figures/Noise_M3_k0.pdf}
    \caption{SINDy-PI and implicit-SINDy are compared on the Michaelis-Menten kinetics, where the structure error quantifies the number of terms in the model that are incorrectly added or deleted, compared with the true model. The derivative is computed by the total-variation regularization difference~(TVRegDiff)~\cite{Chartrand2011TVReg} on noisy state measurements. The violin plots show the cross-validated distribution of the number of incorrect terms across $30$ models. The green region indicates no structural difference between the identified model and the ground truth model. Details are provided in Appendix~\ref{App:NumericalExperiments}.}
    \label{TVReg_k0}
\end{figure}

\subsection{Data Usage}
\label{sec3.4}
The data required to correctly identify a model is a critical aspect when comparing SINDy-PI and implicit-SINDy.
Many experimental data sets are limited in volume, and thus our goal is to identify a model with as little data as possible.  In this section, we compare the SINDy-PI and implicit-SINDy methods on the challenging yeast glycolysis model~\cite{Schmidt2011YeastGlycolysis,Mangan2016ImSINDY} given by
\begin{subequations}
    \label{s4eq4}
    \begin{align}
       \dot{x}_{1} & =c_{1}+\frac{c_{2} x_{1} x_{6}}{1+c_{3} x_{6}^{4}}\label{s4eq4sub1},\\
       \dot{x}_{2} & =\frac{d_{1} x_{1} x_{6}}{1+d_{2} x_{6}^{4}}+d_{3} x_{2}-d_{4} x_{2} x_{7}\label{s4eq4sub2},\\
       \dot{x}_{3} & =e_{1} x_{2}+e_{2} x_{3}+e_{3} x_{2} x_{7}+e_{4} x_{3} x_{6}+f_{5} x_{4} x_{7}\label{s4eq4sub3},\\
       \dot{x}_{4} & =f_{1} x_{3}+e_{2} x_{4}+f_{3} x_{5}+f_{4} x_{3} x_{6}+f_{5} x_{4} x_{7}\label{s4eq4sub4},\\
       \dot{x}_{5} & =g_{1} x_{1}+g_{2} x_{5}\label{s4eq4sub5},\\
       \dot{x}_{6} & =h_{3} x_{3}+h_{5} x_{6}+h_{4} x_{3} x_{6}+\frac{h_{1} x_{1} x_{6}}{1+h_{2} x_{6}^{4}}\label{s4eq4sub6},
       \\ \dot{x}_{7} & =j_{1} x_{2}+j_{2} x_{2} x_{7}+j_{3} x_{4} x_{7}\label{s4eq4sub7}.
    \end{align}
\end{subequations}
 Equation~\eqref{s4eq4sub6} is the most challenging equation to discover in this system, and Fig.~\ref{Data_Length_Compare} compares the success rate of SINDy-PI and implicit-SINDy in identifying this equation.  SINDy-PI uses about $12$ times less data than the implicit-SINDy when identifying~\eqref{s4eq4sub6}. Details are provided in Appendices~\ref{B1} and~\ref{E1}.

\begin{figure}[t]
    \vspace{-.05in}
    \centering
    \includegraphics[width=1\textwidth]{Figures/DataLength.pdf}
    \caption{Success rate of SINDy-PI and implicit-SINDy identifying yeast glycolysis~\eqref{s4eq4sub6} with different percentage of training data. Each data usage percentage is randomly sampled from the entire data set composed of all trajectories.  The success rate is calculated by averaging the results of $20$ runs.}
    \label{Data_Length_Compare}
\end{figure}

\subsection{Comparison for Implicit PDE Identification}
\label{sec3.5}

\begin{figure}[t]

    \centering
    \includegraphics[width=1\textwidth]{Figures/Cpmpare_PDE_Find.pdf}
    \caption{Comparison of SINDy-PI and PDE-FIND on an implicit PDE problem given by the modified KdV equation \eqref{s5.3.1eq2}. As we increase $g_0$, the rational term begins to play a significant role in the system behavior. For small $g_0$, PDE-FIND compensates for the effect of the rational term by tuning the other coefficients. When $g_0$ is  large, PDE-FIND overfits the library. SINDy-PI, on the other hand, correctly identifies the rational term.}
    \label{PDE_KdV}
\end{figure}

We now investigate the ability of SINDy-PI to discover a PDE with rational terms, given by a modified KdV equation
\begin{equation}
    \label{s5.3.1eq2}
    u_t=-u_{xxx}-6uu_{x}-\gamma u+\frac{2g_0}{1+u},
\end{equation}
where $\gamma u$ is a loss term and $2g_0/(1+u)$ is a gain term.
We fix $\gamma=0.1$ and vary the value of $g_0$ from $0$ to $1$.
As $g_0$ increases, the implicit term gradually dominates the dynamics.
Figure~\ref{PDE_KdV} shows the results of SINDy-PI and PDE-FIND~\cite{Rudy2017PDE-find} for different values of $g_0$.
For large $g_0$, SINDy-PI is able to accurately identify the rational function term, while this is not possible for PDE-FIND, since this term is not in the library.
Details of the identification process are given in App.~\ref{D1}.

\section{Advanced Examples}
\label{sec5}
We will now demonstrate the SINDy-PI framework on several challenging examples, including the double pendulum, an actuated single pendulum on a cart, the Belousov-Zhabotinsky PDE, and the identification of conserved quantities.  All examples are characterized by rational nonlinearities, and we were unable to identify them using SINDy or implicit-SINDy, even in the absence of noise.

\subsection{Mounted Double Pendulum}
\label{secDoublePen}

\begin{figure}[t]
\vspace{-.25in}
    \centering
    \includegraphics[width=1\textwidth]{Figures/Double_Pendulum_DL.pdf}
    \vspace{-.15in}
    \caption{Schematic illustration of SINDy-PI identifying a mounted double pendulum system.}
    \label{DoublePenFig}
    \vspace{-.075in}
\end{figure}

In our first example, we use SINDy-PI to discover the equations of motion of a mounted double pendulum, shown in Fig.~\ref{DoublePenFig}.
The double pendulum is a classic example of chaotic dynamics~\cite{Graichen2007}, and was an original challenging example used to demonstrate the capabilities of genetic program for model discovery~\cite{Schmidt2009science}.
Correctly modeling the nonlinear dynamics is vital for accurate control~\cite{Graichen2007}.

We simulate the double pendulum dynamics, derived from the Euler-Lagrange equations, and use SINDy-PI to re-discover the dynamics from noisy measurements of the trajectory data.
The governing equations and SINDy-PI models are provided in App.~\ref{App:SINDyDoublePendulum}.
Because these dynamics have rational nonlinearities, the original SINDy algorithm is unable to identify the dynamics, making this a challenging test case.
The state vector is given by $\boldsymbol{x}= [\phi_1, \phi_2,  \dot{\phi}_1, \dot{\phi}_2]^T$, and the parameters of the simulation are given in App.~\ref{E1}.  The training data is generated from an initial condition $x_{\text{train}}=[\pi+1.2,\ \pi-0.6,\ 0,\ 0]^T$, simulated for $10$ seconds using a time step of $dt=0.001$ seconds. The validation data is generated from an initial condition $x_{\text{val}}=[\pi-1,\ \pi-0.4,\ 0.3,\ 0.4]^T$, simulated  for $3$ seconds with time step $dt=0.001$ seconds.

To test the robustness of SINDy-PI, we add Gaussian noise to both the training and validation data.  We test the resulting models using a new testing initial condition $x_{\text{test}}=[\pi+0.3,\ \pi-0.5,\ 0,\ 0]^T$.
We construct our library $\boldsymbol{\Theta}$ to include over 40 trigonometric and polynomial terms.
The most challenging part of this example is building a library with the necessary terms, without it growing too large. The library cannot be too extensive, or else the matrix $\boldsymbol{\Theta}$ becomes ill conditioned, making it sensitive to noise.
To reduce the library size, we use one piece of expert knowledge: the trigonometric terms should only consist of $\phi_1$ and $\phi_2$, the rotational angles of the pendula.

The candidate functions are chosen as a combination of state derivatives and trigonometric functions. Fig.~\ref{DoublePenFig} shows that SINDy-PI can identify the equations of motion for low noise.  For larger noise, SINDy-PI misidentifies the dynamics, although it still has short term prediction ability.

\subsection{Single Pendulum on a Cart}
\label{secSinglePen}
\begin{figure}[t]
    \vspace{-.25in}
    \centering
    \includegraphics[width=1\textwidth]{Figures/Single_Pendulum_DL.pdf}
    \vspace{-.175in}
    \caption{SINDy-PI is used to identify the single pendulum on a cart system. Control  is applied to the cart, and both the cart and pendulum states are measured. When the measurement noise is small, SINDy-PI can identify the correct structure of the model. }
    \vspace{-.075in}
    \label{Single_Pen}
\end{figure}

We now apply SINDy-PI to identify a fractional ODE problem with control input, given by the single pendulum on a cart in Fig.~\ref{Single_Pen}. SINDy has already been extended to include control inputs~\cite{Kaiser2018prsa}, although the original formulation doesn't accommodate rational functions.

The dynamics are derived from the Euler-Lagrange equations. All system parameters except for gravity are chosen to be $1$, as summarized in App.~\ref{E1}; the governing equations and SINDy-PI models are shown in App.~\ref{App:SINDyPendulumCart}.
The cart position is denoted by $s$.
The state vector is given by $\boldsymbol{x}= [\phi, s,  \dot{\phi}, \dot{s}]^T$.
The equations of motion are given by
\begin{subequations}
    \label{s5.2eq2}
    \begin{align}
        \frac{d}{dt}\phi&=\dot{\phi},\\
        \frac{d}{dt}s&=\dot{s},\\
        \frac{d}{dt}\dot{\phi}&=-\frac{(M+m)g\sin{(\phi)}+FL_1\cos{(\phi)}+mL_1^2\sin{(\phi)}\cos{(\phi)}\dot{\phi}^2}{L_1^2(M+m-m\cos{(\phi)}^2)},\\
        \frac{d}{dt}\dot{s}&=\frac{mL_1^2\sin{(\phi)}\dot{\phi}^2+FL_1+mg\sin{(\phi)}\cos{(\phi)}}{L_1(M+m-m\cos{(\phi)}^2)},
    \end{align}
\end{subequations}

Eq.~\eqref{s5.2eq2} is simulated with a time step of $dt=0.001$ to generate the training and testing data for model selection. The training data is generated using an initial condition $x_{\text{train}}=[0.3,\ 0,\ 1,\ 0]^T$ with the control input chosen as $F_{\text{train}}=-0.2+0.5\sin{(6t)}$, for time $t=0$ to $t=16$.
Similarly, the validation data is generated using an initial condition $x_{\text{val}}=[0.1,\ 0,\ 0.1,\ 0]^T$ with the control input chosen as $F_{\text{val}}=-1+\sin{(t)}+3\sin{(2t)}$, for time  $t=0$ to $t=2$.

 The library is constructed using a combination of trigonometric and polynomial terms. Around $50$ different basis functions are used for the library, and around $10$ terms are tested as candidate functions. We add Gaussian noise to all system states. We then test the SINDy-PI model on a testing initial condition $x_{\text{test}}=[\pi,\ 0,\ 0,\ 0]^T$ with control input $F_{\text{test}}=-0.5+0.2\sin{(t)}+0.3\sin{(2t)}$ for time $t=0$ to $t=2$.
 Fig.~\ref{Single_Pen} shows the resulting SINDy-PI models. The structure of the model is correctly identified up to a noise magnitude of $0.01$. Beyond this noise level, the SINDy-PI identified model only has short term prediction ability.
\subsection{Simplified Model of the Belousov–Zhabotinsky Reaction}
\label{secBZ}
\begin{figure}[t]
    \vspace{-.2in}
    \centering
    \includegraphics[width=1\textwidth]{Figures/BZ_Reaction.pdf}
    \caption{SINDy-PI is able to identify the simplified Belousov–Zhabotinsky reaction model.}
    \label{BZ_Reaction}
\end{figure}

We now apply SINDy-PI to a challenging PDE with rational nonlinearities, a simplified model of the Belousov\hyp{}Zhabotinsky~(BZ) reaction. The simplified BZ reaction model is given by~\cite{Vanag2004BZ_Reaction}
\begin{subequations}
    \label{s5.3.2eq1}
    \begin{align}
       \frac{\partial x}{\partial \tau}&=\frac{1}{\varepsilon}\left(\frac{fz(q-x)}{q+x}+x-x^2-\beta x+s\right)+\frac{D_x}{D_u}\Delta x,
       \label{s5.3.2eq1sub1}
       \\
       \frac{\partial z}{\partial \tau}&=x-z-\alpha z+\gamma u+\frac{D_z}{D_u}\Delta z,
       \label{s5.3.2eq1sub2}
       \\
       \frac{\partial s}{\partial \tau}&=\frac{1}{\varepsilon_2}(\beta x-s+\chi u)+\frac{D_s}{D_u}\Delta s,
       \label{s5.3.2eq1sub3}
       \\
       \frac{\partial u}{\partial \tau}&=\frac{1}{\varepsilon_3}[\alpha z -(\gamma +\frac{\chi}{2})u]+\frac{D_u}{D_u}\Delta u,
       \label{s5.3.2eq1sub4}
    \end{align}
\end{subequations}
where $x$, $z$, $s$, and $u$ are dimensionless variables and $\Delta=\frac{\partial^2}{\partial x_s^2}+\frac{\partial^2}{\partial y_s^2}$ denotes the Laplacian operator.

The strong coupling dynamics and implicit behavior in~\eqref{s5.3.2eq1sub1} make the data-driven discovery of the simplified BZ reaction challenging when using implicit-SINDy and PDE-FIND. However, SINDy-PI correctly identifies the simplified dynamics of the BZ-Reaction, as shown in Fig.~\ref{BZ_Reaction}.
To generate the simplified BZ reaction data, we use a spectral method~\cite{Trefethen_Book,kutz2013data} with time horizon $T=1$ and time step of $dt=0.001$.
We use $n=128$ discretization points with spatial domain ranging from $-10$ to $10$. The initial condition is chosen to be a mixture of Gaussian functions.
$80\%$ of the data is used for training, and the remaining $20\%$ is used for model selection.
The right-hand side library is normalized during the sparse regression process. A range of sparsity parameters $\lambda$ are tested from $0.1$ to $1$, with increments of $0.1$
The other system parameters in~\eqref{s5.3.2eq1} are given in App.~\ref{E1} and the SINDy-PI model is given in App.~\ref{App:SINDyBZ}.

\subsection{Extracting Physical Laws and Conserved Quantities}
\label{sec6}
In this final example, we demonstrate how to use SINDy-PI to extract governing physical laws and conserved quantities from data.
Many systems of interest are governed by Hamiltonian or Lagrangian dynamics.
Instead of identifying the ODE or PDE equations of motion, it might be possible to extract the physical laws directly.
These equations contain important information about the system and may be more concise, useful, and straightforward than the underlying ODE or PDE.
For example, given a Lagrangian, we can  derive the equations of motion.

The most difficult aspect of using SINDy-PI to identify a physical law is how to build the library. Conservation laws may contain higher-order derivatives, such as $\ddot{x}$. To include all possible terms, the library may become exceedingly large.
The library size will also increase if the system has many states.
Large libraries make the sparse regression sensitive to noise. Thus, extracting the physical law from data using SINDy-PI is still challenging due to the lack of constraints when constructing the library function. We only show one example in our paper to demonstrate that it is possible to achieve this using SINDy-PI, but further work is required to reduce the library size so that the sparse regression is robust.

\begin{figure}[t]
    \vspace{-.2in}
    \centering
    \includegraphics[width=1\textwidth]{Figures/Extract_Law.pdf}
    \caption{SINDy-PI is used to extract the conserved quantity for a double pendulum.}
    \label{ExtractLaw}
\end{figure}

As an example, we consider the double pendulum shown in Fig.~\ref{ExtractLaw}, with the system parameters given in App.~\ref{E1}. In this case, we also account for the friction in the pendulum joint, with friction constants of $k_1=7.2484\times 10^{-4}$ and $k_2=1.6522\times 10^{-4}$ for the pendulum arms, respectively.
In this case, we extract the Lagrangian of the double pendulum~\cite{Graichen2007} using SINDy-PI.
To extract this Lagrangian, we simulate the system with initial condition $x_{\text{train}}=[\pi-0.6,\ \pi-0.4,\ 0,\ 0]^T$ from $t=0$ to $t=15$  with time step $dt=0.001$. The resulting model is shown in Fig.~\ref{ExtractLaw}.

\section{Conclusions and Future Work}
\label{sec7}
In this paper, we develop SINDy-PI (parallel,implicit), a robust variant of the SINDy algorithm to identify implicit dynamics and rational nonlinearities.
SINDy-PI overcomes the sensitivity of the previous implicit-SINDy approach, which is based on a null-space calculation, making it highly sensitive to noise.
Instead, we introduce both parallel and constrained optimizations to test candidate terms in the dynamics, making the new SINDy-PI algorithm as robust as the original SINDy algorithm.
We also extend the algorithm to incorporate external forcing and actuation, making it more applicable to real-world systems.
We demonstrate this approach on several challenging systems with implicit and rational dynamics, including ODEs, actuated systems, and PDEs.
In particular, we discover the implicit dynamics for a simplified model for the BZ chemical reaction PDE, the double pendulum mechanical system, and the yeast glycolisis model, which have all been challenging test cases for advanced identification techniques.
Throughout these examples, we demonstrate considerable noise robustness and reductions to the data required, over the previous implicit-SINDy algorithm.

Despite the advances outlined here, there are still many important avenues of future work.
One limitation of this approach, and of SINDy in general, is in the design of the library of candidate functions.  The goal is a descriptive library, but the library size grows rapidly, which in turn makes the sparse regression ill-conditioned; other issues effecting robustness are discussed in App.~\ref{App:Robustness}. Recently, tensor approaches have been introduced to alleviate this issue, making libraries both descriptive and tractable~\cite{Gelss2019mindy}, and this is a promising approach that may be incorporated in SINDy-PI as well.
More generally, automatic library generation, guided by expert knowledge, is an important topic.
Other research directions will involve parameterizing elements of the library, so that the algorithm simultaneously identifies the model structure and the parameters of the sparsely selected terms.
Recent unified optimization frameworks, such as SR3~\cite{ZhengSR3,Champion2019SINDySR3}, may make this possible.
Model selection is another key area that will required focused attention.
Balancing accuracy on test data, sparsity of the model, and the potential for overfitting are all serious concerns.
The sparse regression and optimization may also be improved for better noise robustness.
Finally, modifying SINDy-PI to incorporate prior physical knowledge and to only model the discrepancy with an existing model~\cite{KK2019Discrepancy} will be the focus of ongoing work.

\section*{Acknowledgments}
SLB acknowledges support from the Army Research Office (ARO W911NF-19-1-0045) and the Air Force Office of Scientific Research (AFOSR FA9550-18-1-0200).
JNK acknowledges support from the Air Force Office of Scientific Research (AFOSR FA9550-17-1-0329).
We also acknowledge valuable discussions with Aditya Nair, Eurika Kaiser, Brian DeSilva, Tony Piaskowy, Jared Callaham, and Benjamin Herrmann. We thank Ariana Mendible for reviewing the manuscript and providing useful suggestions.

 \begin{spacing}{.8}
 \small{
 \setlength{\bibsep}{1.5pt}
\bibliographystyle{IEEEtran}

redpajama_set_name: RedPajamaArXiv