matlab - 一种更好的在MATLAB中实现alpha突触功能的方法

Question

我已经实现了一个泄漏集成和激发神经元系统，它以 alpha 函数 alpha = t/tau * exp(1 - (t/tau)) 的形式提供输出。但是，我用来实现它的代码片段占用了至少 80% 的运行时间(总共 5 秒中大约有 4 秒)。在程序运行过程中，此 prt 至少被调用 30000 次，而 computealphasynapseoutput 函数至少被调用 150 万次。所以我想减少这部分的运行时间。我曾尝试使用 arrayfun 来实现它，但这比这要花费更多的时间。

有人可以建议更有效地实现此代码吗？

为了实现 alpha 突触，我使用了以下代码:

% Get the identity of the currently active neurons
idAllActiveNeuron = idAllActiveNeuron > 0;
if any(idAllActiveNeuron) % Only run if atleast one neuron is active
    for iActiveNeuron = find(idAllActiveNeuron).' % run for each active neuron
        %% synapticOutputArray stores the synaptic output for each neuron for each time instant
        %% iIntegration is the time instant
        %% spikeTimesArray is a cell array that is composed of spike times for 
        % each neuron. So if I have 5 neurons and neuron 4 spikes twice 
        % while neuron 5 spikes once, spikeTimesArray would be something 
        % like {[], [], [], [0.0023, 0.0034], [0.0675]}
        %% integrationInstant would be a time value like 0.0810
        %% tau_syn stores the value of tau for each neuron

        synapticOutputArray(iActiveNeuron, iIntegration) = computealphasynapseoutput(spikeTimesArray{iActiveNeuron}, integrationInstant, tau_syn(iActiveNeuron));

    end % iActiveNeuron
end

函数computealphasynapse实现如下:

function synapticOutput = computealphasynapseoutput(firingTime, integrationInstant, tauSyn)
%%COMPUTEALPHASYNAPSEOUTPUT Calculates the synaptic output over all
%previous spikes of the neuron at a particular time instant using the
%alpha synapse function
%
% Usage:
%   synapticOutput = computealphasynapseoutput(firingTime, integrationInstant, tauSyn)
%
% Inputs:
%           firingTime: vector of previous firing times (in seconds)
%   integrationInstant: current integration time instant (in seconds)
%               tauSyn: synaptic time constant of the neuron (in seconds)
%
% Output:
%   synapticOutput: Synaptic output of the neuron at current time summed
%                   over all firing instances
%

% Calculate the time difference of firing from current time
timeDifference = (integrationInstant - firingTime) / tauSyn;
% Calculate and sum the synaptic output for each spike
synapticOutput = sum(timeDifference .* exp(1 - timeDifference));

end % computealphasynapseoutput

编辑:

I have finally awarded the bounty for this question to gnovice for his awesome answer. It helped me to shave off an entire 40 seconds from my simulation time (from 68 to 28 s). I hope it works for people in the future too. I would also like to acknowledge MrAzzaman and qbzenker for taking their time to answer the question and teaching me some cool new approaches. Also, others who commented and helped me. Thanks

Answer 1

这个解决方案可能有点奇怪，而且非常针对手头的问题，但非常有效。这个想法是重新组织用于计算总和 alpha 突触输出的公式，以最小化完成的计算，特别是减少使用 exp 计算的求幂次数。 ，这在计算上是昂贵的。这是单个神经元的重新表述:

对于结果，我们有一个时间函数 f(t) (计算当前时间点的一个指数 t )乘以 t 中的线性函数具有不随时间变化的参数 A和 B .注意 A和 B仅取决于尖峰时间，对所有先前尖峰出现的指数求和。我们可以存储这两个参数 A和 B对于每个神经元，当神经元中出现新的尖峰时，我们只需计算一对新的项并将它们添加到 A 的现有值中。和 B .

注意事项:注意指数的值。对于 A和 B值为 ts/tau , 如果它变得太大，它可能会溢出计算，导致 Inf .对于我们在这里使用的参数，这不会发生。 ts 的值最大值为 1(因为我们只模拟一秒)，最小值为 tau。是0.01秒。指数最多为 100，给出的值很大(大约 10^43)，但很容易被 double 处理。多变的。同样，f(t) 中的指数将有一个最大的负值 -99，给出非常小的值(大约 10^-43)，但仍然很容易被 double 处理。变量不会下溢到 0。但是，请记住，使用较小的 tau值或更长的仿真时间可能会导致问题。

那么，我们如何实现呢？以下是您必须添加/修改的相关代码片段(请注意，您甚至不再需要 spikeTimesArray):

% You'll have to initialize storage for A and B:
A = zeros(1, nNeurons);
B = zeros(1, nNeurons);
idAllActiveNeuron = false(1, nNeurons);

...

% When you are generating new spikes, modify A and B like so:
index = ...;  % A vector of neuron indices where new spikes occur this integration step
if ~isempty(index)
  expTerm = exp(integrationInstant./tau_syn(index));
  A(index) = A(index)+expTerm;
  B(index) = B(index)-integrationInstant.*expTerm;
  idAllActiveNeuron(index) = true;
end

...

% Updating the synaptic output no longer requires a loop:
if any(idAllActiveNeuron)
  synapticOutputArray(idAllActiveNeuron, iIntegration) = ...
    (integrationInstant.*A(idAllActiveNeuron)+B(idAllActiveNeuron)).*...
    exp(1-integrationInstant./tau_syn(idAllActiveNeuron))./tau_syn(idAllActiveNeuron);
end

它的表现如何？以下是我对 1101 个神经元进行的一些测量(使用 timeit )，使用不同的平均尖峰发射频率模拟，时间步长为 0.1 毫秒一秒(10001 个时间点):

出现的峰值越多，完成原始解决方案所需的时间就越多。然而，重新配制溶液所需的时间几乎恒定，从未超过 0.6 秒。我还可以确认这两种方法的输出是相等的(突触输出波形的最大差异约为 10^-13)。