SpikingNN_on_trial

refs

1. Error-backpropagation in temporally encoded networks of spiking neurons, 10.1016/S0925-2312(01)00658-0
2. Spatio-Temporal Backpropagation for Training High-Performance Spiking Neural Networks, 10.3389/fnins.2018.00331

神经元的脉冲在信息编码中的特点

1. [29] W. Maass, Fast sigmoidal networks via spiking neurons, Neural Comput. 9 (2) (1997) 279–304.

   · 在上述文章中，已经给出证明，脉冲神经网络可以模拟任意形式的前馈sigmoid神经元组成的神经网络。 基于这个特性，SNN可以近似任意的连续函数。SNN中单个尖峰时间传递信息的神经元比具有sigmoid激活函数的神经元计算能力更强。

2. 脉冲，在数学上表征为一个“活人/游戏事件”，即“地点，时刻”的坐标。而且活跃（尖峰）神经元的数量通常比较稀少。基于这个特性，SNN可以允许神经网络大模型在VLSI上高效部署。

3. SNN会产生与Kohonen’s SOM类似的自组织，无监督聚类网

针对时空编码范式的有监督的单脉冲学习算法：基于误差-反向传播的有监督学习算法

概念：阈值函数？（threshold function）
速率编码网络？（rate-coded network）活/人游戏？
突触后电位（呃呃呃死去的回忆突然攻击我）

SNN的网络结构

脉冲响应模型：spike response model (SRM)，表征输入脉冲和神经元内部的状态的关系。

单个神经元(第j个)的 内在状态（inertial state） 的度量:

$$x_j(t)=\sum_{i \in \Gamma_j}w_{ij}\epsilon(t-t_i)$$

式中，$w_{ij}$代表了单个突触的效能（synaptic eJcacy (“weight”)）

这个式子也表征了一个神经元的未加权的post-synaptic potential（突触后电位）。

至此，我们有了一个基本的认知：神经元的突触前膜产生的若干个电位信号作用在神经元，其和作用达到一个阈值之后，就会产生动作电位，并向下传导。这和高中学过的神经生理学的知识很相似。

每个突触与前后两个神经元的连接建模如下：

$d_k =$
突触前神经元的激活时刻（动作电位刚刚发出的时刻） - 突触后电位刚开始上升的时刻

单个突触末端对神经元状态变量的非加权贡献值为

$$y_i^k(t) = \epsilon(t-t_i-d^k)$$

$\epsilon(t)$表示PSP脉冲响应函数，是一个因果信号。

$$\epsilon(t) = \frac{t}{\tau}e^{1-\frac{t}{\tau}}$$

式中$\tau$用来表征膜电位延迟常数。

至此，我们对单个突出的前后电位传递进行了建模。那么重写一下$x_j(t)$:

$$x_j(t) = \sum_{i \in \Gamma_j}\sum_{k=1}^{m} w_{ij}^k y_i^k(t)$$

$t_j$是神经元的激活时刻，为状态变量（膜电位）第一次超过阈值电位（$\vartheta$）的时间.阈值电位在整个神经网络里都是相等的。

SNN的误差反向传播

假设整个网络包括H（输入层）、I（隐藏层）、J（输出层）：

前馈算法的目的是学习一组目标激活时刻/放电时刻（firing times）：$t_j^d$

对于神经元$h \in H$, 给定一组输入范式$\{P[t_1, ..., t_h]\}$, 输出的每一个神经元$j \in J$都存在一组输出(firing time)。以最小均方误差为代价函数：(实际的放电时刻$t_j^a$和预估的/期待的放电时刻)

$$E = \frac{1}{2}\sum_{j \in J} (t_j^a-t_j^d)^2$$

对于误差反向传播，我们将每一个突触连接看作一个单独的作用单元k, 其权重是$w_{ij}^k$:

$$\Delta w_{ij}^k = -\eta \frac{\partial E}{\partial w_{ij}^k}$$

此步就是正常的反向传播 ?思路。$\eta$即学习率。

$$\frac{\partial E}{\partial w_{ij}^k}=\frac{\partial E}{\partial t_j}(t_j^\mathrm{a})\frac{\partial t_j}{\partial w_{ij}^k}(t_j^\mathrm{a})=\frac{\partial E}{\partial t_j}(t_j^\mathrm{a})\frac{\partial t_j}{\partial x_j(t)}(t_j^\mathrm{a})\frac{\partial x_j(t)}{\partial w_{ij}^k}(t_j^\mathrm{a})$$

1.

$t_j$作为一个突触后阈值在$t=t_j^a$附近的输入$x_j(t)$. 在$t=t_j^a$的邻域内，我们假设$x_j$近似关于t线性(很正常的假设)。

$$\left.\frac{\partial t_j}{\partial x_j(t)}(t_j^\mathrm{a})=\left.\frac{\partial t_j(x_j)}{\partial x_j(t)}\right|_{x_j=\vartheta}=\frac{-1}\alpha=\frac{-1}{\partial x_j(t)/\partial t(t_j^\mathrm{a})}=\frac{-1}{\sum_{i,l}w_{ij}^l(\partial y_i^l(t)/\partial t)(t_j^\mathrm{a})}.\right.$$

其中，$\alpha$为邻域内$x_j(t)$关于$t$的斜率。$x_j(t)$越大，到达电位阈值就会更容易，表征为$\delta t_j$越小。因此公式中出现负号。

2.

$$\frac{\partial E(t_j^\mathrm{a})}{\partial t_j^\mathrm{a}}=(t_j^\mathrm{a}-t_j^\mathrm{d}).$$

3.

$$\frac{\partial x_j(t_j^\mathrm{a})}{\partial w_{ij}^k}=\frac{\partial\{\sum_{n\in\Gamma_j}\sum_lw_{nj}^ly_n^l(t_j^\mathrm{a})\}}{\partial w_{ij}^k}=y_i^k(t_j^\mathrm{a}).$$

将上面三个子项代入原始的更新权重中：

$$\Delta w_{ij}^k(t_j^\mathrm{a})=-\eta\frac{y_i^k(t_j^\mathrm{a})(t_j^\mathrm{d}-t_j^\mathrm{a})}{\sum_{i\in\Gamma_j}\sum_lw_{ij}^l(\partial y_i^l(t_j^\mathrm{a})/\partial t_j^\mathrm{a})}.$$

为便于表示，定义超参数$\delta_j$

$$\delta_j\equiv\frac{\partial E}{\partial t_j^\mathrm{a}}\mathrm{~}\frac{\partial t_j^\mathrm{a}}{\partial x_j(t_j^\mathrm{a})}=\frac{(t_j^\mathrm{d}-t_j^\mathrm{a})}{\sum_{i\in\Gamma_j}\sum_lw_{ij}^l(\partial y_i^l(t_j^\mathrm{a})/\partial t_j^\mathrm{a})}$$

$$\Delta w_{ij}^k=-\eta y_i^k(t_j^\mathrm{a})\delta_j.$$

对于第J层之前的隐藏层，同理：对于第i个神经元（$i \in I$），其点火时刻为$t_i^a$,

$$\begin{aligned} \delta_{i}& \equiv\frac{\partial t_i^\mathbf{a}}{\partial x_i(t_i^\mathbf{a })}\frac{\partial E}{\partial t_i^\mathbf{a}} \\ &=\frac{\partial t_i^\mathrm{a}}{\partial x_i(t_i^\mathrm{a})}\sum_{j\in\Gamma^i}\frac{\partial E}{\partial t_j^\mathrm{a}}\frac{\partial t_j^\mathrm{a}}{\partial x_j(t_j^\mathrm{a})}\frac{\partial x_j(t_j^\mathrm{a})}{\partial t_i^\mathrm{a}} \\ &=\frac{\partial t_i^\mathrm{a}}{\partial x_i(t_i^\mathrm{a})}\sum_{j\in\Gamma^i}\delta_j\frac{\partial x_j(t_j^\mathrm{a})}{\partial t_i^\mathrm{a}}, \end{aligned}$$

$$\begin{gathered} \frac{\partial x_j(t_j^\mathrm{a})}{\partial t_i^\mathrm{a}} =\frac{\partial\sum_{l\in I}\sum_kw_{lj}^ky_l^k(t_j^\mathrm{a})}{\partial t_i^\mathrm{a}} \\ =\sum_kw_{ij}^k\frac{\partial y_i^k(t_j^\mathrm{a})}{\partial t_i^\mathrm{a}}. \end{gathered}$$

$$\delta_i=\frac{\sum_{j\in\Gamma^i}\delta_j\{\sum_kw_{ij}^k(\partial y_i^k(t_j^\mathrm{a})/\partial t_i^\mathrm{a})\}}{\sum h\in\Gamma_i\sum_lw_{hi}^l(\partial y_h^l(t_i^\mathrm{a})/\partial t_i^\mathrm{a})}.$$

$$\Delta w_{hi}^k=-\left.\eta y_h^k(t_i^\mathrm{a})\delta_i=-\eta\frac{y_h^k(t_i^\mathrm{a})\sum_j\{\delta_j\sum_kw_{ij}^k(\partial y_i^k(t_j^\mathrm{a})/\partial t_i^\mathrm{a})\}}{\sum_{n\in\Gamma_i}\sum_lw_{ni}^l(\partial y_n^l(t_i^\mathrm{a})/\partial t_i^\mathrm{a})}.\right.$$

（以上的所有推导全部基于2002年的误差-反向传播思路，但是太™麻烦，也很难理清思路，不急，有的是时间理清）

tutor：上面的模型有些许误导性，若干个突触连接前后两个神经元，看似是把网络维度扩大了，网络更复杂了。但是如今的SNN网络架构已经摒弃了这种想法。类似ANN，前后神经元仅通过一个权重参数进行前向传播。

我：huh？

error-bp SNN 为什么不好？

SNN发展的几个派系（ε=ε=ε=┏(゜ロ゜;)┛）

1. ANN-to-SNN：采用近似的思想，实际工程中会造成精度丢失。近似过后的SNN性能也不如ANN。

2. 膜电位驱动的学习算法(membrane potential-driven learning algorithms)：空间需求大。

3. 峰值驱动学习算法(spike-driven learning algorithms): 包括SpikeProp及其衍生算法。这类算法依赖神经元膜电位在点火时间邻域内增长近似线性，导数计算变得更容易。仍存在坏死的神经元和梯度爆炸问题（版本答案）

DeepSNN的建模过程以及其训练弊端

考虑一个全连接DeepSNN网络，为便于分析，假设所有神经元最多进行一次点火。$V_j^l$表示第$l$层的第$j$个神经元的膜电位（内在状态/inertial state）, 其突触后连接表示为：

$$V_j^l(t)=\sum_i^N\omega_{ij}^l\varepsilon\left(t-t_i^{l-1}\right)-\eta(t-t_j^l)$$

参数解释：

Parameter	Description
$t_i^{l-1}$	第$l-1$层的第$i$个神经元的点火时刻（spike）
$\omega_{ij}^l$	连接神经元$N_i^{l-1}$和神经元$N_j^{l}$的突触权重
$\epsilon(\cdot)$	PSP（突触后电位）函数,以$\alpha-$函数最为常用
$\epsilon(t-t_i^{l-1})$	在$t_i^{l-1}$时刻产生的PSP
$t_j^l=\mathcal{F}\{t|V_j^l(t)=\vartheta,t\geq0\}$	尖峰发生函数，在膜电位等于$\vartheta$的时刻产生脉冲的时间
$\eta(t-t_j^l)$	膜电位恢复函数（refractory kernel）

根据上面对SNN网络的建模，我们可以通过BP算法更新参数：

$$\begin{aligned}\frac{\partial t_j^l}{\partial\omega_{ij}^l}&=\frac{\partial t_j^l}{\partial V_j^l\big(t_j^l\big)}\frac{\partial V_j^l\big(t_j^l\big)}{\partial\omega_{ij}^l},&\mathrm{if~}t_j^l>t_i^{l-1}\\\frac{\partial t_j^l}{\partial t_i^{l-1}}&=\frac{\partial t_j^l}{\partial V_j^l\big(t_j^l\big)}\frac{\partial V_j^l\big(t_j^l\big)}{\partial t_i^{l-1}},&\mathrm{if~}t_j^l>t_i^{l-1}.\end{aligned}$$

计算$\frac{\partial t_j^l}{\partial V_j^l}(t_j^l)$, 假设膜电位$V_j^l(t)$在脉冲产生时间$t_j$的邻域内认为线性增加：

$$\frac{\partial t_j^l}{\partial V_j^l(t_j^l)}=\frac{-1}{\partial V_j^l(t_j^l)/\partial t_j^l}=\frac{-1}{\sum_i^N\omega_{ij}^l\frac{\partial \epsilon\left(t_j^l-t_i^{l-1}\right)}{\partial t_j^l}}$$

$$with\space\frac{\partial\varepsilon\left(t_j^l-t_i^{l-1}\right)}{\partial t_j^l}=\frac{\exp\left(1-\left(t_j^l-t_i^{l-1}\right)/\tau\right)}{\tau^2}\big(\tau+t_i^{l-1}-t_j^l\big).$$

问题一：由于脉冲函数的离散性，脉冲函数不可导

问题二：梯度爆炸发生在$\frac{\partial V_j^l(t_j^l)}{\partial t_j}$,比如说即将到达阈值发射脉冲的时刻（如下图所示）

问题三：坏死的神经元（dead neuron）：如果突触前电位不足以使神经元点火，即探察不到有$t_j$点火时刻的出现，就更无从谈起计算误差了。误差反向传播不了。PSP尖峰产生机制具有泄漏特性（？），坏死现象更容易产生（???）

三个问题怎么克服？

1. ReL-PSP-Based Spiking Neuron Model

索性重新建模神经元：

$$new: \space V_j^l(t) = \sum_{i}^{N}\omega_{ij}^lK(t-t_i^{l-1})$$

$$old: \space V_j^l(t)=\sum_i^N\omega_{ij}^l\varepsilon\left(t-t_i^{l-1}\right)-\eta(t-t_j^l)$$

对比一下，不难发现，不仅舍去了恢复函数$\eta(\cdot)$，还更改了核函数$\epsilon(\cdot)$为$K(\cdot)$.

核函数K的定义如下：

$$K\left(t-t_i^{l-1}\right)=\begin{cases}t-t_i^{l-1},&\mathrm{if~}t>t_i^{l-1}\\0,&\text{otherwise.}&\end{cases}$$

眼熟吗？很像一位老朋友：ReLU

$$ReLU(x) = \begin{cases}x,&\mathrm{if~}x \gt 0\\0,&\text{otherwise.}&\end{cases}$$

不能说极度相似，简直是一模一样。（难怪作者起这个名字）

ReL-PSP的线性特性，使得膜电位在尖峰时刻达到之前（prior to）线性增长. 计算$\frac{\partial t_j^l}{\partial V_j^l(t_j^l)}$也变得很简单：

$$\begin{aligned} & new: \space \\ \frac{\partial t_j^l}{\partial V_j^l(t_j^l)}& =-\frac1{\partial V_j^l(t_j^l)/\partial t_j^l} \\ &=\frac{-1}{\sum_i^N\omega_{ij}^l\frac{\partial K\left(t_j^l-t_i^{l-1}\right)}{\partial t_j^l}} \\ &=\frac{-1}{\sum_i^N\omega_{ij}^l},\quad\mathrm{if~}t_j^l>t_i^{l-1}. \end{aligned}$$

$$old: \space \frac{\partial t_j^l}{\partial V_j^l(t_j^l)}=\frac{-1}{\partial V_j^l(t_j^l)/\partial t_j^l}=\frac{-1}{\sum_i^N\omega_{ij}^l\frac{\partial \epsilon\left(t_j^l-t_i^{l-1}\right)}{\partial t_j^l}}$$

有了ReL-PSP核函数，我们甚至不用假设近似线性。逐层的累计误差相当于没有了。

梯度爆炸问题中，由于分母中的$\omega_{ij} ^l$依然有可能接近于零，这个问题没有得到完全解决。但由于$t_j^l$是输入脉冲和突触权重有关的函数，可以根据以下公式计算：

$$t_j^l=\frac{\vartheta+\sum_i^N\omega_{ij}^lt_i^{l-1}}{\sum_i^N\omega_{ij}^l},\quad\mathrm{if~}t_j^l>t_i^{l-1}.$$

即便是$\sum_{i}^{N}\omega_{ij}^l$接近于0，对膜电位影响很小，产生点火的时刻就可能会更晚，表征为$t_j^l$更大，且也可能对下一层的点火贡献很小或者没有。因此神经元$N_j^i$在上述情况下不参与误差反向传播，也不会引起梯度爆炸问题。

（好的这段话已经把我绕晕了）

坏死神经元问题：ReL-PSP摒弃了恢复函数，即膜电位不会恢复，只会不断积累。也就是说无论如何神经元都会发出脉冲，只是时间早晚的问题。

Error-BP

问题背景：n分类问题。

损失函数：交叉熵

$$L(g,\mathbf{t^o})=-\ln\frac{\exp(-\mathbf{t^o}[g]))}{\sum_i^n\exp(-\mathbf{t^o}[i])}$$

输出层函数：softmax

$$p_j=\frac{\exp(-t_j)}{\sum_i^n\exp(-t_i)}$$

参数解释：

Parameter	Description
$ \bold {t^o}$	输出层脉冲发射时刻的向量
$g$	目标类别索引

反向传播的偏导计算：

$$\begin{aligned}\frac{\partial t_j^l}{\partial\omega_{ij}^l}&=\frac{\partial t_j^l}{\partial V_j^l(t_j^l)}\frac{\partial V_j^l(t_j^l)}{\partial\omega_{ij}^l}=\frac{t_i^{l-1}-t_j^l}{\sum_i^N\omega_{ij}^l},&\text{if }t_j^l>t_i^{l-1}\\\frac{\partial t_j^l}{\partial t_i^{l-1}}&=\frac{\partial t_j^l}{\partial V_j^l(t_j^l)}\frac{\partial V_j^l(t_j^l)}{\partial t_i^{l-1}}=\frac{\omega_{ij}^l}{\sum_i^N\omega_{ij}^l},&\text{if }t_j^l>t_i^{l-1}.\end{aligned}$$

遗留问题：

1. 为什么反向传播更新的过程中要更新$t_i^{l-1}$?

2. 关于梯度消失推导的最后一段描述，和坏死神经元的关系？

3. OpenSMART NoC?

经典算法STBP：

SNN训练的三种方式：

1. 无监督学习：STDP

2. 间接监督学习：ANN-to-SNN

3. 直接监督学习：梯度下降

尖峰神经网络中的迭代泄漏整合与发射模型（Iterative Leaky Integrate-And-Fire （LIF） Model in Spiking Neural Networks）

LIF用于对SNN网络中神经元行为进行描述：

参数解释：

Parameter	Description
$u(t)$	t时刻的神经元膜电位
$\tau$	时间常数
$I(t)$	突触前活动状态，外部刺激和突触权重
$V_{th}$	膜电位阈值

$$\tau\frac{du(t)}{dt}=-u(t)+I(t)$$

直观感受这个公式，它在描述膜电位随时间的变化状态。变化快慢由$\tau$调控，外部没有刺激的时候，膜电位随时间在不断减少，表现为$u(t)$存在负号；外部刺激或者突触前馈对膜电位的变化起到的是正向的激励作用($I(t)$为正号)

为了解决SNN网络在时间域上动态变化复杂，不便于进行反向传播，首先对上述的线性微分方程在$u(t) |_{t=t_i-1}$的初始条件下求解，得到如下的迭代规律：

$$u(t) = u(t_{i-1})e^{\frac{t_{i-1}-t}{\tau}} + \hat{I}(t)$$

即：t时刻的膜电位与前一时刻的膜电位（Temporal Dynamics, TD）以及t时刻突触前输入$\hat{I}(t)$（Spatial Dynamics, SD）有关。

通过DNN中误差反向传播的训练方法，研究者从中获得启发，DNN的误差反向传播是在空间域上进行的，如何拓展到基于LIF模型的SNN上呢？

$$\begin{aligned} &x_i^{t+1,n} =\sum_{j=1}^{l(n-1)} w_{ij}^no_j^{t+1,n-1} \\ &u_{i}^{t+1,n} =u_i^{t,n}f(o_i^{t,n})+x_i^{t+1,n}+b_i^n \\ &o_{i}^{t+1,n} =g(u_i^{t+1,n}) \end{aligned}$$

参数解释：

Parameter	Description
$f(x) = \tau e^{-\frac{x}{\tau}}$	forget gate(from LSTM ,我去太机智了 )，用于控制TD记忆泄露程度
$g(x)=\begin{cases}1,&x\ge V_{th}\\0,&x<V_{th}\end{cases}$	output gate(from LSTM)，用于控制脉冲释放
上标t	t时刻
$n$	第n层
$l(n)$	（第n层中）第l个神经元
$w_{ij}$	突触前层第 j 个神经元到突触后层第 i 个神经元的突触权重
$o_j \in \{0,1\}$	神经元输出脉冲与否，1代表输出，0代表不输出
$u_i^{t+1,n}$	第n层第i个神经元在t+1时刻的膜电位
$x_i$	$\hat{I}(t)$的简化表示，即第i个神经元的突触前输入
$b_i$	表征点火阈值的变量, 利用可调偏置 b 来模拟阈值行为

对于一个很小的时间之内（又要进行必要的假设了），$f(\cdot)$可以近似认为：

$$f(o_i^{t,n})\approx\begin{cases}\tau,&o_i^{t,n}=0\\0,&o_i^{t,n}=1\end{cases}$$

$\tau e^{-\frac{1}{\tau}} \approx 0$

(额？)

STBP训练框架

定义损失函数为均方误差累积：

$$L=\frac{1}{2S}\sum_{s=1}^{S}\parallel y_s-\frac{1}{T}\sum_{t=1}^{T}\boldsymbol{o}_s^{t,N}\parallel_2^2$$

参数解释：

Parameter	Description
$y_s$	第s个训练样本标签
$o_s$	第s个训练样本在特定时刻的输出

损失函数$L$是一个关于权重$\boldsymbol{W}$和偏置$\boldsymbol{b}$有关的函数式。即每次需计算$\frac{\partial L}{\partial \boldsymbol{W}}$ 和 $\frac{\partial L}{\partial \boldsymbol{b}}$并迭代更新。

假设我们已经得到了每一层n在特定时刻t的输出偏导数：$\frac{\partial L}{\partial o_i}$和$\frac{\partial L}{\partial u_i}$

对于单个神经元，传播被分解为垂直路径 SD 和水平路径 TD。SD 中的误差传播数据流类似于 DNN 的典型 BP，即每个神经元累积来自上层的加权误差信号，并迭代更新不同层的参数。而 TD 中的数据流共享相同的神经元状态，这使得直接获得解析解变得相当复杂。为了解决这个问题，研究者利用提出的迭代 LIF 模型在 SD 和 TD 方向上展开状态空间，从而可以区分不同时间步长的 TD 中的状态，这就实现了迭代传播的链式规则。Werbos (1990) 中用于训练 RNN 的 BPTT 算法也有类似的思想。

denote:

$$\delta_i^{t,n} = \frac{\partial L}{\partial o_i^{t,n}}$$

Case 1: $t=T$ at output layer $n=N$

$$\frac{\partial L}{\partial o_i^{T,N}}=-\frac1{TS}(y_i-\frac1T\sum_{k=1}^To_i^{k,N}) = \delta_i^{T,N}.$$

$$\frac{\partial L}{\partial u_i^{T,N}}=\frac{\partial L}{\partial o_i^{T,N}}\frac{\partial o_i^{T,N}}{\partial u_i^{T,N}}=\delta_i^{T,N}\frac{\partial o_i^{T,N}}{\partial u_i^{T,N}}.$$

Case 2: $t=T$ at output layer $n \lt N$ (非输出层，隐藏层中)

$$\frac{\partial L}{\partial o_i^{T,n}}=\sum_{j=1}^{l(n+1)}\delta_j^{T,n+1}\frac{\partial o_j^{T,n+1}}{\partial o_i^{T,n}}=\sum_{j=1}^{l(n+1)}\delta_j^{T,n+1}\frac{\partial g}{\partial u_i^{T,n}}w_{ji}.$$

(本人对上式深表怀疑)

$$\frac{\partial L}{\partial u_i^{T,n}}=\frac{\partial L}{\partial o_i^{T,n}}\frac{\partial o_i^{T,n}}{\partial u_i^{T,n}}=\delta_i^{T,n}\frac{\partial g}{\partial u_i^{T,n}}.$$

Case 3: $t \lt T$ at output layer $n=N$

$$\begin{aligned} \frac{\partial L}{\partial o_i^{t,N}}& =\delta_i^{t+1,N}\frac{\partial o_i^{t+1,N}}{\partial o_i^{t,N}}+\frac{\partial L}{\partial o_i^{T,N}} \\ &=\delta_i^{t+1,N}\frac{\partial g}{\partial u_i^{t+1,N}}u_i^{t,N}\frac{\partial f}{\partial o_j^{t,N}}+\frac{\partial L}{\partial o_i^{T,N}}, \end{aligned}$$

$$\frac{\partial L}{\partial u_i^{t,N}}=\frac{\partial L}{\partial u_i^{t+1,N}}\frac{\partial u_i^{t+1,N}}{\partial u_i^{t,N}}=\delta_i^{t+1,N}\frac{\partial g}{\partial u_i^{t+1,N}}f(o_i^{t,n}),$$

where $\frac{\partial L}{\partial o_{i}^{T,N}}=-\frac{1}{TS}(y_{i}-\frac{1}{T}\sum_{k=1}^{T}o_{i}^{k,N})$

Case 4: $t \lt T$ at output layer $n \lt N$

$$\begin{aligned} &\frac{\partial L}{\partial o_i^{t,n}} =\sum_{j=1}^{l(n+1)}\delta_j^{t,n+1}\frac{\partial o_j^{t,n+1}}{\partial o_i^{t,n}}+\frac{\partial L}{\partial o_i^{t+1,n}}\frac{\partial o_i^{t+1,n}}{\partial o_i^{t,n}} \\ &=\sum_{j=1}^{l(n+1)}\delta_j^{t,n+1}\frac{\partial g}{\partial u_i^{t,n}}w_{ji}+\delta_i^{t+1,n}\frac{\partial g}{\partial u_i^{t,n}}u_i^{t,n}\frac{\partial f}{\partial o_i^{t,n}}, \\ &\frac{\partial L}{\partial u_i^{t,n}} =\frac{\partial L}{\partial o_i^{t,n}}\frac{\partial o_i^{t,n}}{\partial u_i^{t,n}}+\frac{\partial L}{\partial o_i^{t+1,n}}\frac{\partial o_i^{t+1,n}}{\partial u_i^{t,n}} \\ &=\delta_i^{t,n}\frac{\partial g}{\partial u_i^{t,n}}+\delta_i^{t+1,n}\frac{\partial g}{\partial u_i^{t+1,n}}f(o_i^{t,n}). \end{aligned}$$

最终得到关于$\boldsymbol{W,b}$的偏导数表达式：

$$\frac{\partial L}{\partial\boldsymbol{b}^n}=\sum_{t=1}^T\frac{\partial L}{\partial\boldsymbol{u}^{t,n}}\frac{\partial\boldsymbol{u}^{t,n}}{\partial \boldsymbol{b}^n}=\sum_{t=1}^T\frac{\partial L}{\partial\boldsymbol{u}^{t,n}},$$

$$\begin{aligned} \frac{\partial L}{\partial W^n}& =\sum_{t=1}^T\frac{\partial L}{\partial u^{t,n}}\frac{\partial u^{t,n}}{\partial W^n} \\ &=\sum_{t=1}^T\frac{\partial L}{\partial u^{t,n}}\frac{\partial u^{t,n}}{\partial x^{t,n}}\frac{\partial x^{t,n}}{\partial W^n}=\sum_{t=1}^T\frac{\partial L}{\partial u^{t,n}}o^{t,n-1},T, \end{aligned}$$

(上面这个式子应该是错误的)

对尖峰函数的近似表达：

$$\begin{aligned} &h_{1}(u) =\frac1{a_1}sign(|u-V_{th}|<\frac{a_1}2), \\ &h_2(u) =(\frac{\sqrt{a_2}}2-\frac{a_2}4|u-V_{th}|)sign(\frac2{\sqrt{a_2}}-|u-V_{th}|), \\ &h_3(u) =\frac1{a_3}\frac{e^{\frac{V_{th}-u}{a_3}}}{(1+e^{\frac{V_{th}-u}{a_3}})^2}, \\ &h_4(\mathcal{U}) =\frac1{\sqrt{2\pi\textit{a}_4}}e^{-\frac{(u-V_{th})^2}{2a_4}}, \end{aligned}$$

STBP代码复现1：MNIST分类任务

ref: https://github.com/yjwu17/STBP-for-training-SpikingNN

a. spiking_model.py

1. 包依赖

import torch
import torch.nn as nn
import torch.nn.functional as F

2. 基本参数设定

device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
thresh = 0.5  # 阈值
lens = 0.5  # 冲激函数近似值的超参数
decay = 0.2  # 衰减常数(decay constants)
num_classes = 10
batch_size = 100
learning_rate = 1e-3
num_epochs = 100  # max epoch

3. 近似点火函数

class ActFun(torch.autograd.Function):

    @staticmethod
    def forward(ctx, input):
        ctx.save_for_backward(input)
        return input.gt(thresh).float()

    @staticmethod
    def backward(ctx, grad_output):
        input, = ctx.saved_tensors
        grad_input = grad_output.clone()
        temp = abs(input - thresh) < lens
        return grad_input * temp.float()

该部分重写了两个方法：前向传播和反向传播。首先，函数传进torch.autograd.Function，forward 方法接受输入张量 input，将其与预定义的 thresh 进行比较，然后通过 input.gt(thresh).float() 返回一个新的张量，其中每个元素都是对应输入张量元素是否大于阈值的浮点值。ctx.save_for_backward(input) 用于保存输入张量，以便在反向传播时使用。

在反向传播中，backward 方法接受梯度输出 grad_output，然后通过 ctx.saved_tensors 获取保存的输入张量。梯度计算采用了一些条件判断，具体来说，计算了绝对值小于 lens 的部分。最终，返回的梯度是 grad_input * temp.float()，其中 temp 是一个布尔张量，指示输入张量元素是否满足条件。

最后，通过 act_fun = ActFun.apply 创建了一个可以在模型中使用的激活函数。使用时，可以通过调用 act_fun(input) 进行前向传播。

4. 膜电位更新函数

def mem_update(ops, x, mem, spike):
    mem = mem * decay * (1. - spike) + ops(x)
    spike = act_fun(mem)  # act_fun : approximation firing function
    return mem, spike

膜电位 * 衰减系数 * (1 - 上一次的脉冲) + 对x运算过后的结果（$ops(x)$）

当前膜电位决定是否点火：
$act \_ fun(men)$

5. CNN网络层配置参数

# cnn_layer(in_planes, out_planes, stride, padding, kernel_size)
cfg_cnn = [(1, 32, 1, 1, 3),
           (32, 32, 1, 1, 3), ]
# kernel size
cfg_kernel = [28, 14, 7]
# fc layer
cfg_fc = [128, 10]

参数解释：

Parameters	Description
in_planes	输入通道数
out_planes	输出通道数
stride	步长
padding	填充
kernel_size	卷积核大小

第一层：输入通道数为 1，输出通道数为 32，步幅为 1，填充为 1，卷积核大小为 3。

第二层：输入通道数为 32，输出通道数为 32，步幅为 1，填充为 1，卷积核大小为 3。

cfg_kernel: 整数列表，存有不同层次的卷积核大小；

cfg_fc: 整数列表，表示全连接层的配置，一边是128个神经元，一边是10个神经元。

6. 学习率调度器函数

# Dacay learning_rate
def lr_scheduler(optimizer, epoch, init_lr=0.1, lr_decay_epoch=50):
    """Decay learning rate by a factor of 0.1 every lr_decay_epoch epochs."""
    if epoch % lr_decay_epoch == 0 and epoch > 1:
        for param_group in optimizer.param_groups:
            param_group['lr'] = param_group['lr'] * 0.1
    return optimizer

在每个 lr_decay_epoch 周期后将学习率衰减为原来的 0.1 倍。

optimizer.param_groups是一个包含参数组的列表，每个参数组是一个字典，包含一组参数和这组参数对应的优化选项，如学习率、权重衰减等。

7. SCNN网络架构设计

class SCNN(nn.Module):
    def __init__(self):
        super(SCNN, self).__init__()
        in_planes, out_planes, stride, padding, kernel_size = cfg_cnn[0]
        self.conv1 = nn.Conv2d(in_planes, out_planes, kernel_size=kernel_size, stride=stride, padding=padding)
        in_planes, out_planes, stride, padding, kernel_size = cfg_cnn[1]
        self.conv2 = nn.Conv2d(in_planes, out_planes, kernel_size=kernel_size, stride=stride, padding=padding)

        self.fc1 = nn.Linear(cfg_kernel[-1] * cfg_kernel[-1] * cfg_cnn[-1][1], cfg_fc[0])
        self.fc2 = nn.Linear(cfg_fc[0], cfg_fc[1])

    def forward(self, input, time_window=20):
        # 初始化膜电位和激活脉冲
        c1_mem = c1_spike = torch.zeros(batch_size, cfg_cnn[0][1], cfg_kernel[0], cfg_kernel[0], device=device)
        # 存储第一个卷积层的膜电位和脉冲的张量,它们的形状应该与第一个卷积层的输出形状相同
        c2_mem = c2_spike = torch.zeros(batch_size, cfg_cnn[1][1], cfg_kernel[1], cfg_kernel[1], device=device)
        # 同理
        h1_mem = h1_spike = h1_sumspike = torch.zeros(batch_size, cfg_fc[0], device=device)
        h2_mem = h2_spike = h2_sumspike = torch.zeros(batch_size, cfg_fc[1], device=device)

        for step in range(time_window):  # 模拟时间步长
            x = input > torch.rand(input.size(), device=device)  # 概率发放，x为布尔张量，每个元素都是由与对应位置上 input 元素相比较的随机数是否大于该位置的 input 元素得到的，用于模拟脉冲神经网络中的神经元的随机脉冲发放行为

            c1_mem, c1_spike = mem_update(self.conv1, x.float(), c1_mem, c1_spike) # 第一层对x卷积之后，根据输出和前一时刻的膜电位更新电位状态，以及确定要不要点火。

            x = F.avg_pool2d(c1_spike, 2)
            # c1_spike作为卷积层的输出，进入avg_pool过程，池化核大小是2，等价于二倍降采样。

            c2_mem, c2_spike = mem_update(self.conv2, x, c2_mem, c2_spike)

            x = F.avg_pool2d(c2_spike, 2)
            x = x.view(batch_size, -1)
            # view函数改变张量形状，第二个参数“-1”代表自动计算输入张量大小，相当于展成一个大的向量。这种操作通常在卷积层和全连接层之间进行，因为全连接层期望的输入是一维的，而卷积层的输出是多维的。

            h1_mem, h1_spike = mem_update(self.fc1, x, h1_mem, h1_spike)
            h1_sumspike += h1_spike
            h2_mem, h2_spike = mem_update(self.fc2, h1_spike, h2_mem, h2_spike)# ？为什么不是h1_sumspike?(mem_update 函数可能是用来更新神经元的状态，包括膜电位和脉冲。这个更新过程通常只依赖于当前时间步的输入和当前的膜电位，而不需要知道过去的脉冲总和。h1_sumspike 可能用于其他目的，例如在训练过程中监控脉冲的总数，或者在网络的其他部分使用。)
            h2_sumspike += h2_spike 

        outputs = h2_sumspike / time_window
        return outputs

b. main.py

1. 包依赖

from __future__ import print_function
import torchvision
import torchvision.transforms as transforms
import os
import time
from spiking_model import *

2. 设定训练任务的基本配置

names = 'spiking_model'
data_path = './MNIST/'  # todo: input your data path
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
train_dataset = torchvision.datasets.MNIST(root=data_path, train=True, download=True, transform=transforms.ToTensor())
train_loader = torch.utils.data.DataLoader(train_dataset, batch_size=batch_size, shuffle=True, num_workers=0)

test_set = torchvision.datasets.MNIST(root=data_path, train=False, download=True, transform=transforms.ToTensor())
test_loader = torch.utils.data.DataLoader(test_set, batch_size=batch_size, shuffle=False, num_workers=0)

best_acc = 0  # best test accuracy
start_epoch = 0  # start from epoch 0 or last checkpoint epoch
acc_record = list([])
loss_train_record = list([])
loss_test_record = list([])

snn = SCNN()
snn.to(device)
criterion = nn.MSELoss()
optimizer = torch.optim.Adam(snn.parameters(), lr=learning_rate)

3. SCNN训练循环

for epoch in range(num_epochs):
    running_loss = 0 # 初始化损失为 0
    start_time = time.time() # 记录开始时间。
    for i, (images, labels) in enumerate(train_loader): # 在训练数据加载器（train_loader）中的每一批数据（images 和 labels）进行循环
        snn.zero_grad() # 清零神经网络（snn）的梯度
        optimizer.zero_grad() # 清零优化器（optimizer）的梯度
        # 调用 .backward() 方法计算梯度时，梯度会累积在叶子节点的 .grad 属性中，而不是被替换。这意味着每次你计算梯度，新的梯度值会被添加到已有的梯度值上，而不是替换它。因此，如果你在一个训练循环中多次计算梯度，可能会得到意外的结果，因为梯度是从多个反向传播过程中累积起来的。optimizer.zero_grad() 就足够了，因为它会清零所有模型参数的梯度，这包括 snn 的梯度。snn.zero_grad() 可能是多余的，除非 snn 中有一些参数没有被 optimizer 管理。（事实证明，snn.zero_grad()确实多余）

        images = images.float().to(device) # 将图像数据转换为浮点数并移动到设备（device）上。
        outputs = snn(images) # outputs = h2_sumspike / time_window
        labels_ = torch.zeros(batch_size, 10).scatter_(1, labels.view(-1, 1), 1) # 标签转为独热码
        loss = criterion(outputs.cpu(), labels_)
        running_loss += loss.item() # 累加损失到运行损失（running_loss）中
        loss.backward() # 反向传播损失
        optimizer.step() # 使用优化器更新神经网络的参数。
        if (i + 1) % 100 == 0:
            print('Epoch [%d/%d], Step [%d/%d], Loss: %.5f'
                  % (epoch + 1, num_epochs, i + 1, len(train_dataset) // batch_size, running_loss))
            running_loss = 0
            print('Time elapsed:', time.time() - start_time)
    correct = 0
    total = 0
    optimizer = lr_scheduler(optimizer, epoch, learning_rate, 40)

4. 测试阶段：

with torch.no_grad(): # 在测试模型时，我们不需要计算梯度，因此可以关闭梯度计算以节省内存和计算资源。
    for batch_idx, (inputs, targets) in enumerate(test_loader):
        inputs = inputs.to(device)
        optimizer.zero_grad()
        outputs = snn(inputs)
        # outputs 是神经网络的输出，它是一个形状为 [batch_size, num_classes] 的张量，其中 batch_size 是批次大小，num_classes 是类别数量。每一行都包含了一个样本对每个类别的预测分数。
        labels_ = torch.zeros(batch_size, 10).scatter_(1, targets.view(-1, 1), 1)
        loss = criterion(outputs.cpu(), labels_)
        _, predicted = outputs.cpu().max(1) 
        # .max(1) 是在这个张量的第二个维度（索引从 0 开始）上找最大值。这会返回两个张量：一个包含每行的最大值，一个包含最大值的索引。在分类问题中，最大值的索引就是模型预测的类别。
        # .cpu() 将张量从 GPU 移动到 CPU。这通常在你需要使用 numpy 或者打印张量的值时需要，因为这些操作只能在 CPU 上进行。
        total += float(targets.size(0)) # 得到目前为止处理过的所有样本的数量。
        correct += float(predicted.eq(targets).sum().item())
        # 比较预测的类别（predicted）和真实的标签（targets）是否相等。这会返回一个布尔张量，其中 True 表示预测正确，False 表示预测错误。
        # .sum() 是将所有的 True（被视为 1）加起来，得到预测正确的样本数。
        # .item() 是将一个只有一个元素的张量转换为一个 Python 数字。
        if batch_idx % 100 == 0:
            acc = 100. * float(correct) / float(total)
            print(batch_idx, len(test_loader), ' Acc: %.5f' % acc)

5. 输出其他内容和保存模型状态

print('Iters:', epoch, '\n\n\n')
print('Test Accuracy of the model on the 10000 test images: %.3f' % (100 * correct / total))
acc = 100. * float(correct) / float(total)
acc_record.append(acc)
if epoch % 5 == 0:
    print(acc)
    print('Saving..')
    state = {
        'net': snn.state_dict(),
        'acc': acc,
        'epoch': epoch,
        'acc_record': acc_record,
    }
    if not os.path.isdir('checkpoint'):
        os.mkdir('checkpoint')
    torch.save(state, './checkpoint/ckpt' + names + '.t7')
    best_acc = acc

STBP代码复现2：CIFAR-10分类任务

ref:

1. Deep Residual Learning in Spiking Neural Networks, Wei Fang et al., Neural Information Processing Systems (NeurIPS) 2021 (SEW/Spike-Element-Wise block)

2. Differentiable Spike: Rethinking Gradient-Descent for Training Spiking Neural Networks, Yuhang Li et al., Neural Information Processing Systems (NeurIPS 2021). (a new family of Differentiable Spike (Dspike))

1. SEW block的结构：

神经元行为建模：

$$\begin{aligned} &H[t] =f(V[t-1],X[t]), \\ &S[t] =\Theta(H[t]-V_{th}), \\ &V[t] =H[t](1-S[t])+V_{reset}S[t], \end{aligned}$$

参数解释：

Parameter	Description
$H[t]$	神经元在进行运算活动之后的膜电位
$f(\cdot)$	神经元动态活动（运算方式），可以是卷积，可以是其他形式
$V[t]$	t时刻发出脉冲之后的膜电位
$\Theta(\cdot)$	阶跃函数
$V_{th}$	点火阈值
$V_{reset}$	膜电位恢复函数
$S[t]$	t时刻发射的脉冲

例子：两个神经元行为的建模方式：IF（Integrate-and-Fire，式一）和LIF（Leaky-Integrate-and-Fire，式二）

$$\begin{aligned}&H[t]=V[t-1]+X[t],\\&H[t]=V[t-1]+\frac{1}{\tau}(X[t]-(V[t-1]-V_{reset})),\end{aligned}$$

Identity Mapping(恒等映射): ResNet的核心思想

当$\mathcal{F}^l(X^l)\equiv0,Y^l=\mathrm{ReLU(X^l)}$，多数情况下$X^l$是上一个ReLU层的输出（非负）,所以$Y^l = \mathrm{ReLU(X^l)} = X^l$

但对于从ANN迁移到SNN的Spiking Resnet块，就无法满足Identity Mapping：

$$\begin{aligned}Y^l&=\mathrm{ReLU}(\mathcal{F}^l(X^l)+X^l),\\O^l[t]&=\mathrm{SN}(\mathcal{F}^l(S^l[t])+S^l[t]).\end{aligned}$$

当 $\mathcal{F}^l(S^l[t])\equiv0,O^{l}[t]=\mathrm{SN}(S^{l}[t])\neq S^{l}[t].$。满足 $\mathrm{SN}(S^{l}[t]) = S^{l}[t]$ 的条件是后一个神经元在t时刻也发出脉冲，在 t 时刻没有接收到尖峰脉冲后保持沉默。

mad先不读了，spikingjelly是真好用，文档可视化做的真的不错，非常适合入门

Continue:

对于IF神经元，恒等映射可以满足。设定 $0 \lt V_{th} \leq 1$, $V[t-1] = 0$。进而保证当$X[t] = 1$时, $H[t] \gt V_{th}$；$X[t] = 0$时, $H[t] \lt V_{th}$.

但对于更复杂的神经元行为建模方法，很难保证这一条件成立。例如LIF神经元，存在一个电位恢复常数$\tau$。它作为一个超参参与信息传递，当设定 $0 \lt V_{th} \leq 1$, $V[t-1] = 0$。$X[t] = 1$时, $H[t] = \frac{1}{\tau}$.此时很难确定$V_{th}$的值。

Vanishing & Exploding Gradient in Spiking Resnet blocks

当在t时刻，有k个Spiking Resnet块参与神经网络的前馈传递$S^l[t]$时，且假设满足恒等映射，有

$$S^l[t]=S^{l+1}[t]=...=S^{l+k-1}[t]=O^{l+k-1}[t].$$

输出对第l层的输入求偏导进行反向传播时有：

$$\frac{\partial O_j^{l+k-1}[t]}{\partial S_j^l[t]}=\prod_{i=0}^{k-1}\frac{\partial O_j^{l+i}[t]}{\partial S_j^{l+i}[t]}=\prod_{i=0}^{k-1}\Theta'(S_j^{l+i}[t]-V_{th})\to\begin{cases}0,\mathrm{if}0<\Theta'(S_j^l[t]-V_{th})<1\\1,\mathrm{if}\Theta'(S_j^l[t]-V_{th})=1\\+\infty,\mathrm{if}\Theta'(S_j^l[t]-V_{th})>1\end{cases}$$

式中，$\Theta(x)$是阶跃函数，$\Theta'(x)$是替代梯度函数。

逐元素点火的Resnet块（Spike-Element-Wise Resnet）

$$O^l[t]=g(\mathrm{SN}(\mathcal{F}^l(S^l[t])),S^l[t])=g(A^l[t],S^l[t]),$$

$g(\cdot)$是一个逐元素函数，输入是代表两个冲激的张量。$A^l[t] = \mathrm{SN}(\mathcal{F}^l(S^l[t]))$.

SEW-Resnet的优点：

1. 更易于实现恒等映射：利用尖峰的二进制属性，我们可以找到满足恒等映射的不同元素函数$g(\cdot)$

Name	Expression of $g(\cdot)$
ADD(加)	${A^{l}[t]+S^{l}[t]}$
AND(与)	$\begin{aligned}A^l[t]\wedge S^l[t]=A^l[t]\cdot S^l[t]\end{aligned}$
IAND	$(\neg A^l[t])\wedge S^l[t]=(1-A^l[t])\cdot S^l[t]$

e.g.: 选择ADD函数或者IAND函数作为Element-Wise Function，设定$A^l[t] = 0$，即$\mathrm{SN}(\mathcal{F}^l(S^l[t])) = 0$, 最后一个BN层的所有权重和直流偏置都置为0，使得主干路上不能使神经元点火。

相反地，选择AND函数，设定$A^l[t] = 1$，最后一层的权重和偏置足够大，使得后一个神经元绝对可以产生脉冲。但是不能时时刻刻都保证权重足够大，毕竟梯度爆炸问题不容易发生。

拓展：SEW-Resnet的下采样块设计：用于输入输出维度不一致的时候。Conv层的步长大于1

SEW-Resnet与ReLU-before-Addition（RBA）：

某种程度上，SEW的行为是RBA的行为拓展，但不一样的是，RBA会导致输出层不断积累，而SEW中，使用AND和IAND作为g会输出尖峰(即二元张量)，这意味着ANNs中的无限输出问题不会发生在具有SEW块的SNNs中，因为所有尖峰都小于或等于1。当选取ADD为g时，由于k个顺序SEW块的输出不会大于k + 1，可以缓解无穷输出问题。此外，当g为ADD时，一个下采样的SEW块将调节输出不大于2。

SEW对于梯度爆炸/消失问题的缓解

逐层计算( l + k - 1) - th SEW块的输出相对于第l - th SEW块输入的梯度：

$$\frac{\partial O_j^{l+k-1}[t]}{\partial S_j^l[t]}=\prod_{i=0}^{k-1}\frac{\partial g(A_j^{l+i}[t],S_j^{l+i}[t])}{\partial S_j^{l+i}[t]}=\begin{cases}\prod_{i=0}^{k-1}\frac{\partial(0+S_j^{l+i}[t])}{\partial S_j^{l+i}[t]}\text{,if }g=ADD\\\prod_{i=0}^{k-1}\frac{\partial(1\cdot S_j^{l+i}[t])}{\partial S_j^{l+i}[t]}\text{,if }g=AND\\\prod_{i=0}^{k-1}\frac{\partial((1-0)\cdot S_j^{l+i}[t])}{\partial S_j^{l+i}[t]}\text{,if }g=IAND\end{cases}=1.$$

补习填坑：复变函数