A Differentiable Point Process with Its Application to Spiking Neural Networks
June 02, 2021
Abstract
This paper is concerned about a learning algorithm for a probabilistic model of spiking neural networks (SNNs). [1] ([1]) proposed a stochastic variational inference algorithm to train SNNs with hidden neurons. The algorithm updates the variational distribution using the score function gradient estimator, whose high variance often impedes the whole learning algorithm. This paper presents an alternative gradient estimator for SNNs based on the path-wise gradient estimator. The main technical difficulty is a lack of a general method to differentiate a realization of an arbitrary point process, which is necessary to derive the path-wise gradient estimator. We develop a differentiable point process, which is the technical highlight of this paper, and apply it to derive the path-wise gradient estimator for SNNs. We investigate the effectiveness of our gradient estimator through numerical simulation.
1 Introduction↩︎
A spiking neural network (SNN) is an artificial neural network (ANN) where neurons communicate with each other using spikes rather than real values as the conventional ANNs do. The conventional ANN is a special case of SNN where information is encoded into the firing rate of neurons (which we call the rate coding) and the rate serves as the communication currency. This specification facilitates developing learning algorithms for ANNs, leading to the recent great success of deep neural networks. On the other hand, in the community of neuroscience, experimental evidence on biological neurons indicates that the rate coding alone cannot explain the whole brain activity [2] and more precise modeling of neural coding is anticipated. Since there still exist performance gaps between the rate-based ANNs and biological neural networks (i.e., brains) in terms of inference capability and energy efficiency, this raises the following question: how much of the current performance gaps can be attributed to this difference on neural coding? This open problem motivates us to study SNNs.
One of the major obstacles towards answering it is a lack of practical learning algorithms for SNNs, which discourages us from empirical investigation. While there exist a number of attempts to develop learning algorithms, most of them have more or less limited applicability. We consider a practical learning algorithm should at least be (i) theoretically grounded, (ii) empirically confirmed to work well, and (iii) easy to simulate (fewer hyperparameters, less computation time, etc.)1. For example, theoretical aspects of the algorithms based on spike-timing-dependent plasticity (Chapter 19 [3]) are not well understood. For another example, simulating learning algorithms for continuous-time deterministic SNNs requires the step-size parameter of time-axis discretization when the dynamics of a neuron is described by differential equations (e.g., [4]). The step-size parameter brings about the trade-off between the simulation quality and computation time, which makes the simulation more intricate. These examples illustrate that even major approaches do not satisfy all the requirements above, and therefore, there still exists much room for improvement.
Among a number of approaches, we employ as a foundation a probabilistic formulation of SNNs [5], which models spike trains (temporal sequence of spikes emitted from neurons) as a realization of a multivariate point process. It is easier for us to start from it than others because it already satisfies requirements (i) and (iii), which are more intrinsic properties than requirement (ii). In fact, learning algorithms are formalized by maximum likelihood estimation, and its exact simulation has no trade-off hyperparameter as will be explained in Section 2.2. Therefore, the remaining concern is its empirical performance.
One of the state-of-the-art learning algorithms for probabilistic SNNs is the work by [1] ([1]). The authors propose a stochastic variational inference algorithm for SNNs with hidden neurons. Since spike trains of hidden neurons are unobservable and it is intractable to compute the marginal likelihood, an evidence lowerbound (ELBO) is instead used as the objective function (Section 4.2). The key factor for optimizing ELBO is the way we estimate the gradient of ELBO. The authors employed the score function gradient estimator, also known as the REINFORCE estimator, which is widely applicable but is often reported to suffer from its high variance.
Our main idea is to substitute a path-wise gradient estimator for the score function gradient estimator. The path-wise gradient estimator tends to have lower variance than the score function gradient estimator [6], but it is not widely applicable (and is not applicable to SNNs) because it requires a sample from the variational distribution to be differentiable. Our contribution is that we develop a differentiable point process (Section 3) and apply it to derive the path-wise gradient estimator for SNNs (Section 4.2.2). \(\partial\mathcal{PP}\) is a point process whose realization is differentiable with respect to the model parameter,
We empirically investigate the effectiveness of the proposed learning algorithm in Section 5. We will confirm that (i) the proposed gradient estimator has lower variance than the existing one and (ii) this lower variance contributes to improve the performance of the learning algorithm. By comparing the performance of the proposed and existing ones, we obtain experimental results supporting these hypotheses. Therefore, we conclude that our path-wise gradient estimator improves empirical performance of SNNs.
One of the limitations of our learning algorithm as compared to the existing algorithm [1] is computation time. Since our algorithm generates more hidden spikes than the existing one does, our algorithm requires more computation time. We empirically examine the computational overhead of our algorithm against the existing one, and find that our algorithm requires 2.8 times more computation time than the existing one.
Notation. Let \([N]=\{1,2,\dots,N\}\). For any vector \(\mathbf{x}\), its \(d\)-th element is represented by \(x_d\). \(\begin{bmatrix}x_d\end{bmatrix}_{d\in[D]}\) denotes a \(D\)-dimensional vector whose \(d\)-th element is \(x_d\). For any vector \(\mathbf{x}\in\mathbb{R}^D\) and scalar \(c\in\mathbb{R}\), \(\begin{bmatrix}\mathbf{x}^\top & c\end{bmatrix}^\top\) denotes the \((D+1)\)-dimensional vector concatenating \(\mathbf{x}\) and \(c\). Let \(\mathbb{R}_{\geq 0}=\{x \geq 0 \}\) and \(\mathbb{R}_{> 0}=\{x > 0 \}\). Let \(\mathbb{1}^D=\{\mathbf{1}_d \}_{d\in[D]}\) be the set of \(D\)-dimensional one-hot vectors, where \(\mathbf{1}_d\in\{0,1\}^D\) is the one-hot vector whose \(d\)-th element is \(1\) and the others are \(0\). For any set \(A\), let \(\mathrm{conv}(A)\) be its convex hull, let \(\mathrm{conv}_{0}(A):=\mathrm{conv}(A\cup\{0\})\). Let \(\mathrm{Cat}(\mathbf{p})\) be the categorical distribution with parameter \(\mathbf{p}\in\mathrm{conv}(\mathbb{1}^D)\), whose random variable takes \(\mathbf{1}_d\in\mathbb{1}^D\) with probability \(p_d\). Let \(U[a,b]\) denote the uniform distribution over \([a,b]\). For any expectation operator \(\mathbb{E}_p\), let \(\hat{\mathbb{E}}_p\) be its Monte-Carlo approximation using an i.i.d. sample from \(p\).
2 Preliminaries↩︎
This section introduces temporal point processes along with their parameter estimation and sampling methods.
2.1 Point Processes↩︎
A point process [7] is a probabilistic model of an event collection. It is called a temporal point process when the event collection evolves in time. This paper only deals with a temporal point process, and therefore, we refer to it as a point process. We assume that point processes are simple, i.e., no events coincide.
2.1.1 Univariate Point Process↩︎
Assume we observe a sequence of \(N\in\mathbb{N}\) discrete events during time interval \([0, T]\), and let \(\mathcal{T}\) denote such an observation. \(\mathcal{T}\) can be represented by a series of event time stamps \(\{t_n\in[0,T]\}_{n\in [N]}\) as well as the information that we observe no event during \([0,t_1)\), \(\{(t_n, t_{n+1})\}_{n=1}^{N-1}\), and \((t_N, T]\). Let \(\mathcal{T}^{\leq t_n}\) represents a partial observation of \(\mathcal{T}\) up to and including time \(t_n\). One way of modeling \(\mathcal{T}\) is to specify the probability density function of the event time stamp \(t_{n+1}\) given the collection of its past events \(\mathcal{T}^{\leq t_n}\), which we describe, \(f(t\mid \mathcal{T}^{\leq t_n})\). Note that the probability density function must satisfy \(f(t \mid \mathcal{T}^{\leq t_n}) = 0\) for \(t\leq t_n\) and \(\int_{t_n}^\infty f(t \mid \mathcal{T}^{\leq t_n}) \dl{t} = 1\). The cumulative distribution function can be defined accordingly: \(F(t \mid \mathcal{T}^{\leq t_n}) = \int_{t_n}^{t} f(s \mid \mathcal{T}^{\leq t_n}) \dl{s} = \Pr[t_{n+1}\in(t_n, t)\mid \mathcal{T}^{\leq t_n}]\).
Another way of modeling it is to specify the conditional intensity function, which is related to the distributions as, \[\begin{align} \label{eq:5} \lambda(t \mid \mathcal{T}^{\leq t_n}) = \begin{cases} \dfrac{f\left(t \mid \mathcal{T}^{\leq t_n}\right)}{1 - F\left(t \mid \mathcal{T}^{\leq t_n}\right)} & (t > t_n),\\ 0 & (t \leq t_n). \end{cases} \end{align}\tag{1}\] In the following, let \(t_n\) denote an arbitrary event time stamp and we only specify the conditional intensity function for \(t>t_n\), because its value for \(t\leq t_n\) is trivially \(0\). Observing that \(\lambda(t \mid \mathcal{T}^{\leq t_n}) \dl{t} = \Pr[t_{n+1}\in[t, t+\dl{t}] \mid t_{n+1}\notin (t_n, t), \mathcal{T}^{\leq t_n}]\) holds as \(\dl{t}\rightarrow +0\) [8], the conditional intensity function represents how likely the event occurs at time \(t\) given that we have observed \(n\) events so far and no event has been observed during \((t_n, t)\).
A point process is more often specified by the conditional intensity function than the time interval distribution. Let \(\mathcal{PP}(\lambda)\) be the point process with the conditional intensity function \(\lambda\). Corollary 1, which is an immediate consequence of Proposition 2, states the conditions under which the conditional intensity function uniquely specifies a point process.
Corollary 1. A conditional intensity function \(\lambda\) uniquely defines a point process if it satisfies the following conditions for any observation of discrete events \(\mathcal{T}^{\leq t_n}\) and any \(t>t_n\):
\(\lambda(t \mid \mathcal{T}^{\leq t_n})\) is non-negative and integrable on any interval starting at \(t_n\),
\(\int_{t_n}^t \lambda(s\mid \mathcal{T}^{\leq t_n})\dl{s}\rightarrow\infty\) as \(t\rightarrow\infty\), and
\(\int_{t_n}^t \lambda(s\mid \mathcal{T}^{\leq t_n})\dl{s}\) is right continuous w.r.t. \(t\).
The log-likelihood of observation \(\mathcal{T}\) on \(\mathcal{PP}(\lambda)\) is given as, \[\begin{align} \label{eq:7} \log p(\mathcal{T}) = \sum_{t\in\mathcal{T}} \log \lambda\left(t\mid \mathcal{T}^{\leq t_{n(t)}}\right) - \Lambda^{[0,T]}(\mathcal{T}), \end{align}\tag{2}\] where let \(\Lambda^{[0,T]}(\mathcal{T})=\int_0^T \lambda(t\mid \mathcal{T}^{\leq t_{n(t)}})\dl{t}\) be the integrated conditional intensity function, also known as the compensator, which accounts for no-event periods, and let \(n(t):\mathbb{R}_{\geq 0}\rightarrow\mathbb{Z}_{\geq 0}\) be the left-continuous2 counting process of the observation \(\mathcal{T}\), which counts the number of events up to but not including time \(t\). The latest event time stamp at time \(t\) can be denoted by \(t_{n(t)}\in[0,t)\). Figure 1 illustrates a realization of a point process and its counting representation. A typical procedure of modeling \(\mathcal{T}\) is to design a parametric model of the conditional intensity function that satisfies the conditions of Corollary 1 and train it by maximizing the log-likelihood function (Equation 2 ).
2.1.2 Multivariate Point Process↩︎
A multivariate point process is a set of mutually dependent point processes and can be defined via a marked point process, in which each event is associated with a mark. We call a marked point process whose mark belongs to set \(X\), an \(X\)-marked point process. Let \(\mathcal{T}_X\) denote an observation of an \(X\)-marked point process, which contains a series of event time stamps and marks, \(\{(t_n, \mathbf{p}_n)\in[0, T]\times X\}_{n\in[N]}\). As illustrated in Figure 2, a \(D\)-variate point process can be defined by a \(\mathbb{1}^D\)-marked point process, where each mark \(\mathbf{p}_n\) indicates which dimension the event belongs to. For example, if \(\mathbf{p}_n=\mathbf{1}_1\), the \(n\)-th event occurs at the first dimension. In Figure 2, blue-circle and red-diamond marks correspond to the first and the second dimensions respectively.
Letting \(f(t, \mathbf{p}\mid \mathcal{T}_{\mathbb{1}^D}^{\leq t_{n}})\) be the probability density function of each event \((t_{n+1},\mathbf{p}_{n+1})\) given its past events \(\mathcal{T}_{\mathbb{1}^D}^{\leq t_n}\), the conditional intensity function can be defined similarly: \[\begin{align} \label{eq:8}\lambda\left(t, \mathbf{p}\mid \mathcal{T}_{\mathbb{1}^D}^{\leq t_n}\right) = \dfrac{f\left(t, \mathbf{p}\mid \mathcal{T}_{\mathbb{1}^D}^{\leq t_n}\right)}{1 - F\left(t \mid \mathcal{T}_{\mathbb{1}^D}^{\leq t_n}\right)}, \end{align}\tag{3}\] where \(F(t \mid \mathcal{T}_{\mathbb{1}^D}^{\leq t_{n}}) = \int_{t_{n}}^t \dl{s} \sum_{\mathbf{p}\in \mathbb{1}^D} f(s, \mathbf{p}\mid \mathcal{T}_{\mathbb{1}^D}^{\leq t_{n}})\). The conditional intensity function represents how likely event \((t,\mathbf{1}_d)\) occurs: \(\lambda(t,\mathbf{1}_d \mid \mathcal{T}_{\mathbb{1}^D}^{\leq t_n}) \dl{t} = \Pr[t_{n+1}\in[t, t+\dl{t}], \mathbf{p}_{n+1}=\mathbf{1}_d \mid t_{n+1}\notin (t_n, t), \mathcal{T}_{\mathbb{1}^D}^{\leq t_n}]\). Proposition 2 states conditions under which the conditional intensity function uniquely specifies a marked point process. See Appendix 9 for its proof. Let \(\mathcal{MPP}(\lambda)\) be the multivariate point process with the conditional intensity function \(\lambda\).
Proposition 2. Let \(X\) be a set. A conditional intensity function \(\lambda\) uniquely defines an \(X\)-marked point process if it satisfies the following conditions for any \(\mathcal{T}_X^{\leq t_n}\) and \(t>t_n\):
\(\lambda(t,\mathbf{p}\mid\mathcal{T}_{X}^{\leq t_n})\geq 0\) and integrable w.r.t. \(\mathbf{p}\) and w.r.t. \(t\) on any interval starting at \(t_n\),
\(\int_{t_n}^t \dl{s} \int_{X}\dl{\mathbf{p}} \lambda(s, \mathbf{p}\mid \mathcal{T}_X^{\leq t_n}) \rightarrow\infty\) as \(t\rightarrow\infty\), and
\(\int_{t_n}^t \dl{s} \int_{X}\dl{\mathbf{p}} \lambda(s, \mathbf{p}\mid \mathcal{T}_X^{\leq t_n})\) is right continuous in \(t\).
The log-likelihood of observation \(\mathcal{T}_{\mathbb{1}^D}\) is written as: \[\begin{align} \label{eq:12} \begin{aligned} & \log p(\mathcal{T}_{\mathbb{1}^D}) \\ = &\sum_{(t,\mathbf{p})\in\mathcal{T}_{\mathbb{1}^D}} \log \lambda\left(t,\mathbf{p}\mid \mathcal{T}_{\mathbb{1}^D}^{\leq t_{n(t)}}\right) - \Lambda^{[0,T]}(\mathcal{T}_{\mathbb{1}^D}), \end{aligned} \end{align}\tag{4}\] where let \(\Lambda^{[0,T]}(\mathcal{T}_{\mathbb{1}^D})=\int_{0}^T \sum_{\mathbf{p}\in\mathbb{1}^D}\lambda(t, \mathbf{p}\mid \mathcal{T}_{\mathbb{1}^D}^{\leq t_n(t)}) \dl{t}\) be the compensator. Since its analytical form is not available for a general conditional intensity function, we resort to Monte-Carlo approximation to estimate the compensator. In specific, we draw \(M\) examples, \(\{t_m\}_{m\in[M]}\), from \(U[0,T]\) and approximate it as, \[\begin{align} \label{eq:1}\Lambda^{[0,T]}(\mathcal{T}_{\mathbb{1}^D})\approx &\frac{T}{M} \sum_{m=1}^M \sum_{\mathbf{p}\in\mathbb{1}^D} \lambda\left(t_m, \mathbf{p}\mid \mathcal{T}_{\mathbb{1}^D}^{\leq t_{n(t_m)}}\right). \end{align}\tag{5}\]
2.2 Sampling Algorithms↩︎
This section introduces sampling algorithms for a point process given a conditional intensity function. A notable feature of the algorithms is that they can exactly simulate point processes without any approximation. This indicates that there exists no hyperparameter controling the trade-off between computational cost and accuracy of the simulation, which greatly facilitates simulating SNNs.
2.2.1 Homogeneous Poisson Process↩︎
The simplest point process is the homogeneous Poisson process whose conditional intensity function is constant; \(\lambda(t\mid\mathcal{T}^{\leq t_n})=\lambda\) for any \(\mathcal{T}^{\leq t_n}\). It is straightforward to sample from it because the interval between two successive events \(\tau\) is independently and identically distributed according to the exponential distribution, \(f(\tau; \lambda)=\lambda \mathrm{e}^{-\lambda \tau}\).
2.2.2 General Point Process↩︎
It is not straightforward to sample from a general point process when a closed-form expression of the inter-event time distribution is not available. This is true for many point processes including SNNs. Among several sampling methods, the thinning algorithm [9], [10] allows us to sample from such a point process without knowing the closed-form expression. For other sampling algorithms, please refer to Section 6.
The main idea is to generate a sequence of time stamps from a homogeneous Poisson process with sufficiently high intensity (which we call the base process) and then to thin some of the events so that the sequence follows the given point process. Algorithm 3 describes it for the multivariate case, where let \(\boldsymbol{\lambda}\left(t \mid \mathcal{T}\right) = \begin{bmatrix}\lambda\left(t, \mathbf{1}_d \mid \mathcal{T}\right) \end{bmatrix}_{d\in[D]}\), and let \(\boldsymbol{\pi}_{\bar{\lambda}}\circ\) be an operator that receives a \(D\)-dimensional vector \(\boldsymbol{\lambda}\) and returns \(\frac{1}{\bar{\lambda}}\begin{bmatrix}\boldsymbol{\lambda} \\ \bar{\lambda} - \left\|\boldsymbol{\lambda}\right\|_1 \end{bmatrix}\).
It first generates a new time stamp \(s\) from the homogeneous Poisson process with intensity \(\bar{\lambda}\) (line 3). Then it decides whether or not to accept the event, and if accepting, decides which dimension the event is assigned to (lines 6-8); \(s\) is rejected if \(r=1\), i.e., with probability \(1-{\frac{1}{\bar{\lambda}}}{\sum_{\mathbf{p}\in \mathbb{1}^D}\lambda(s, \mathbf{p}\mid \mathcal{T})}\), and \(s\) is accepted as the event from the \(d\)-th dimension (\(d\in[D]\)) if \(p_d=1\), i.e., with probability \(\lambda\left(s, \mathbf{1}_d \mid \mathcal{T}\right)/\bar{\lambda}\),
Intuitively, the correctness of Algorithm 3 is understood as follows. Assuming we have sampled \(\mathcal{T}_{\mathbb{1}^D}^{\leq t_n}\), at any time \(t > t_n\), the probability that the algorithm generates the event with mark \(\mathbf{1}_d\) in interval \([t, t+\dl{t}]\) is, \[\begin{align} &\Pr\left[t_{n+1}\in[t, t+\dl{t}], \mathbf{p}_{n+1}=\mathbf{1}_d \mid \mathcal{T}_{\mathbb{1}^D}^{<t}\right]\\ =& \underbrace{\bar{\lambda} \dl{t}}_{\substack{\text{Prob. that the base process}\\ \text{generates the event in } [t, t+\dl{t}]}} \cdot \underbrace{\lambda\left(t, \mathbf{1}_d \mid \mathcal{T}_{\mathbb{1}^D}^{\leq t_n}\right)/\bar{\lambda}}_{\substack{\text{Prob. that }t\text{ is}\\\text{assigned the }d\text{-th mark}}}\\ =& \lambda\left(t, \mathbf{1}_d \mid \mathcal{T}_{\mathbb{1}^D}^{\leq t_n}\right) \dl{t}, \end{align}\] where let \(\mathcal{T}_{\mathbb{1}^D}^{<t}\) denote the event \(t_{n+1}\notin (t_n,t)\) and \(\mathcal{T}_{\mathbb{1}^D}^{\leq t_n}\). This shows that the output follows \(\mathcal{MPP}(\lambda)\). For its formal proof, please refer to Reference [10].
3 Differentiable Point Process↩︎
We present the key building block of our method called a differentiable point process, whose realization is differentiable with respect to its parameters. Differentiability plays an essential role when designing a learning algorithm for latent variable models as will be discussed in Section 4.2.
The key idea is that the output of Algorithm 3 becomes differentiable if we replace the categorical distribution in line 6 with a reparameterizable distribution such as the concrete distribution, also known as the Gumbel-softmax distribution [11], [12]. We first review the concrete distribution (Section 3.1), and then we present the differentiable point process (Section 3.2).
3.1 Concrete Distribution↩︎
The concrete distribution has been developed as a reparameterizable substitute for the categorical distribution. The idea comes from the Gumbel-max trick, which enables us to sample from the categorical distribution. Letting \(\boldsymbol{\pi}\in\mathbb{R}_{\geq 0}^D\) be an unnormalized parameter of the categorical distribution, the Gumbel-max trick first samples \(u_d\sim U[0,1]\) for each \(d\in[D]\), and then outputs \(\mathbf{1}_{d^\star}\) where \(d^\star=\arg\max_{d\in[D]} \log \pi_d - \log(-\log u_d)\). The output is known to be distributed according to \(\mathrm{Cat}(\boldsymbol{\pi}/\|\boldsymbol{\pi}\|_1)\). While the Gumbel-max trick successfully divides the sampling procedure into random sampling from the fixed distribution and a parameterized transformation of it, which is necessary to be differentiable, the gradient of its realization with respect to \(\boldsymbol{\pi}\) is non-informative, because a small variation to \(\boldsymbol{\pi}\) does not change the gradient.
The concrete distribution is defined by relaxing the range of the random variable from \(\mathbb{1}^D\) to its convex hull \(\mathrm{conv}(\mathbb{1}^D)\) so that its gradient is more informative. Accordingly, the argmax operator in the Gumbel-max trick is replaced with the softmax operator with temperature \(\tau>0\). Since softmax becomes equivalent to argmax as \(\tau\rightarrow 0\), the concrete distribution also becomes equivalent to the categorical distribution as \(\tau\rightarrow 0\). Let \(g_\tau(\mathbf{p};\boldsymbol{\pi})\) denote the probability density function of the concrete distribution with temperature \(\tau\) and unnormalized parameter \(\boldsymbol{\pi}\in\mathbb{R}^{D}_{>0}\).
3.2 Multivariate Differentiable Point Process↩︎
We present a constructive definition of a differentiable point process in Definition 3.
Definition 3. Assume the conditional intensity function \(\lambda(t, \mathbf{p}\mid \mathcal{T}_{\mathbb{1}^D}^{\leq t_n})\) can be computed with an observation of a \(\mathrm{conv}_{0}(\mathbb{1}^D)\)-marked point process. Let \(\bar{\lambda}\) be a constant satisfying \(\bar{\lambda}>\sum_{\mathbf{p}\in\mathbb{1}^D} \lambda(t, \mathbf{p}\mid \mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{\leq t_n})\) for any \(\mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{\leq t_n}\) and \(t > t_n\), and \(\tau > 0\) be temperature. The differentiable point process \(\mathcal{\partial PP}\left(\lambda, \bar{\lambda}, \tau\right)\) is defined as a \(\mathrm{conv}_{0}(\mathbb{1}^D)\)-marked point process constructed by Algorithm 4.
The main idea is to substitute the concrete distribution for the categorical distribution in line 5 of Algorithm 3, resulting in Algorithm 4. In Algorithm 3, \(r_k\) indicates whether \(s_k\) should be discarded or not, and if not, mark \(\mathbf{p}_k\in\mathbb{1}^D\) is assigned to the event time stamp \(s_k\). In Algorithm 4, each event is always accepted, and \(\mathbf{p}_k\in\mathrm{conv}_{0}(\mathbb{1}^{D})\) is used as a mark of the event \(s_k\).
Algorithms 3 and 4 are different in two ways. First, all events from the base process are accepted in Algorithm 4, while some are rejected in Algorithm 3. Second, in Algorithm 3, the mark is defined over \(\mathbb{1}^D\), while in Algorithm 4, it is defined over \(\mathrm{conv}_{0}(\mathbb{1}^D)\); each mark is associated with amplitude that is continuous w.r.t. the model parameter.
The differentiable point process as defined above can be understood as a marked point process (Proposition 4).
Proposition 4. The differentiable point process \(\mathcal{\partial PP}(\lambda,\bar{\lambda},\tau)\) is a \(\mathrm{conv}_{0}(\mathbb{1}^D)\)-marked point process with conditional intensity function, \[\begin{align} &\lambda_{\partial}\left(t, \mathbf{p}\mid \mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^{D})}^{\leq t_n};\boldsymbol{\lambda}, \bar{\lambda}, \tau\right) \\ = & \bar{\lambda} \cdot g_\tau\left( \begin{bmatrix} \mathbf{p}\\ 1- \|\mathbf{p}\|_1 \end{bmatrix}; \boldsymbol{\pi}_{\bar{\lambda}}\circ\boldsymbol{\lambda}\left(t \mid \mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{\leq t_n}\right)\right). \end{align}\] where let, \[\begin{align} \boldsymbol{\lambda}\left(t \mid \mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{\leq t_n}\right) &= \begin{bmatrix}\lambda\left(t, \mathbf{1}_d \mid \mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{\leq t_n}\right) \end{bmatrix}_{d\in[D]},\\ \boldsymbol{\pi}_{\bar{\lambda}}\circ\boldsymbol{\lambda} &= \frac{1}{\bar{\lambda}}\begin{bmatrix}\boldsymbol{\lambda} \\ \bar{\lambda} - \left\|\boldsymbol{\lambda}\right\|_1 \end{bmatrix}. \end{align}\]
We can confirm the differentiability of a realization of \(\partial\mathcal{PP}\) (Proposition 8). We can also confirm that in the limit of \(\tau\rightarrow 0\), the differentiable point process becomes equivalent to the original point process (Proposition 9). See Appendix 10 for their formal statements and proofs.
As discussed by [11] ([11]), the concrete distribution often suffers from underflow and we have to implement it in the logarithmic scale. Our implementation also suffers from the same issue, and we provide a numerically stable implementation idea in Appendix 12.
4 Learning Algorithm for SNNs↩︎
We present a learning algorithm for spiking neural networks (SNNs) based on the differentiable point process. We first define a probabilistic model of SNNs (Section 4.1) and then will present our learning algorithm, highlighting the difference from the existing one (Section 4.2).
4.1 Probabilistic Model of Spiking Neural Networks↩︎
We employ the standard probabilistic model in the literature [5]. Let \(D\) be the number of neurons, let \(\mathcal{N}=\mathbb{1}^D\) be the set of neurons, each of which is indexed by a one-hot vector, and let \(\mathcal{T}_{\mathcal{N}}\) be spike trains emitted from SNN during time interval \([0, T]\). We assume that \(\mathcal{T}_{\mathcal{N}}\) is a realization of an \(\mathcal{N}\)-marked point process. We define the conditional intensity function based on a spike response model (SRM) [3]. SRM assumes that the \(d\)-th spiking neuron is driven by its internal state called a membrane potential, \[\begin{align} \label{eq:11} \begin{aligned} u_d\left(t \mid \mathcal{T}_{\mathcal{N}}^{\leq t_n}\right) = \bar{u}_d + \sum_{(t^\prime, \mathbf{p})\in\mathcal{T}_{\mathcal{N}}^{\leq t_n}} \mathbf{f}_{d}(t-t^\prime)\cdot \mathbf{p}, \end{aligned} \end{align}\tag{6}\] where \(\mathbf{f}_d(s) = \begin{bmatrix}f_{d^\prime,d}(s) \end{bmatrix}_{d^\prime\in[D]}\) is a vector of filter functions from all of the neurons to the \(d\)-th neuron. In specific, \(f_{d^\prime, d}(s)\) describes the time course of the membrane potential of neuron \(d\) in response to a spike emitted by neuron \(d^\prime\) at time \(s=0\). We assume \(f_{d, d}(s)\leq 0\) for all \(d\in\mathcal{N}\). This assumption allows us to reproduce the resetting behavior of a biological neuron; the membrane potential is reset to a lower level after the neuron fires. We also assume that \(f_{d^\prime, d}(s)=0\) for \(s<0\). This assumption ensures that future events have no influence on past events.
Then, the conditional intensity function is defined by, \[\begin{align} \label{eq:10}\lambda^{\mathrm{SNN}}(t, \mathbf{p}\mid \mathcal{T}_{\mathcal{N}}^{\leq t_n}) = \mathbf{p}\cdot \sigma(\mathbf{u}(t \mid \mathcal{T}_{\mathcal{N}}^{\leq t_n})), \end{align}\tag{7}\] where \(\sigma\colon\mathbb{R}^D\rightarrow\mathbb{R}_{\geq 0}^D\) is element-wisely non-decreasing and differentiable3 and let \(\mathbf{u}(t \mid \mathcal{T}_{\mathcal{N}}^{\leq t_n}) = \begin{bmatrix} u_d(t \mid \mathcal{T}_{\mathcal{N}}^{\leq t_n}) \end{bmatrix}_{d\in[D]}\). As the membrane potential of one neuron increases, the neuron is more likely to fire and generate a spike.
For numerical simulation, we assume that the filter functions are parameterized by weights \(\{w_{d^\prime,d,l}\in\mathbb{R}\}_{l=1}^L\) as, \[\begin{align} \label{eq:21}f_{d^\prime,d}(s) = \begin{cases} \sum_{l=1}^L w_{d^\prime, d, l} \cdot \kappa(s-s_l) & (s\geq 0),\\ 0 & (s < 0), \end{cases} \end{align}\tag{8}\] where \(\{s_l\in\mathbb{R}\}_{l=1}^L\) are fixed and \(\kappa(s)=\max\{\frac{3}{4}(1-s^2), 0\}\) is the Epanechnikov kernel. We chose this kernel because the bounded support of the kernel allows us to ignore events that occurred more than a certain period ago for membrane potential computation. Let \(\theta=\{\bar{u}_d\in\mathbb{R}\}_{d=1}^D\cup\{w_{d^\prime, d, l}\in\mathbb{R}\mid l\in[L]\}_{d, d^\prime=1}^D\) denote the set of model parameters.
4.2 Learning Algorithms↩︎
Assume some of the neurons are hidden and their spike trains are unobservable. Let \(\mathcal{O}\subset\mathcal{N}\) and \(\mathcal{H}=\mathcal{N}\backslash\mathcal{O}\) be the sets of observable and hidden neurons, respectively. Accordingly, the spike trains of all of the neurons are divided into observable and hidden ones: \(\mathcal{T}_{\mathcal{N}}=\mathcal{T}_{\mathcal{O}}\cup\mathcal{T}_{\mathcal{H}}\). We consider an estimation procedure for the model parameters of SNN, \(\theta\), given a set of observed spike trains \(\left\{\mathcal{T}_{\mathcal{O}, n}\right\}_{n=1}^N\).
Letting \(p(\mathcal{T}_\mathcal{N};\theta)=p(\mathcal{T}_{\mathcal{O}}, \mathcal{T}_{\mathcal{H}};\theta)\) be the joint distribution of the observable and hidden spike trains, the parameter \(\theta\) is estimated by maximum likelihood estimation:
|l|[3] _n=1^N (; _,n)
where \(\ell(\theta;\mathcal{T}_{\mathcal{O}}) = \log \int p(\mathcal{T}_{\mathcal{O}}, \mathcal{T}_{\mathcal{H}}; \theta) \dl \mathcal{T}_\mathcal{H}\) is the marginalized log-likelihood function. Since it is intractable to compute it, we substitute its lower bound called an evidence lower bound (ELBO) for the marginalized log-likelihood function as the objective function: \[\begin{align} \nonumber \underline{\ell}(\theta,\phi;\mathcal{T}_{\mathcal{O}}) &= \mathbb{E}_{q(\mathcal{T}_{\mathcal{H}};\phi)}\left[\log p(\mathcal{T}_{\mathcal{O}}, \mathcal{T}_{\mathcal{H}}; \theta) - \log q(\mathcal{T}_{\mathcal{H}};\phi)\right],\\ \label{eq:14} &\equiv \mathbb{E}_{q(\mathcal{T}_{\mathcal{H}};\phi)} \underline{\ell}(\theta, \phi;\mathcal{T}_{\mathcal{O}}, \mathcal{T}_{\mathcal{H}}), \end{align}\tag{9}\] where \(q(\mathcal{T}_{\mathcal{H}};\phi)\) is an arbitrary distribution called a variational distribution, parameterized by \(\phi\). We specifically assume that the variational distribution is modeled by SNN driven by both observable and hidden spike trains. In the following, we omit the index of data \(n\) for ease of presentation and consider ELBO using a single observation \(\mathcal{T}_{\mathcal{O}}\).
Since there exists no closed-form solution to the maximization problem, we resort to stochastic gradient ascent methods, resulting in Algorithm 5. The basic procedure to train SNN is to choose one realization \(\mathcal{T}_{\mathcal{O}}\) from the data set randomly, and update \(\theta\) and \(\phi\) so as to maximize Equation 9 . In the following, we present both an existing approach and our novel approach to compute the gradients, \(\frac{\partial\underline{\ell}}{\partial\theta}\) and \(\frac{\partial\underline{\ell}}{\partial\phi}\).
4.2.1 Gradient with respect to \(\theta\)↩︎
The gradient with respect to \(\theta\) is straightforwardly computed by applying Monte-Carlo approximation: \(\diffp*{\underline{\ell}(\theta,\phi;\mathcal{T}_\mathcal{O})}{\theta} \approx \hat{\mathbb{E}}_{q(\mathcal{T}_{\mathcal{H}};\phi)}\left[\diffp*{\log p(\mathcal{T}_{\mathcal{O}}, \mathcal{T}_{\mathcal{H}}; \theta)}{\theta}\right]\). This can be numerically calculated with the help of automatic differentiation tools. Letting \(\mathcal{T}_{\mathcal{H},m}^{<T}\) (\(m\in[M]\)) be i.i.d. sample of size \(M\) from \(q(\mathcal{T}_{\mathcal{H}}^{<T};\phi)\), the gradient is approximated as,
4.2.2 Gradient with respect to \(\phi\)↩︎
The gradient with respect to \(\phi\) is more involved. In Equation 9 , the expectation operator depends on \(\phi\) and we cannot exchange \(\frac{\partial}{\partial\phi}\) and \(\mathbb{E}_{q(\mathcal{T}_\mathcal{H};\phi)}\). There are at least two approaches to computing the gradient in this situation [6]. One approach is to rely on the score function gradient estimator, also known as the REINFORCE estimator [13]. While it is widely applicable to a variety of models, it is often reported that the gradient estimator has high variance. Another approach is the path-wise gradient estimator, which makes use of the reparameterization trick [14]. While its variance is often reported to be lower than that of the score function gradient estimator [6], its application is limited because the probability distribution \(q\) must be reparameterizable.
In the literature of SNNs, the score function gradient estimator with respect to \(\phi\) has been developed by [1] ([1]). Our contribution is to develop a path-wise gradient estimator for SNNs based on a differentiable point process presented in Section 3.
Score function gradient estimator. [1] ([1]) used the score function gradient estimator for computing the gradient with respect to \(\phi\): \(\frac{\partial\underline{\ell}}{\partial\phi}(\theta, \phi;\mathcal{T}_\mathcal{O}) \approx\hat{\mathbb{E}}_{q(\mathcal{T}_\mathcal{H};\phi)}[\frac{\partial\log q(\mathcal{T}_\mathcal{H};\phi)}{\partial\phi} (\underline{\ell}(\theta, \phi;\mathcal{T}_{\mathcal{O}}, \mathcal{T}_{\mathcal{H}}) - 1)]\). While this is an unbiased estimator of the gradient, its high variance is often problematic. We employ the variational distribution with the conditional intensity function, \[\begin{align} \label{eq:19}\lambda_q(t, \mathbf{p}\mid \mathcal{T}_{\mathcal{N}}^{\leq t_n};\phi) = \mathbf{p}\cdot \sigma(\mathbf{u}(t \mid \mathcal{T}_{\mathcal{N}}^{\leq t_n};\phi)), \end{align}\tag{10}\] for any \(\mathbf{p}\in\mathcal{H}\). In particular, we use shared parameters for the model and the variational distribution, i.e., we set \(\phi=\theta\) as we observe it improves the performance.
Path-wise gradient estimator. We propose a path-wise gradient estimator for SNNs. Our main idea is to employ the differentiable point process, \(\mathcal{\partial PP}(\lambda_q(t, \mathbf{p}\mid \mathcal{T}_\mathcal{N};\phi); \bar{\lambda}, \tau)\), as the variational distribution, where \(\lambda_q\) is defined in Equation 10 . This allows us to differentiate a Monte-Carlo approximation of ELBO (Equation 9 ) using automatic differentiation tools: \[\begin{align} \label{eq:20}\diffp{\underline{\ell}(\theta,\phi;\mathcal{T}_{\mathcal{O}})}{\phi} \approx \diffp{\hat{\mathbb{E}}_{\partial\mathcal{PP}} \underline{\ell}(\theta,\phi;\mathcal{T}_{\mathcal{O}}, \mathcal{T}_{\mathrm{conv}_{0}(\mathcal{H})}(\phi))}{\phi}. \end{align}\tag{11}\]
The main technical issue in applying the differentiable point process is that its realization \(\mathcal{T}_{\mathrm{conv}_{0}(\mathcal{H})}(\phi)\) is incompatible with the SNN model defined by Equations 6 and 7 . The model assumes that a mark \(\mathbf{p}\) is a one-hot vector, while a mark of a differentiable point process belongs to \(\mathrm{conv}_{0}(\mathcal{H})\). We address this by devising a differentiable spiking neural network (\(\partial\)SNN), which can handle a mark in \(\mathrm{conv}_{0}(\mathcal{H})\), while keeping the conditional intensity function proper.
Let \(\bar{\mathcal{N}}=\mathcal{O}\cup\mathrm{conv}_{0}(\mathcal{H})\) be the set of marks for \(\partial\)SNN. We define the membrane potential of neuron \(d\in\mathcal{N}\) as, \[\begin{align} \label{eq:18} & u_d\left(t \mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_n}\right) = \bar{u}_d + \sum_{(t^\prime, \mathbf{p})\in\mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_n}} \mathbf{f}_{d}(t-t^\prime)\cdot \mathbf{p}, \end{align}\tag{12}\] and the conditional intensity of \(\partial\)SNN for \(\mathbf{p}\in\bar{\mathcal{N}}\) as, \[\begin{align} \label{eq:17}&\lambda^{\partial\mathrm{SNN}}\left(t, \mathbf{p}\mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_n}; \bar{\lambda}, \tau\right)\\ \nonumber =& \sum_{\mathbf{1}_d\in\mathcal{O}}\delta(\mathbf{p}-\mathbf{1}_d) \lambda^{\mathrm{SNN}}\left(t, \mathbf{p}\mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_n}\right)\\ \nonumber &+ \mathbb{I}[\mathbf{p}\in\mathrm{conv}_{0}(\mathcal{H})] \lambda_{\partial}\left(t, \mathbf{p}_\mathcal{H}\mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_n};\boldsymbol{\lambda}_{\mathcal{H}}, \bar{\lambda}, \tau\right) \end{align}\tag{13}\] where \(\boldsymbol{\lambda}_{\mathcal{H}}\left(t \mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_n}\right) = \sigma\left(\begin{bmatrix}u_{d}(t \mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_n})\end{bmatrix}_{d\in\mathcal{H}}\right)\), \(\mathbb{I}[\cdot]\) is the indicator function, and \(\mathbf{p}_\mathcal{H}= \begin{bmatrix}p_d \end{bmatrix}_{d\in\mathcal{H}}\).
It is necessary to confirm that (i) the conditional intensity function can be calculated using past events whose marks are in \(\bar{\mathcal{N}}\) and (ii) the conditional intensity function satisfies all of the conditions listed in Proposition 2 for \(\mathcal{X}=\bar{\mathcal{N}}\). The first requirement immediately follows from Equations 12 and 13 . In Appendix 11, we provide the formal statement and proof of the second requirement (Proposition 5). We also confirm that ELBO is differentiable (Proposition 6) and that the differentiable SNN becomes equivalent to the vanilla SNN in the limit of \(\tau\rightarrow 0\) (Proposition 7).
Proposition 5. Assume that the filter functions are parameterized as Eq. 8 . Then, the conditional intensity function (Eq. 13 ) uniquely defines an \(\bar{\mathcal{N}}\)-marked point process.
Finally, let us discuss two properties of the proposed method. Proposition 6 states that the objective function is differentiable with respect to \(\phi\). Proposition 7 states that the extended SNN model converges to the vanilla SNN as \(\tau\rightarrow+0\).
Proposition 6. The Monte-Carlo approximation of ELBO (Eq. 11 ) is differentiable with respect to \(\phi\) if we employ \(\mathcal{\partial PP}(\lambda_q(t, \mathbf{p}\mid \mathcal{T}_\mathcal{O}\cup\mathcal{T}_\mathcal{H};\phi); \bar{\lambda}, \tau)\) as the variational distribution.
Proposition 7. In the limit of \(\tau\rightarrow +0\), a realization of the extended SNN (Eq. 13 ) is distributed according to the vanilla SNN (Eq. 7 ) if we discard events with mark \(\mathbf{p}=\mathbf{0}\).
To summarize our contribution, we develop a path-wise gradient estimator for SNNs, whose variance is often reported to be lower than the existing estimator. We develop a differentiable point process, whose realization is differentiable (Section 3) and apply it to SNNs to derive the gradient estimator.
5 Empirical Studies↩︎
| Network size | \(D=6\), \(|\cO|=2\), \(|\cH|=4\) |
| Activation/filter functions | \(a=5\), \(L=2\), \(s_1=0\), \(s_2=10\) |
| \(\dpp\) | \(\tau=0.3\), \(\bar{\lambda}=20\) |
| # of samplings | \(100\) (Eq. [eq:1]), \(1\) (Eq. [eq:14]) |
Let us investigate the effectiveness of our gradient estimator through numerical simulation. Our hypothesis is that (i) the path-wise gradient estimator will have lower variance than the score function estimator and (ii) lower variance will improve the predictive performance. We design two experiments (Sections 5.1 and 5.2) to verify these two hypotheses. We additionally compare computation cost of the learning algorithms using each of the gradient estimators in Section 5.3. All the experiments are conducted on IBM Cloud4, and the code is publicly available [15].
Data set. We use a synthetic data set generated by the vanilla SNN (Equation 7 ). Table 1 summarizes its configuration. We set \(\bar{\lambda}=a|\mathcal{H}|=20\), which is the tightest upperbound because we use the sigmoid activation function with amplitude \(a\). The weights are randomly sampled: biases from \(U[-1,1]\), off-diagonal kernel weights from \(U[-5, 5]\), and diagonal kernel weights from \(U[-5, -0.1]\).
Methods compared. Since our objective is to highlight the performance gap between our path-wise gradient estimator (\(\partial\)SNN) and the score function gradient estimator (SNN), we use the same hyperparameters and initialization for both of them as much as possible. We initialize their parameters randomly using the same random seed so that both of them have random but the same initial parameters. We also set their hyperparameters as Table 1. The temperature is the only hyperparameter that impacts the performance gap. In preliminary experiments, we observe no significant impact for \(\tau\in[0.1, 0.5]\), and we only report the result at \(\tau=0.3\).
For the ground-truth network topology, letting the first two dimensions be observable and the last two dimensions be hidden, we consider a fully-connected architecture with the following weight and bias parameters: \[\begin{align} \bar{\mathbf{u}} &= \begin{bmatrix} 0 & 0 & 0 & 0 \end{bmatrix}^\top,\\ W_1 &= \begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 \end{bmatrix},\\ W_2 &= \begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 \end{bmatrix}, \end{align}\] where the \((d^\prime, d)\)-th element of \(W_l\) corresponds to \(w_{d^\prime, d, l}\) for \(d, d^\prime\in[D]\) and \(l=1, 2\).
For SNN to be trained, we initialize the bias term uniform randomly from \([0, 1]\) and the diagonal component of \(W_l\) from \([-1, -0.1]\). We set the upperbound of the diagonal components to be \(-0.1\) and at each parameter update we clip them so that they satisfy the constraint.
5.1 Variance of the Gradient Estimators↩︎
First, let us study the variance of the gradient estimators.
Protocol. We generate a single random parameter setting and use it to generate a synthetic data set consisting of \(10\) examples of length \(50\). Then, we compute the gradient estimators using the whole data set \(1000\) times, which yields \(1000\) gradient estimates for each method. Finally, we compute the standard deviations of each element of the gradients, and report the mean of the standard deviations.
Result. The mean standard deviation of \(\partial\)SNN was \(66.3\), whereas that of SNN was \(2.49\times 10^3\). This clearly demonstrates that the variance of our estimator tends to be lower than that of the existing estimator.
5.2 Predictive Performance↩︎
The second experiment studies the predictive performance of the models learned by each of the methods compared.
Protocol. We generate \(24\) random parameter settings, and consistently use them in this experiment. We aim to evaluate the performance gap between SNN and \(\partial\)SNN in different sizes of training sets. To this end, we execute the following, varying the size as \(N_{\mathrm{train}}=10\), \(20\), \(30\), \(40\), \(50\), \(75\) \(100\), \(200\), and for each parameter setting.
We generate training/test sets consisting of \(N_{\mathrm{train}}\)/\(100\) examples of length \(50\) respectively. SNN and \(\partial\)SNN are trained on the training set using AdaGrad [16] with initial learning rate \(0.05\) for \(10\) epochs. We evaluate the predictive performance by computing ELBO (Equation 9 ) on the test set. For fair comparison, we evaluate the performance of \(\partial\)SNN by transferring its parameters to the vanilla SNN5. By repeating this over \(24\) parameter settings, we obtain \(24\) ELBO scores. We report their mean as the performance of each method for each \(N_{\mathrm{train}}\).
Result. Figure 6 summarizes the experimental results. It clearly shows that \(\partial\)SNN consistently outperforms SNN especially in the small-sample regime, which supports the benefit of our low-variance estimator.
5.3 Computational Overhead↩︎
The last experiment studies computation overhead of \(\partial\)SNN over SNN. The computation time depends on the number of spikes, and the number of (hidden) spikes is proportional to \(a\), the amplitude of the non-linearlity \(\sigma\) that maps the membrane potential into the conditional intensity function. In general, \(\partial\)SNN generates more hidden spikes than SNN because the thinning algorithm for the differentiable point process does not reject any of the candidate spikes. Therefore, we expect that \(\partial\)SNN requires more computation time than SNN. The purpose of this experiment is to measure the computational overhead of \(\partial\)SNN over SNN.
Protocol. We generate a single parameter setting, and generate a training set of \(10\) examples of length \(50\). We then set up both SNN and \(\partial\)SNN with amplitude \(a=1,2,\dots,20\), resuting in \(40\) models to be trained. For each model, we measure the computation time of running \(100\) epochs, and obtain per-epoch computation time by averaging them.
Result. Figure 7 summarizes the experimental results. As is expected, \(\partial\)SNN requires \(2.8\) times more computation time than SNN on average. This result can be used as a reference for users to decide which gradient estimator to be employed. If a user can afford this overhead, our path-wise gradient estimator is recommended; otherwise, please consider to use the score function gradient estimator.
Note that we can improve the computation time of our method by introducing an adaptive upperbound \(\bar{\lambda}\) in Algorithm 5, if it is a tigher upperbound than the fixed upperbound. We leave this improvement as future work.
6 Related Work↩︎
The present work is related to the communities of SNNs and point processes. Let us discuss our contributions to them.
6.1 Spiking Neural Networks↩︎
The most relevant work is the stochastic variational learning algorithm for SNNs [1]. As discussed in Section 4.2.2, the difference is the gradient estimator. The authors used the score function gradient estimator, because the path-wise gradient estimator (which became popular by VAE [14]) was not popular at that time and the reparameterization trick for point processes was not trivial. Our contribution is to develop a differentiable point process that enables us to derive the path-wise gradient estimator. Since it is widely acknowledged that the reparameterization trick has often empirically less variance than the score function estimator, we develop a reparameterizable point process and apply it to the variational learning algorithm for SNNs. In our empirical study, we empirically validate that our algorithm achieves better generalization ability than that reported in Reference [1].
Less relevant but still worth mentioning are the line of work in learning algorithms for deterministic SNNs, where a neuron fires when the membrane potential exceeds a threshold. Although our technique cannot directly contribute to them, we believe it is worthwhile to compare the pros and cons of these different approaches for further development. Of a number of approaches proposed so far [17], we introduce two inspiring studies.
SpikeProp [18] is one of the earliest attempts to develop a learning algorithm for deterministic SNNs. SpikeProp uses backpropagation to minimize the difference between the target firing times \(\{t_n^{\star}\}_{n=1}^N\) and the actual firing times \(\{t_n\}_{n=1}^N\) of the network, i.e., \(\sum_{n=1}^N |t_n^\star - t_n|^2\). The gradient is approximated by assuming a linear relationship between the firing time and the membrane potential, which is valid only for a small learning rate.
[4] ([4]) propose a differentiable alternative to the threshold-based spike generation, which facilitates gradient computation. They employ a soft-threshold mechanism, and therefore, is differentiable without approximation. Another important contribution is that their model can handle not only spike trains but also a real-valued time-series. They use a readout network that maps spike trains from/into a real-valued time-series. This end-to-end formulation is significant towards practical applications of SNNs, and probabilistic SNNs should be equipped with this feature.
One interesting feature of probabilistic SNNs including our method is that both inference and learning algorithms can be executed naturally in an event-based manner without any discretization of time axis. This is in contrast to deterministic SNNs, where many learning algorithms require us to discretize the continuous-time dynamics for simulation.
6.2 Differentiable Point Processes↩︎
Our differentiable point process is significant in the community of point processes in that it largely expands the applicability of the reparameterization trick for point processes. Let us review the approaches to differentiable point processes, and discuss their pros and cons.
There are mainly three approaches to sample from point processes, and each of them can be used as a basis of differentiable point processes. The first approach [19] is to model the inter-event time conditioned on the past history by a log-normal mixture model, instead of modeling the conditional intensity function. Since it is straightforward to develop a reparameterizable sampling algorithm for the mixture model, the resultant point process is also reparameterizable. The second one is the inverse method [8], which utilizes the fact that the inverse of the compensator \(\Lambda^{[0,t]}\) can convert a unit-rate Poisson process into the point process with the corresponding conditional intensity function. [20] ([20]) propose a reparameterization trick based on the inverse method. The third one is the thinning algorithm, as we presented.
Of these three approaches, it is interesting to compare the second and the third approaches. When applying the inverse method [20] to computing ELBO, it is reported that the objective function contains discontinuous points, making optimization difficult. The discontinuity arises because time stamps of a realization are parameterized, and the algorithm involves a discrete decision whether a time stamp is less than \(T\) or not for termination. In contrast, Our differentiable point process does not suffer from it because not time stamps but marks are parameterized. In this sense, these two approaches are complementary.
When developing a path-wise gradient estimator for SNNs, only the third approach is feasible. The first approach is difficult to be applied because SNNs are modeled via the conditional intensity function, and the inter-event time distribution is not available in a closed form. The second approach is also difficult due to the lack of a closed-form expression of the inverse of the compensator. Our approach only assumes the existence of an upperbound of the conditional intensity function, and therefore, can be applied to SNNs. The assumption on the existence of a constant upperbound can be relaxed in the same way as Ogata’s method [10], which determines \(\bar{\lambda}\) adaptively.
7 Conclusion and Future Work↩︎
We develop a path-wise gradient estimator for SNNs based on a differentiable point process. Given the experimental results in Section 5, we conclude that our estimator has lower variance than the existing one, which contributes to improve the learning capability.
Throughout this paper, we only focus on the dependency of the gradient estimator on learning capability, and we have not discussed about its practical applications. In the community of SNNs, however, an increasing number of studies have started to apply SNNs to real-world tasks [21], [22]. One of the major concerns towards applying our method to real-world tasks is a method to convert real-valued data into/from spike trains. While there are a number of information encoding methods for spike trains, it is still an open problem which encoding is preferred. One interesting direction is to empirically and theoretically investigate the performance of different encoding methods and to understand their pros and cons.
Another limitation is the computational overhead as discussed in Section 5.3. While the probabilistic formulation can be simulated by an event-based manner, the gradient computation involves backpropagation through time (BPTT), whose complexity increases proportionally to the number of spikes. In addition to relying on the adaptive upperbound \(\bar{\lambda}\), applying online BPTT calculation and its approximation techniques [23] to SNNs may be an interesting research direction.
8 Relation to the Standard Notation↩︎
Our notation is different from the standard one employed by many others including a textbook [7]. A major difference is the conditional intensity function as illustrated in Figure 8. Our conditional intensity function \(\lambda(t\mid \mathcal{T}^{\leq t_n})\) is defined for \(t>t_n\) and is consistently conditioned on the history of events up to \(t_n\) (the top panel of Figure 8). On the other hand, the standard conditional intensity function \(\lambda^\ast(t)\) is defined for all \(t\), and the history that conditions it is ambiguous and depends on the context; it is sometimes conditioned on the history of events that occured before (and not including) \(t\), which is represented as \(\mathcal{T}^{\leq t_{n(t)}}\) (the bottom panel of Figure 8), and it is sometimes equivalent to our definition of the conditional intensity function.
While such an ambiguity helps to simplify equations (e.g., the compensator \(\Lambda^{[0,T]}=\int_0^T \lambda^\ast(t)\dl{t}\) can be represented by a single integral), it is sometimes very confusing especially for those who are not familiar with temporal point processes. Therefore, we employ a less ambiguous definition. In order to represent the standard conditional intensity function by our conditional intensity function, we introduced the left-continuous counting process \(n(t)\) as illustrated in Figure 1; with this counting process, we obtain the relationship between the standard and our conditional intensity functions, \(\lambda^\ast(t)=\lambda(t\mid\mathcal{T}^{\leq t_{n(t)}})\).
9 Conditions for Conditional Intensity Function↩︎
This section provides a proof of Proposition 2, which states several conditions under which the conditional intensity function can specify a marked point process uniquely.
Proof of Proposition 2. Since Equation 3 reads as, \[\begin{align} \int_X \dl{\mathbf{p}} \lambda(t,\mathbf{p}\mid \mathcal{T}_X^{\leq t_n}) &= \dfrac{f\left(t \mid \mathcal{T}_{X}^{\leq t_{n}}\right)}{1 - F\left(t \mid \mathcal{T}_{X}^{\leq t_{n}}\right)}\\ &= - \diff*{\log(1 - F(t \mid \mathcal{T}_{X}^{\leq t_n}))}{t}, \end{align}\] we can represent the cumulative distribution function by the conditional intensity function as, \[\begin{align} \label{eq:22}F(t \mid \mathcal{T}_X^{\leq t_n}) = 1 - \exp\left({-\int_{t_n}^t \dl{s} \int_{X}\dl{\mathbf{p}} \lambda(s,\mathbf{p}\mid \mathcal{T}_{X}^{\leq t_n})}\right). \end{align}\tag{14}\] In order for the conditional intensity function to define a proper cumulative distribution function, the function \(F(t \mid \mathcal{T}_{X}^{\leq t_n})\) defined as Equation 14 must satisfy the following four conditions:
\(\displaystyle\lim_{t\rightarrow\infty} F(t \mid \mathcal{T}_{X}^{\leq t_n}) = 1\),
\(\displaystyle\lim_{t\rightarrow -\infty} F(t \mid \mathcal{T}_{X}^{\leq t_n}) = 0\),
\(F(t \mid \mathcal{T}_{X}^{\leq t_n})\) is non-decreasing in \(t\), and
\(F(t \mid \mathcal{T}_{X}^{\leq t_n})\) is right-continuous.
Condition 1 is satisfied by Assumption 2. Condition 2 is satisfied because \(F(t \mid \mathcal{T}_X^{\leq t_n})=0\) holds for any \(t\leq t_n\). Condition 3 holds because \(\lambda(t, \mathbf{p})\geq 0\) implies that the integral in the exponential function in Equation 14 is non-decreasing. Condition 4 holds because Assumption 3 ensures that the exponent in Equation 14 is right-continuous, which indicates that Equation 14 itself is right-continuous. ◻
10 Properties of Differentiable Point Processes↩︎
This section provides several properties of differentiable point processes and their proofs. First, let us prove Proposition 4, which clarifies the conditional intensity function of a differentiable point process.
Proof of Proposition 4. The following calculus clarifies the relationship: \[\begin{align} \nonumber & \Pr\left[t_{n+1}\in[t, t+\dl{t}], \mathbf{p}_{n+1}=\mathbf{p}\mid \mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{<t}\right]\\ \nonumber = & \Pr\left[t_{n+1}\in[t, t+\dl{t}]\mid \mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{<t}\right] \\ \nonumber &\cdot p\left(\mathbf{p}_{n+1}=\mathbf{p}\mid t_{n+1}\in[t, t+\dl{t}], \mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{< t}\right)\\ = & \bar{\lambda} \dl{t} \cdot g_\tau\left(\begin{bmatrix} \mathbf{p}\\ 1- \|\mathbf{p}\|_1 \end{bmatrix}; \boldsymbol{\pi}_{\bar{\lambda}}\circ\boldsymbol{\lambda}\left(t \mid \mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{\leq t_n}\right)\right), \end{align}\] where let \(\mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{<t}\) denote the event \(t_{n+1} \notin (t_n, t)\) and \(\mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{\leq t_n}\). This suggests that in a differentiable point process, time stamps are distributed according to the homogeneous Poisson process with intensity \(\bar{\lambda}\), and each mark is distributed according to the concrete distribution. ◻
We then investigate two properties of the differentiable point process. Proposition 8 states that a realization of the differentiable point process is differentiable with respect to model parameters under mild conditions. Proposition 9 states that the differentiable point process becomes equivalent to the original point process as temperature goes zero.
Proposition 8. Let \(\lambda_\theta(t, \mathbf{p}\mid \mathcal{T}_{\mathbb{1}^D}^{\leq t_n})\) be a conditional intensity function of a multivariate point process parameterized by \(\theta\). Assume that the conditional intensity function can be calculated with \(\mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{\leq t_n}\) and is differentiable with respect to any mark in \(\mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{\leq t_n}\) and \(\theta\). Then, the marks of a realization of the corresponding differentiable point process is differentiable with respect to \(\theta\).
Proof. We prove Proposition 8 by induction. The first mark is distributed according to \(\boldsymbol{\pi}_{\bar{\lambda}}\circ\boldsymbol{\lambda}_\theta\left(t \mid \emptyset\right)\), which is differentiable with respect to \(\theta\). Assume that marks observed up to (but not including) time \(t\) are differentiable with respect to \(\theta\). A mark \(\mathbf{p}\) at time \(t\) is a realization of the concrete distribution with parameter \(\boldsymbol{\pi}:=\boldsymbol{\pi}_{\bar{\lambda}}\circ\boldsymbol{\lambda}_\theta\left(t \mid \mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{\leq t_n}\right)\), and thus, is differentiable with respect to \(\boldsymbol{\pi}\). \(\boldsymbol{\pi}\) is differentiable with respect to \(\boldsymbol{\lambda}_\theta\), and \(\boldsymbol{\lambda}_\theta\) is assumed to be differentiable with respect to \(\theta\) and the past marks, which are differentiable by assumption, and this completes the proof. ◻
Proposition 9. Assume that \(\lambda\left(t, \mathbf{p}\mid \mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{\leq t_n}\right)=\lambda\left(t, \mathbf{p}\mid \mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{\leq t_n}\backslash\{(t_k, \mathbf{p}_k)\}\right)\) holds for any \(\mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{\leq t_n}\) and any \(k\in[n]\) such that \(\mathbf{p}_k=\mathbf{0}\). Then, in the limit of \(\tau\rightarrow +0\), the output of Algorithm 4 is distributed according to \(\mathcal{MPP}(\lambda)\) if we discard the event with mark \(\mathbf{p}=\mathbf{0}\).
Proof. As proven by [11] [11], the random variable following the concrete distribution converges to the one-hot representation of the categorical variate in the limit of \(\tau\rightarrow +0\). The random variable in line 6 of Algorithm 4 satisfies the following, \[\begin{align} \Pr\left[\lim_{\tau\rightarrow +0} \begin{bmatrix}\mathbf{p}\\ r\end{bmatrix}=\mathbf{1}_d\mid s, \mathcal{T}\right] = \dfrac{\lambda(s, \mathbf{1}_d \mid \mathcal{T})}{\bar{\lambda}}, \end{align}\] for \(d\in[D+1]\), where let \(\lambda(s, \mathbf{1}_{D+1} \mid \mathcal{T}) \equiv \bar{\lambda} - \sum_{d=1}^{D}\lambda(s,\mathbf{1}_d\mid \mathcal{T})\). The above expression states that if the value of the conditional intensity function is the same, the random variable in line 6 of Algorithm 4 is equivalent to that in line 6 of Algorithm 3. The only difference between these algorithms is whether the event with zero mark \(\mathbf{p}=\mathbf{0}\) is discarded (Algorithm 3) or not (Algorithm 4). If the assumption made in Proposition 9 is satisfied, zero marks do not affect the conditional intensity function, and thus, the output of Algorithm 4 is distributed according to \(\mathcal{MPP}(\lambda)\) if we discard the events with zero marks. ◻
11 Properties of \(\partial\)SNNs↩︎
This section describes the properties of \(\partial\)SNNs. First, Proposition 5 states that the conditional intensity function (Equation 13 ) satisfies all of the conditions listed in Proposition 2, and thus, it defines an \(\bar{\mathcal{N}}\)-marked point process uniquely.
Proposition 10. Assume that the filter functions \(\{f_{d^\prime,d}(s)\}_{d,d^\prime\in[D]}\) are continuous with respect to \(s\). Then, the conditional intensity function (Equation 13 ) uniquely defines an \(\bar{\mathcal{N}}\)-marked point process.
Proof. We will confirm the assumptions of Proposition 2. Observing that, \[\begin{align} &\int_{t_n}^{t}\dl{s} \int_{\bar{\mathcal{N}}} \dl{\mathbf{p}} \lambda\left(s, \mathbf{p}\mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_n}\right) \\ = & \int_{t_n}^{t}\dl{s} \left[\int_{\mathrm{conv}_{0}(\mathcal{H})} \dl{\mathbf{p}} \lambda\left(s, \mathbf{p}\mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_n}\right) \right.\\ & \left.+ \sum_{\mathbf{p}\in\mathcal{O}} \lambda\left(s, \mathbf{p}\mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_n}\right) \right]\\ = & \int_{t_n}^{t}\dl{s} \left[\bar{\lambda} + \sum_{\mathbf{p}\in\mathcal{O}} \lambda\left(s, \mathbf{p}\mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_n}\right) \right]\\ = & \bar{\lambda}(t-t_n) + \int_{t_n}^{t}\dl{s} \sum_{\mathbf{p}\in\mathcal{O}} \lambda\left(s, \mathbf{p}\mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_n}\right), \end{align}\] the first condition is satisfied. By taking \(t\rightarrow\infty\) in the above expression, we can confirm that the second condition is satisfied (the first term \(\bar{\lambda}(t-t_n)\) goes to infinity and the second term is guaranteed to be non-negative). The third condition is satisfied because \(\lambda(s,\mathbf{p}\mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_n})\) is a continuous function with respect to \(s>t_n\), which can be guaranteed by the continuity of the filter functions. ◻
Then, let us discuss two properties of \(\partial\)SNN. Proposition 6 states that the objective function is differentiable with respect to \(\phi\). Proposition 7 states that \(\partial\)SNN becomes equivalent to the vanilla SNN in the limit of \(\tau\rightarrow +0\).
Proposition 11. Assume that the filter functions \(\{f_{d^\prime,d}(s)\}_{d,d^\prime\in[D]}\) are differentiable with respect to their parameters \(\phi\). The Monte-Carlo approximation of ELBO (Equation 9 ) is differentiable with respect to \(\phi\) if we employ \(\mathcal{\partial PP}(\lambda_q(t, \mathbf{p}\mid \mathcal{T}_{\bar{\mathcal{N}}};\phi); \bar{\lambda}, \tau)\) as the variational distribution.
Proof. We first show that marks of a realization of the variational distribution \(\mathcal{T}_{\mathcal{H}}(\phi)\) are differentiable with respect to \(\phi\). We then show that both \(\log p(\mathcal{T}_{\mathcal{O}}, \mathcal{T}_{\mathrm{conv}_{0}(\mathcal{H})})\) and \(\log q(\mathcal{T}_{\mathrm{conv}_{0}(\mathcal{H})})\) are differentiable with respect to the marks of \(\mathcal{T}_{\mathcal{H}}(\phi)\). We finally show that \(\log q(\mathcal{T}_{\mathrm{conv}_{0}(\mathcal{H})};\phi)\) is differentiable with respect to \(\phi\). By the fact that the composition of differentiable functions is differentiable, the proposition is implied by these three statements.
\(\mathcal{T}_{\mathcal{H}}(\phi)\) is sampled by Algorithm 4 using the conditional intensity function \(\lambda_q(t,\mathbf{p}\mid \mathcal{T}_{\bar{\mathcal{N}}};\phi)\), which is differentiable with respect to any mark in \(\mathcal{T}_{\bar{\mathcal{N}}}\) and \(\phi\) (by the assumption). Therefore, by Proposition 8, the marks of \(\mathcal{T}_{\mathcal{H}}(\phi)\) are differentiable with respect to \(\phi\).
The logarithm of the joint distribution \(\log p(\mathcal{T}_{\mathcal{O}},\mathcal{T}_{\mathrm{conv}_{0}(\mathcal{H})})\) can be written as, \[\begin{align} &\log p(\mathcal{T}_{\mathcal{O}},\mathcal{T}_{\mathrm{conv}_{0}(\mathcal{H})}; \theta) \\ = &\sum_{(t,\mathbf{p})\in\mathcal{T}_{\mathcal{O}}} \log \lambda^{\mathrm{SNN}}\left(t, \mathbf{p}\mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_{n(t)}}\right) \\ & + \sum_{(t,\mathbf{p})\in\mathcal{T}_{\mathrm{conv}_{0}(\mathcal{H})}} \log \lambda_{\partial}\left(t, \mathbf{p}_{\mathcal{H}} \mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_{n(t)}};\boldsymbol{\lambda}_\mathcal{H}, \bar{\lambda}, \tau\right)\\ & - \bar{\lambda} T - \int_{0}^T \dl{s}\sum_{\mathbf{p}\in\mathcal{O}}\lambda^{\mathrm{SNN}}\left(s,\mathbf{p}\mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_{n(s)}}\right). \end{align}\] The first and the last terms are differentiable with respect to marks of \(\mathcal{T}_{\mathrm{conv}_{0}(\mathcal{H})}\), and the third term does not depend on \(\mathcal{T}_{\mathrm{conv}_{0}(\mathcal{H})}\). For any \((t,\mathbf{p})\in\mathcal{T}_{\mathrm{conv}_{0}(\mathcal{H})}\), the summand of the second term, \[\begin{align} & \log \lambda_{\partial}\left(t, \mathbf{p}_{\mathcal{H}} \mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_{n(t)}};\boldsymbol{\lambda}_\mathcal{H}, \bar{\lambda}, \tau\right)\\ = & \log\bar{\lambda} + \log g_\tau\left( \begin{bmatrix} \mathbf{p}_{\mathcal{H}}\\ 1- \|\mathbf{p}_{\mathcal{H}}\|_1 \end{bmatrix}; \boldsymbol{\pi}_{\bar{\lambda}}\circ\boldsymbol{\lambda}_{\mathcal{H}}\left(t \mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_n(t)}\right)\right), \end{align}\] is differentiable with respect to \(\mathbf{p}_\mathcal{H}\) because the probability density function of the concrete distribution is differentiable with respect to \(\mathbf{p}_\mathcal{H}\). It is also differentiable with respect to the past marks in \(\mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_{n(t)}}\) because the probability density function \(g_\tau\) is differentiable with respect to its parameter \(\boldsymbol{\pi}:=\boldsymbol{\pi}_{\bar{\lambda}}\circ\boldsymbol{\lambda}_\mathcal{H}\left(t \mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_n(t)}\right)\), which is differentiable with respect to the marks in \(\mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_n(t)}\). Therefore, \(\log p(\mathcal{T}_\mathcal{O}, \mathcal{T}_{\mathrm{conv}_{0}(\mathcal{H})};\theta)\), is differentiable with respect to the marks of \(\mathcal{T}_{\mathrm{conv}_{0}(\mathcal{H})}\).
The logarithm of the variational distribution \(\log q(\mathcal{T}_{\mathrm{conv}_{0}(\mathcal{H})})\) can be written as, \[\begin{align} \nonumber &\log q(\mathcal{T}_{\mathrm{conv}_{0}(\mathcal{H})}) \\ \label{eq:2} =& \sum_{(t,\mathbf{p})\in\mathcal{T}_{\mathrm{conv}_{0}(\mathcal{H})}} \log \lambda_{\partial}\left(t, \mathbf{p}_{\mathcal{H}} \mid \mathcal{T}_{\bar{\mathcal{N}}}^{\leq t_{n(t)}};\boldsymbol{\lambda}_\mathcal{H}, \bar{\lambda}, \tau\right) - \bar{\lambda}T, \end{align}\tag{15}\] which is also differentiable with respect to the marks of \(\mathcal{T}_{\mathrm{conv}_{0}(\mathcal{H})}\) in the same way as the above discussion.
Finally, \(\log q(\mathcal{T}_{\mathrm{conv}_{0}(\mathcal{H})};\phi)\) is differentiable with respect to \(\phi\), because the first term of Equation 15 is differentiable with respect to \(\boldsymbol{\lambda}_{\mathcal{H}}\), which is differentiable with respect to \(\phi\) (partially by the assumption that the filter functions are differentiable with respect to \(\phi\)).
The proposition follows from the three statements above. ◻
Proposition 12. In the limit of \(\tau\rightarrow +0\), a realization of \(\partial\)SNN (Equation 13 ) is distributed according to the vanilla SNN (Equation 7 ) if we discard events with mark \(\mathbf{p}=\mathbf{0}\).
Proof. Since for any event \((t_k,\mathbf{p}_k)\) (\(k\in[n]\)) such that \(\mathbf{p}_k=\mathbf{0}\), \[\begin{align} \lambda^{\mathrm{SNN}}\left(t,\mathbf{p}\mid \mathcal{T}_{\bar\mathcal{N}}^{\leq t_n}\right) = \lambda^{\mathrm{SNN}}\left(t,\mathbf{p}\mid \mathcal{T}_{\bar\mathcal{N}}^{\leq t_n}\backslash\{(t_k,\mathbf{p}_k)\}\right), \end{align}\] holds, the event with mark \(\mathbf{p}=\mathbf{0}\) has no influence on the conditional intensity function. In the limit of \(\tau\rightarrow +0\), \[\begin{align} & {\lambda^{\partial\mathrm{SNN}}\left(t,\mathbf{p}\mid \mathcal{T}_{\bar\mathcal{N}}^{\leq t_n}\right)}\\ = & \begin{cases} \lambda^{\mathrm{SNN}}\left(t,\mathbf{p}\mid \mathcal{T}_{\bar\mathcal{N}}^{\leq t_n}\right) & (\mathbf{p}\in\mathcal{O}\cup\mathcal{H}),\\ \displaystyle\left(\bar{\lambda} - \sum_{\mathbf{p}\in\mathcal{H}} \lambda^{\mathrm{SNN}}\left(t,\mathbf{p}\mid \mathcal{T}_{\bar\mathcal{N}}^{\leq t_n}\right)\right) & (\mathbf{p}=\mathbf{0}), \end{cases} \end{align}\] holds. As discussed above, the event with mark \(\mathbf{p}=\mathbf{0}\) has no influence on computing the conditional intensity function, and can be removed without changing the conditional intensity function. Therefore, a realization of \(\mathcal{MPP}(\lambda^{\partial\mathrm{SNN}})\) is equivalent to that of \(\mathcal{MPP}(\lambda^{\mathrm{SNN}})\) if we discard all the events with mark \(\mathbf{p}=\mathbf{0}\). ◻
12 Numerically Stable Implementation↩︎
For numerical stability, we recommend to represent all probability and subprobability vectors and conditional intensity functions in the logarithmic scale. In this section, we describe an accurate computation of the parameter of the concrete distribution.
When computing the parameter of the concrete distribution shown below in the logarithmic scale, \[\begin{align} \begin{aligned} &\boldsymbol{\pi}_{\bar{\lambda}}\circ\boldsymbol{\lambda}\left(t \mid \mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{\leq t_n}\right)\\ = &\frac{1}{\bar{\lambda}} \begin{bmatrix} \boldsymbol{\lambda}\left(t \mid \mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{\leq t_n}\right) & \bar{\lambda} - \left\|\boldsymbol{\lambda}\left(t \mid \mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{\leq t_n}\right)\right\|_1 \end{bmatrix}, \end{aligned} \end{align}\] it is straightforward to compute the first \(D\) elements with the logarithm of the conditional intensity function, \(\log\lambda\left(t, \mathbf{p}\mid \mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{<t}\right)\). However, the last element, \[\begin{align} \label{eq:3} 1-\frac{\sum_{\mathbf{p}^\prime\in \mathbb{1}^D}\lambda\left(t, \mathbf{p}^\prime \mid \mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{<t}\right)}{\bar{\lambda}}, \end{align}\tag{16}\] is not trivial to compute accurately in the logarithmic scale.
We resort to \(\mathtt{log1mexp}\) [24] to compute it, which allows us to compute \(\log(1 - \exp(-|a|))\) accurately for \(a\neq 0\) as follows: \[\begin{align} \mathtt{log1mexp}(a) = \begin{cases} \log(-\mathtt{expm1}(-a)) & (0<a\leq \log 2),\\ \mathtt{log1p}(-\exp(-a)) & (a> \log 2), \end{cases} \end{align}\] where \(\mathtt{expm1(x)}\) and \(\mathtt{log1p}(x)\) approximately compute \(\exp(x)-1\) and \(\log(1+x)\) respectively by using a few terms of their Taylor series. Since we can compute \(\log p := \log\left[{\sum_{\mathbf{p}^\prime\in \mathbb{1}^D}\lambda\left(t, \mathbf{p}^\prime \mid \mathcal{T}_{\mathrm{conv}_{0}(\mathbb{1}^D)}^{<t}\right)}\right] - \log{\bar{\lambda}}\) by \(\mathtt{logsumexp}\), the computation of Equation 16 boils down to the computation of \(\log (1-p)\) given \(\log p\) (\(p\in[0,1)\)). Since \(\log(1-p) = \log(1-\exp(-|\log p|))\) holds for \(p\in[0,1)\), we can utilize \(\mathtt{log1mexp}\) to compute it.
References↩︎
Biological plausibility is of another great interest, but is not mentioned because it is not mandatory for engineering purposes.↩︎
A counting process is usually defined to be right-continuous. We introduce the left-continuous one so as to represent the integrand of the compensator concisely. See Appendix 8 for details.↩︎
We use the sigmoid function multiplied by amplitude \(a>0\) element-wisely as \(\sigma\), for which \(\bar{\lambda}\) is easy to derive.↩︎
Intel Xeon Gold 6248 2.50GHz 48 cores and 192GB memory.↩︎
For better transfer, we decrease \(\tau\) geometrically by ratio \(0.95\) at every epoch, which slightly improves the performance.↩︎