SNN系列｜学习算法篇(1)Tempotron

Tempotron (论文传送门)

Tempotron是一个二层网络学习算法，输入脉冲序列，输出脉冲响应。对于二分类，最重要的是保证以下关系的存在：即应该发放脉冲的神经元的最大膜电势应超过阈值，否则就增加突触传递效率；反之，不该发放脉冲的其最大膜电势就不该超过阈值，否则就减小突触传递效率。
$$V\left(t_{\max }^{\oplus }\right)>V_{\mathrm{t}\mathrm{h}\mathrm{r}}>V\left(t_{\max }^{\ominus }\right)$$

• 神经元模型LIF

这里的LIF神经元模型使用的是连续形式，因为Tempotron是一个基于梯度下降的SNN优化算法，后续可以很方便地直接求导。
$$V(t)=\sum _i{\omega }_i\sum _{t_i}K\left(t-t_i\right)+V_{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}}$$

$$K(t-t_i)=V_0\left(\exp \left[-\left(t-t_i\right)/\tau \right]-\exp \left[-\left(t-t_i\right)/{\tau }_{\mathrm{s}}\right]\right)$$

$t_i$表示脉冲发放时刻，$\tau$表示膜电势整合的延迟参数，${\tau }_s$表示突触电流的延迟参数，$V_0$表示归一化因子，它使得PSP取决于突触效率$\omega$，$V_0$的计算步骤为：设$V_0$为0，求最大$K$所在时刻，然后代入$t_{max}$求方程的倒数即为$V_0$。
$$V_0=\frac{\tau {\tau }_s}{\tau -{\tau }_s}\ln \frac{\tau }{{\tau }_s}$$
结合$y=e^{-x}$函数特性（$x$越大衰减越快）以及脉冲的基本形式（钟形），$\tau$势必是大于${\tau }_s$的，原文取$\tau /{\tau }_s=4$，且比例越大，膜电势上升和下降所需时间比越小（尾巴越长）。

LIF神经元代码如下：

def compute_norm_factor(self, tau, tau_s):
    tmax = (tau * tau_s * np.log(tau/tau_s)) / (tau - tau_s)
    v_max = self.K(1, tmax, 0)
    V_0 = 1/v_max
    return V_0

def K(self, V_0, t, t_i):
		if t < t_i:
    		value = 0
    else:
        value = V_0 * (np.exp(-(t-t_i)/self.tau) - np.exp(-(t-t_i)/self.tau_s))
    return value
  
def compute_spike_contributions(self, t, spike_times):
    N_synapse = len(spike_times)
    spike_contribs = np.zeros(N_synapse)
    for neuron_pos in xrange(N_synapse):
        for spike_time in spike_times[neuron_pos]:
            spike_contribs[neuron_pos] += self.K(self.V_norm, t, spike_time)
    return spike_contribs
  
def compute_membrane_potential(self, t, spike_times):
    spike_contribs = self.compute_spike_contributions(t, spike_times)
    total_incoming = spike_contribs * self.efficacies
		# add sum and add V_rest to get membrane potential
    V = total_incoming.sum() + self.V_rest
		return V

• 网络训练

Tempotron只在错误发生时更新突触传递效率，其更新规则是梯度下降，损失函数如下：
$$\begin{gathered} E_{\pm }=\pm \left(V_{\text{thr }}-V\left(t_{\text{max }}\right)\right)\Theta \left(\pm \left(V_{\text{thr }}-V\left(t_{\text{max }}\right)\right)\right) \\ \Theta (x)=\{\begin{array}{l}1,x\ge 0\\ 0,x<0\end{array} \end{gathered}$$
损失函数对传递效率（权重）求偏导，
$$-\frac{\mathrm{d}{E}_{\pm }}{\mathrm{d}{\omega }_i}=\pm \sum _{t_i<t_{\max }}K\left(\Delta t_i\right)\pm \frac{\partial V\left(t_{\max }\right)}{\partial t_{\max }}\frac{d{t}_{\max }}{d{\omega }_i}$$
因为$t_{max}$的定义是膜电势最大的时刻，因此突触传递效率的更新规则如下：
$$\Delta {\omega }_i=\lambda \sum _{t_i<t_{\text{max }}}K\left(t_{\max }-t_i\right)$$鉴于Tempotron应用于二层网络，其输出希望对$\oplus$模式发放脉冲，而对$\ominus$模式不发放脉冲，因此对应于最大膜电势要大于脉冲发放阈值，而对另一种模式小于脉冲发放阈值。注意：在T时间内，输出神经元只负责计算膜电势，只在最后判断响应对错时才将其与阈值进行比较。也就是说，Tempotron算法的输出层最多只发放一个脉冲（发不发无所谓，我们关心的是最大膜电势），因此它也无法扩展到多层神经网络中使用。

问题就到了如何判断最大膜电势上来，可以将LIF神经元膜电势对时间进行求导，一定在导数为0的那一些时刻序列里，这也就是为什么LIF神经元模型使用连续绝对形式，而不是离散增量形式的原因，求导可得最大膜电势时刻为
$$t_{max}=\frac{\tau {\tau }_s}{\tau -{\tau }_s}\left(\ln \frac{\tau }{{\tau }_s}+\ln \frac{\sum {\omega }_i\exp (\frac{t_i}{\tau })}{\sum {\omega }_i\exp (\frac{t_i}{{\tau }_s})}\right)$$
使用连续LIF形式的另一个好处就是不用在每个仿真时刻都计算膜电势，而只需要计算这些疑似最大膜电势时刻的膜电势即可，能有效加快运行。

参数更新如下：

def adapt_weights(self, spike_times, target, learning_rate):
		tmax = self.compute_tmax(spike_times)
    vmax = self.compute_membrane_potential(tmax, spike_times)
    
    if (vmax >= self.threshold) == target:
        return
      
		dw = self.dw(learning_rate, tmax, spike_times)
    if target is True:
        self.efficacies += dw
    else:
        self.efficacies -= dw

def dw(self, learning_rate, tmax, spike_times):
    spike_contribs = self.compute_spike_contributions(tmax, spike_times)
		update = learning_rate * spike_contribs
		return update

def compute_tmax(self, spike_times):
		spikes_chron = [(time, synapse) for synapse in xrange(len(spike_times)) for time in spike_times[synapse]]
    spikes_chron.sort()
      
    spikes = [(s[0], self.efficacies[s[1]]) for s in spikes_chron]
    times = np.array([spike[0] for spike in spikes])
    weights = np.array([spike[1] for spike in spikes])

    sum_tau = (weights*np.exp(times/self.tau)).cumsum()
    sum_tau_s = (weights*np.exp(times/self.tau_s)).cumsum()

    div = sum_tau_s/sum_tau
    boundary_cases = div < 0
    div[boundary_cases] = 10

    tmax_list = self.tau*self.tau_s*(self.log_tts + np.log(div))/(self.tau - self.tau_s)
    tmax_list[boundary_cases] = times[boundary_cases]

    vmax_list = np.array([self.compute_membrane_potential(t, spike_times) for t in tmax_list])

    tmax = tmax_list[vmax_list.argmax()]
    return tmax

如果觉得求导比较麻烦或者绕，可以使用离散LIF模型，计算仿真时刻内所有的膜电势，然后就可以很轻易地找到最大膜电势并使用梯度下降，但存储量稍微大一些。

受知乎Jay Wang的启发，想到了Tempotron的容量其实不大，不能使用大规模的训练集进行训练，除非有非常高效的编码方法，因为作者提到，Tempotron的类别信息不是嵌入在脉冲数和单个神经元的脉冲时序中，而是在同步性之中。当样本数量增多后，同类输入的产生的输出也会有很大差异。

以上求解过程均为解析过程，同样也可以使用迭代的思想求解，其他实验细节参考原论文或以后补充。

• 附：最大模电势时刻推导

• 这里将一个神经元的多个脉冲理解成多个相同权重的神经元分别发放一个脉冲，从而使得膜电势计算公式有如下变化
  $$\begin{gathered} V(t)=\sum {\omega }_i\sum _{t_i}{V}_0\left(\exp (-\frac{t-t_i}{{\tau }_s})-\exp (-\frac{t-t_i}{\tau })\right) \\ ⟹V(t)=V_0\sum {\omega }_i\left(\exp (-\frac{t-t_i}{{\tau }_s})-\exp (-\frac{t-t_i}{\tau })\right) \end{gathered}$$最大膜电势时刻其导数为0
  $$\begin{gathered} V^{\prime }(t)=V_0\sum {\omega }_i\left(\frac{1}{{\tau }_s}\exp (-\frac{t-t_i}{{\tau }_s})-\frac{1}{\tau }\exp (-\frac{t-t_i}{\tau })\right)=0 \\ ⟹\frac{1}{{\tau }_s}\sum {\omega }_i\exp (\frac{t_i}{{\tau }_s})\exp (-\frac{t_{max}}{{\tau }_s})=\frac{1}{\tau }\sum {\omega }_i\exp (\frac{t_i}{\tau })\exp (-\frac{t_{max}}{\tau }) \\ ⟹\exp (-\frac{t_{max}}{{\tau }_s})\frac{1}{{\tau }_s}\sum {\omega }_i\exp (\frac{t_i}{{\tau }_s})=\exp (-\frac{t_{max}}{\tau })\frac{1}{\tau }\sum {\omega }_i\exp (\frac{t_i}{\tau }) \\ ⟹-\frac{t_{max}}{{\tau }_s}+\ln \frac{1}{{\tau }_s}+\ln \sum {\omega }_i\exp (\frac{t_i}{{\tau }_s})=-\frac{t_{max}}{\tau }+\ln \frac{1}{\tau }+\ln \sum {\omega }_i\exp (\frac{t_i}{\tau }) \\ ⟹\frac{\tau -{\tau }_s}{\tau {\tau }_s}{t}_{max}=\ln \frac{\tau }{{\tau }_s}+\ln \frac{\sum {\omega }_i\exp (\frac{t_i}{{\tau }_s})}{\sum {\omega }_i\exp (\frac{t_i}{\tau })} \\ ⟹t_{max}=\frac{\tau {\tau }_s}{\tau -{\tau }_s}\left(\ln \frac{\tau }{{\tau }_s}+\ln \frac{\sum {\omega }_i\exp (\frac{t_i}{{\tau }_s})}{\sum {\omega }_i\exp (\frac{t_i}{\tau })}\right) \end{gathered}$$
  $\frac{\sum {\omega }_i\exp (\frac{t_i}{{\tau }_s})}{\sum {\omega }_i\exp (\frac{t_i}{\tau })}$是不能小于等于0的，不然$t_{max}$便不会存在，当出现小于0的情况时，手动设置其比值为10。

参考

[1] The tempotron: a neuron that learns spike timing–based decisions

[2] Python代码：dieuwkehupkes/Tempotron

[3] Matlab代码：laurence-lin/tempotron
[4] SNN-Tempotron

​