Matlab实现脉冲神经网络SNN：LIF神经元、STDP学习与Polychronization聚群检测

简介面向SNN脉冲神经网络研究者的Matlab仿真资源基于Matlab 2021a构建覆盖脉冲神经元网络构建、STDP学习规则与多锋电位polychronization仿真等关键环节适合需要快速搭建第三代神经网络实验环境的本科生、研究生或工程师也可作为神经形态计算课程的教学演示素材。资源包共15个文件包含8个.m源码、1个.mat数据文件、5张结果图以及1段avi操作录像压缩包约14.7MB源码与演示数据分离结构清晰便于按模块研究。已有1180人学习下载录像使用Windows Media Player播放详细演示了MATLAB当前文件夹路径设置等易错细节可显著缩短环境配置时间。通过Runme1.m、Runme2.m两个入口可快速启动仿真配合createNetwork、createNeuron、STDP、polychron等多个函数能系统观察脉冲序列的生成、传导与学习过程附带的jpg结果图可直接对比输出特征为复现实验和扩展算法提供了可直接修改的完整工程基础。1. SNN脉冲神经网络为什么值得现在回头看当深度学习社区把注意力集中在 Transformer 和扩散模型上时脉冲神经网络SNN这个“第三代神经网络”反倒在边缘场景里活了下来。它和ANN最本质的区别不是多了一个时间维度而是每个神经元只在膜电位跨过阈值时发出一瞬脉冲——信息分布在“发不发”和“什么时候发”里而不是浮点激活值里。用matlab做SNN仿真有一层特殊价值matlab的矩阵化写法天然适合把突触权重、延迟矩阵和脉冲栅格图组织起来调试门槛比C低得多。这套资源包含从神经元构建、网络连接到STDP学习规则的完整matlab工程还附带一段操作录像覆盖了从Runme1.m跑到观察polychronization聚群现象的全过程。它适合两类人一是课程设计或毕业论文需要SNN实验数据的学生二是想验证“时间编码到底比速率编码多表达什么”的算法工程师。2. LIF神经元与网络构建createNeuron与createNetwork的matlab实现2.1 神经元模型在createNeuron.m里的落点拿到资源后先看createNeuron.m。SNN里最常用的简化模型是LIFLeaky Integrate-and-Fire它在准确性和matlab仿真代价之间取了一个工程上能接受的点。整个过程可以写成一个很紧凑的差分方程[ C_m \frac{dV}{dt} -g_L(V - V_{rest}) I_{syn}(t) ]在离散时间仿真里输入电流用这段时间内所有到达的突触脉冲求和得到。只要膜电位超过阈值V_th就把脉冲写进spikeTrain矩阵然后把电位重置到V_reset。看createNeuron.m时重点看它以结构体数组存了哪些字段常见设计如下表字段名物理含义本仿真里影响V_th发放阈值决定神经元敏感性一般取-55mV级别V_reset重置电位不应期结束后的起点常比静息电位更低tau_m膜时间常数控制电位衰减快慢影响脉冲时序精度t_refr不应期时长ms为单位防止同一神经元连续密集发放V_m当前膜电位主循环里每一ms都会被更新2.2 createNetwork.m如何组织突触与延迟SNN和传统神经网络在结构上的最大区别是突触传导延迟。同一个神经元对权重决定“影响多大”延迟决定“晚到多久”。而polychronization这个“聚群”现象正是特定延迟组合下的产物。在createNetwork.m里常见做法是构建两个N×N矩阵W存权重D存延迟。生成连接的代码逻辑通常是% 常见做法按连接概率创建稀疏突触网络 N size(neuron, 2); % 神经元数量 connProb 0.1; % 连接概率越大网络越密 w0 0.08; % 初始权重 for pre 1:N for post 1:N if rand connProb W(post, pre) w0; % W(post,pre)pre - post D(post, pre) randi([1 20]); % 轴突延迟取1~20ms end end end这段代码的逻辑说明W的列索引是突触前神经元行索引是突触后神经元这种排布方式方便后续直接用W * spikeTrain完成所有输入电流的矩阵化求和。延迟矩阵D不参与权重计算但参与脉冲的时间对齐比如pre神经元在t时刻发放post神经元要等到tD(post,pre)时刻才接收到这个输入。参数connProb取0.1意味着平均每个神经元只有10%的神经元相连这个稀疏度在100~500个神经元的matlab仿真里很常见。2.3 主循环的时序骨架与Runme1/Runme2的分工打开Runme1.m会看到仿真时间T以ms为单位计算摘要里那一段epoch T/ms配合4重循环实际上是在按“每4ms为一个输入周期”的方式构造输入脉冲序列。这种分组的目的很明确让不同输入模式以固定间隔重复出现让STDP有机会在重复时序上强化特定突触路径。主循环的简化写法是dt 1; % 时间步长ms T 18000; % 仿真总时长来自资源中的18000.mat命名 Vm -65 .* ones(N, 1); % 初始膜电位 for t 1:dt:T I W * spk(:, t); % 所有突触输入的加权和 Vm Vm (dt/tau_m) .* (-Vm V_rest I); fired find(Vm V_th); spk(fired, t 1) 1; % 光栅图需要下一列记录当前发放 Vm(fired) V_reset; % 发放后立即重置 % 不应期处理记录每个神经元的最后发放时间t_refr内不参与积分 end说明一下这里的关键点spk(fired, t 1)是因为matlab矩阵索引从1开始如果用第t列表示第t ms的脉冲循环里已经用第t列作为输入就不能在这一列上写发放结果否则相当于让同一时刻的脉冲参与了自己的输入计算。把发放写到下一列只是一种惯用避坑方式。Runme2.m通常用来做多试次统计它会把Runme1.m训练好的W和D加载进来再跑新的输入序列保存成trial结构用于后续聚群检测。3. STDP时序学习规则权重更新与发散控制3.1 为什么第三代神经网络用STDP而不是反向传播如果试图把误差反向传播直接套到SNN上第一个障碍就是脉冲发放函数的梯度几乎处处为0脉冲的不连续性让链式法则失效。STDPSpike-Timing-Dependent Plasticity绕开了这个问题它只看一对神经元的实际发放时间关系不依赖全局误差信号。规则只有两条突触前神经元先发放、突触后神经元后发放这个突触权重就增加反过来突触前神经元后发放权重就减少。这个机制在生物实验里有充分证据在工程里则是一个完全局部的无监督学习规则。3.2 STDP.m里权重更新的常见编写方式在matlab里实现STDP一种直观的方式是让每个神经元记住自己最近一次发放时间然后在检测到某个神经元发放时把它和所有其它神经元的最近发放时间差拿出来更新权重。核心代码长这样% 最简STDP在post神经元j发放时更新它所有输入权重 lastSpike -inf(N, 1); % 记录每个神经元最近发放时间 A_plus 0.1; % 增强幅度 A_minus 0.12; % 削弱幅度略大于增强幅度 tau 20; % STDP时间窗口常数ms for t 1:T firedPost find(spk(:, t) 1); for j 1:length(firedPost) jn firedPost(j); dtPre t - lastSpike; % 其它神经元的“最近发放离现在多久” potent dtPre 0 dtPre 5 * tau; W(jn, :) W(jn, :) A_plus * exp(-dtPre / tau) .* potent; % 对突触前神经元i如果它比post早发放dtPre0权重增加 lastSpike(jn) t; end end这里需要解释几个参数的工程含义。A_minus略大于A_plus是经验做法目的在于抑制过度同步放电避免网络进入所有神经元一起发放的死锁状态。tau取20ms意味着只有时间差在100ms以内的脉冲对才会对权重产生明显影响。代码里的potent是一个0/1掩码把太旧的时间差过滤掉一方面减少计算量另一方面防止无关脉冲对之间产生虚假的权重变化。完整的STDP.m里通常还会加上权重上下界判断用一行W(W w_max) w_max做硬截断防止某个突触被反复增强到N倍。3.3 仿真发散网络行为失控时看哪些指标SNN仿真最常见的失败现象是“发散”表现特征很有辨识度脉冲光栅图上突然出现一条竖直的密集点阵几乎全部神经元在同一时刻发放随后权重矩阵出现Inf值或者全变成w_max。这不是matlab的bug而是学习率与抑制机制不匹配。处理方式有两个方向。第一把A_plus下调一个数量级比如从0.1改到0.01让每次局部更新对网络的影响更小挺过最开始的敏感期。第二检查神经元平均发放率随时间的变化% 检查每个时间窗口的平均发放率判断是否发散 rateWin sum(spk, 1) / N; % 每ms所有神经元的平均发放数 plot(rateWin); % 如果曲线呈上升趋势且持续超过0.1说明抑制不足如果发散窗口发生在某一时刻还可以用find(rateWin 0.5)把那次同步暴发的时间点找出来往前追溯最近一次输入模式切换通常会发现是某个delays组合在特定输入下产生了超过阈值的共振。这恰恰说明polychronization现象本身的形成机制和暴发只有一线之隔。4. Polychronization聚群检测polygroup与polychron的仿真输出4.1 时间编码里的聚群到底是什么前面提到突触延迟矩阵D让同一输入可以在网络里产生“早到晚到”的空间分布。2006年Izhikevich在Cerebral Cortex上发表的论文里给出了一个漂亮的结果在随机脉冲输入刺激下固定延迟的网络加上STDP会自发形成一群神经元这些神经元总是在特定的时间次序上依次发放就像一组预先约定好暗号的成员。他把这种现象命名为polychronization对应的matlab实现函数就是资源里的polygroup.m和polychron.m。两者的分工不同。polygroup.m侧重在空间上把神经元分组——哪些神经元经常形成共同激活的集合polychron.m侧重在时间上检测重复出现的发放次序——集合内神经元的相对时间差在多次试次中是否保持一致。4.2 四次函数在聚群检测中的简化实现把资源跑通后最值得改的代码是聚群检测部分。一个不太复杂但能验证概念的matlab实现思路是对每一对神经元统计它们在不同试次中以固定时间差由延迟矩阵D决定共同发放的次数。如果这个统计量超过阈值就认为这对神经元属于同一个聚群。% 简化版聚群检测统计固定时间差下共同发放的次数 nTrials 20; % 试次数 matchCount zeros(N, N); window 2; % 时间窗 ±2ms容忍抖动 for trial 1:nTrials st spikeData{trial}; % N x T 的稀疏脉冲矩阵 for i 1:N tsi find(st(i, :)); % i神经元所有发放时刻 for j 1:N dDelay D(j, i); % 期望的时间差 for tm tsi seg tm dDelay - window : tm dDelay window; seg seg(seg 1 seg size(st, 2)); if any(st(j, seg)) matchCount(i, j) matchCount(i, j) 1; end end end end end % 重叠发放比例超过80%的神经元对视为同一聚群 [gi, gj] find(matchCount / nTrials 0.8);这段代码的三个嵌套循环在matlab里跑起来不便宜但胜在直观。参数window2表示允许±2ms的抖动这是因为数值积分和STDP更新会带来微小的相位漂移实际检测时不必要求严格对齐到毫秒。注意seg要做边界裁剪否则st(j, seg)很可能因为索引越界而报错。如果希望提速可以先用find(sparse(st), first)把稀疏矩阵转成元胞数组或者用conv2做卷积匹配但验证算法逻辑时这个版本更容易对照论文理解。4.3 plotMe.m如何验证聚群存在脉冲光栅图的读法是判断STDP训练有没有真正起作用的关键。运行plotMe.m后应该看到的是这样一个散点图横轴是时间ms纵轴是神经元编号每个点代表一次脉冲发放。在没有聚群时脉冲点随机分布成均匀的云团出现聚群时云团里会出现明显的重复斜线纹理——一条斜线意味着一组神经元按照相对固定的时间间隔连续发放。斜线出现的频率越高、重复度越好说明STDP强化的突触链越清楚。跑完Runme1.m后如果光栅图一片空白先看输入向量input(1,:)和input(2,:)里是不是确实有1存在再看createNetwork.m里的connProb是不是太小。这两处是最容易让人误以为“网络没动力学”的实际原因。4.4 录像里Runme1到Runme2的操作顺序操作录像里演示的顺序是固定的先运行Runme1.m完成训练等界面输出权重热图后再运行Runme2.m做多试次测试最后运行polygroup.m和plotMe.m查看聚群结果。要注意Runme2.m依赖Runme1.m训练出的18000.mat如果只运行Runme2.m会提示找不到文件。18000这个数字不是固定仿真时长而是训练样本切换周期对应摘要中4ms一个pattern的重复次数。第一次调试时可以把训练时间缩短到四分之一先用最快的路径确认代码能完整走通再跑全量数据避免每改一个参数都要等十几分钟。5. 按录像复现时的路径、报错与可视化验证5.1 左侧当前文件夹路径是第一个坑打开matlab后第一件事不是点Run而是看左侧“当前文件夹”是否已经切换到code目录。因为所有脚本之间都是通过相对路径调用比如createNetwork.m在同目录下才能被Runme1.m直接用函数名调用。如果当前路径不对最常见的报错是Undefined function createNetwork这时候代码一行都没执行和算法没任何关系。建议先执行一次cd(D:\SNN_code\code); % 替换成资源解压后的实际路径 pwd % 确认当前路径 which Runme1.m % 能返回文件路径说明目录切换成功5.2 几类运行期报错的排查顺序第一类报错是矩阵维度不匹配尤其出现在摘要那几行loc fix(4*(i-1)*ms find(input(1,:) 1) * ms)上。fix取整后可能得到0或者超过spikeTraing矩阵列数的值这两个情况都会让索引报错。调试时在loc赋值后加一行loc(loc 1 | loc size(inputSpikes, 2)) [];直接过滤越界索引即可。第二类报错是权重出现NaN多半是exp(-dtPre/tau)里dtPre为负数且绝对值很大或者tau被误改成0导致除零。检查时可以运行dbstop if naninfmatlab会在NaN出现的第一个位置停下来快速定位是哪次循环产生了非法值。5.3 一个快速验证网络确实在学习的技巧不依赖plotMe.m手工检查一条突触路径的权重变化load 18000.mat; % 训练完的结果 wPath W(10, 5); % 神经元5 - 神经元10 的突触权重 fprintf(突触权重为 %f\n, wPath); % 对比 Runme2.m 跑完后再打印一次两个值应当不同如果两次运行这个值完全相同说明STDP没有触发更新请检查STDP.m里的lastSpike初始化是否用了-inf以及spk矩阵是否记录了正确的发放时刻。本文还有配套的精品资源点击获取