第三代神经网络——脉冲神经网络（SNN）的仿真

1.问题描述：

          同传统的人工神经网络一样，脉冲神经网络同样分为三种拓扑结构。它们分别是前馈型脉冲神经网络（feed-forward spiking neural network）、递归型脉冲神经网络（recurrent spiking neural network）和混合型脉冲神经网络（hybird spiking neural network）。

学习是人工智能领域的核心问题，对于SNN来说，基于脉冲时间层次的学习方法研究，对于通过理论模型来验证生物神经系统的信息处理和学习机制是必须的。通过生物可解释的方式建立人工神经系统，科学家希望可以通过神经科学和行为实验来达到预期目的。大脑中的学习可以理解为突触连接强度随时间的变化过程，这种能力称为突触可塑性（synaptic plasticity）。脉冲神经网络的学习方式主要包括无监督学习（unsupervised learning）、监督学习（supervised learning）和强化学习（reinforcement learning）等。

1. 无监督学习算法

无监督学习算法在人类和动物的学习中占据主导地位，人们通过观察能够发现世界的内在结构，而不是被告知每一个客观事物的名称。人工神经网络无监督学习算法的设计主要是针对无标签数据集的训练，要求应用无监督学习规则对神经网络中的连接权值或结构进行自适应的调整。也就是说，在没有“教师”信号的监督下，神经网络必须自己从输入数据中发现规律性（如统计特征、相关性或类别等），并通过输出实现分类或决策。一般来说，只有当输入数据集中存在冗余性时，无监督学习才有意义，否则，无监督学习不能很好地发现输入数据中的任何模式或特征，即冗余性提供了知识。

脉冲神经网络的无监督学习算法大多是借鉴传统人工神经网络的无监督学习算法，是在Hebb学习规则不同变体的基础上提出的。神经科学的研究成果表明，生物神经系统中的脉冲序列不仅可引起神经突触的持续变化，并且满足脉冲时间依赖可塑性（spike timing-dependent plasticity，STDP）机制。在决定性时间窗口内，根据突触前神经元和突触后神经元发放的脉冲序列的相对时序关系，应用STDP学习规则可以对突触权值进行无监督方式的调整。

2. 脉冲神经网络的监督学习

脉冲神经网络的监督学习是指对于给定的多个输入脉冲序列和多个目标脉冲序列，寻找脉冲神经网络合适的突触权值矩阵，使神经元的输出脉冲序列与对应的目标脉冲序列尽可能接近，即两者的误差评价函数最小。对于脉冲神经网络来说，神经信息以脉冲序列的形式表示，神经元内部状态变量及误差函数不再满足连续可微的性质，构建有效的脉冲神经网络监督学习算法非常困难，同时也是该领域的一个重要的研究方向。

根据监督学习所采用的基本思想不同，可以将现有的监督学习算法分为三类：

1. 基于梯度下降的监督学习算法的基本思想是利用神经元目标输出与实际输出之间的误差以及误差反向传播过程，得到梯度下降计算结果作为突触权值调整的参考量，最终减小这种误差。基于梯度下降的监督学习算法是一种数学分析方法，在学习规则的推导过程中，要求神经元模型的状态变量必须是有解析的表达式，主要采用固定阈值的线性神经元模型，如脉冲响应模型（spike response model）和Integrate-and-Fire神经元模型等。
2. 基于突触可塑性的监督学习算法的基本思想是利用神经元发放脉冲序列的时间相关性所引起的突触可塑性机制，设计神经元突触权值调整的学习规则，这是一种具有生物可解释性的监督学习。
3. 基于脉冲序列卷积的监督学习算法通过脉冲序列内积的差异构造脉冲神经网络的监督学习算法，突触权值的调整依赖于特定核函数的卷积计算，可实现脉冲序列时空模式的学习。

3. 脉冲神经网络的强化学习

强化学习是从环境状态到行为映射的学习，以使智能体行为从环境中获得的累积奖赏值最大。基于生物启发的学习机制，人工神经网络强化学习的研究重点在于探索智能体的自适应优化策略，是近年来神经网络和智能控制领域的主要方法之一。强化学习关注的是智能体如何在环境中采取一系列行为，通过强化学习，一个智能体应该知道在什么状态下应该采取什么行为。可以看到，强化学习和监督学习的区别主要在于以下两点：

1. 强化学习是试错学习，由于没有直接的“教师”指导信息，智能体要不断与环境交互，通过试错的方式来获得最佳策略；
2. 延迟回报，强化学习的指导信息很少，而且往往是在事后（最后一个状态）才给出的，这就导致了一个问题，即获得正回报或者负回报以后，如何将汇报分配给前面的状态。

四、脉冲神经网络的进化方法

进化算法（evolutionary algorithm）是模拟生物进化过程的计算模型，是一类基于自然选择和遗传变异等生物进化机制的全局性概率搜索算法，主要包括遗传算法（genetic algorithm）、进化规划（evolutionary programming）和进化策略（evolutionary strategy）等。虽然这些算法在实现方面具有一些差别，但它们具有一个共同的特点，即都是借助生物进化的思想和原理来解决实际问题的。

将进化算法与脉冲神经网络有机结合起来，研究者开辟了进化脉冲神经网络（evolutionary spiking neural network）的研究领域，以提高对复杂问题的求解能力。进化脉冲神经网络可以作为自适应系统的一种一般性框架，在没有人为干预的情况下系统自适应调整神经元的参数、连接权值、网络结构和学习规则。

2.部分程序：

rand('seed',1);
% spnet.m: Spiking network with axonal conduction delays and STDP
% Created by Eugene M.Izhikevich.                February 3, 2004
% Modified to allow arbitrary delay distributions.  April 16,2008
M=100;                 % number of synapses per neuron
D=20;                  % maximal conduction delay 
% excitatory neurons   % inhibitory neurons      % total number 
Ne=800;                Ni=200;                   N=Ne+Ni;
a=[0.02*ones(Ne,1);    0.1*ones(Ni,1)];
d=[   8*ones(Ne,1);    2*ones(Ni,1)];
sm=10;                 % maximal synaptic strength

% post=ceil([N*rand(Ne,M);Ne*rand(Ni,M)]); 
% Take special care not to have multiple connections between neurons
delays = cell(N,D);
for i=1:Ne
    p=randperm(N);
    post(i,:)=p(1:M);
    for j=1:M
        delays{i, ceil(D*rand)}(end+1) = j;  % Assign random exc delays
    end;
end;
for i=Ne+1:N
    p=randperm(Ne);
    post(i,:)=p(1:M);
    delays{i,1}=1:M;                    % all inh delays are 1 ms.
end;

s=[6*ones(Ne,M);-5*ones(Ni,M)];         % synaptic weights
sd=zeros(N,M);                          % their derivatives

% Make links at postsynaptic targets to the presynaptic weights
pre = cell(N,1);
aux = cell(N,1);
for i=1:Ne
    for j=1:D
        for k=1:length(delays{i,j})
            pre{post(i, delays{i, j}(k))}(end+1) = N*(delays{i, j}(k)-1)+i;
            aux{post(i, delays{i, j}(k))}(end+1) = N*(D-1-j)+i; % takes into account delay
        end;
    end;
end;
  

STDP = zeros(N,1001+D);
v = -65*ones(N,1);                      % initial values
u = 0.2.*v;                             % initial values
firings=[-D 0];                         % spike timings

for sec=1:60*60*24                      % simulation of 1 day
  for t=1:1000                          % simulation of 1 sec
    I=zeros(N,1);        
    I(ceil(N*rand))=20;                 % random thalamic input 
    fired = find(v>=30);                % indices of fired neurons
    v(fired)=-65;  
    u(fired)=u(fired)+d(fired);
    STDP(fired,t+D)=0.1;
    for k=1:length(fired)
      sd(pre{fired(k)})=sd(pre{fired(k)})+STDP(N*t+aux{fired(k)});
    end;
    firings=[firings;t*ones(length(fired),1),fired];
    k=size(firings,1);
    while firings(k,1)>t-D
      del=delays{firings(k,2),t-firings(k,1)+1};
      ind = post(firings(k,2),del);
      I(ind)=I(ind)+s(firings(k,2), del)';
      sd(firings(k,2),del)=sd(firings(k,2),del)-1.2*STDP(ind,t+D)';
      k=k-1;
    end;
    v=v+0.5*((0.04*v+5).*v+140-u+I);    % for numerical 
    v=v+0.5*((0.04*v+5).*v+140-u+I);    % stability time 
    u=u+a.*(0.2*v-u);                   % step is 0.5 ms
    STDP(:,t+D+1)=0.95*STDP(:,t+D);     % tau = 20 ms
  end;
  plot(firings(:,1),firings(:,2),'.');
  axis([0 1000 0 N]); drawnow;
  STDP(:,1:D+1)=STDP(:,1001:1001+D);
  ind = find(firings(:,1) > 1001-D);
  firings=[-D 0;firings(ind,1)-1000,firings(ind,2)];
  s(1:Ne,:)=max(0,min(sm,0.01+s(1:Ne,:)+sd(1:Ne,:)));
  sd=0.9*sd;
end;
 

3.仿真结论：

D96