定点神经ODE求解器的精度和效率

在计算神经科学中，大多数神经元模型都是以普通形式指定的微分方程（ODE）或偶尔的偏微分方程（PDE）通常将其本身简化为ODE以获得有效解决方案。中的状态变量该方程模拟了神经元在任何特定时间点的内部状态ODE描述其动态行为。状态变量在模拟，然后通过求解ODE来计算神经元行为的演化每个时间点，考虑状态变量和输入历史。有一个关于初值常微分方程数值解的大量文献，至少可以追溯到1770年欧拉的作品（巴什福斯和亚当斯，1883; 伦格，1895; 希恩，1900; 库塔，1901; 欧拉，1913; 齿轮，1971; 兰伯特，1973，1991; 霍尔和瓦特，1976;屠夫，2003). 有多种方法可用，并且在三个关键性能度量准确性之间存在许多权衡，稳定性和效率（在我们的上下文中，最后一个定义为“每个神经所花费的时间更新”）。

这封信旨在评估各种选择，并找到有偏见的实际解决方案一组特定的情况，同时，我们希望，也强调一些更一般的情况原则和可能性。为了在在空间约束下进行分析，我们选择只研究Izhikevich模型参数化、输入行为和初始条件的有限变化集。我们希望其中的一些想法和方法能更广泛地应用。一系列情况我们关注的是

1. 32位定点算术

2. 主要是运行时效率动机

3. Izhikevich神经模型

状态变量中这种不连续性的存在将对最合适的求解方法和所需的细化在离散时间点计算仿真时的近似连续时间。它应该注意的是，目前许多神经模型都具有类似的不连续性（通常称为一体式模型）。伊日凯维奇(2003)通过选择{一，b条，c（c），d日}适当地，一个人可以复制许多真实的神经行为的类型。其他研究着眼于数值解的选项这类方程是汉弗莱斯和格尼(2007)和图布尔(2010).

在下面的部分中，我们将展示在考虑的上下文中，一些ODE解决方案就效率和准确度和技术可以通过适当考虑说明SpiNNaker仿真系统和算法硬件的特定性质。

考虑可用于求解初值的各种方法的简便方法ODE是根据算法的性质对其进行分类。第一个考虑因素是显式的还是隐式的。隐式解算器需要一组固有的每个时间点的非线性方程，通常通过迭代算法-尽管是直接的在某些情况下，可以使用线性近似求解。它们提供了理想的遇到刚性数值行为时改进稳定性的特性这可能发生在求解神经ODE时。常微分方程背景下的刚度在屠夫(2003)：“这些系统的特点是非常高的稳定性，可以变成非常高的不稳定性什么时候用标准数值方法近似”（第27页）。

这种额外的数值稳定性并不是免费的，因为它比计算成本高，因此考虑到我们的目标计算能力有限的实时计算。相反，显式解算器使用已经可用或易于计算的值，然后预测状态值变量在下一个时间点使用直接代数公式。显式解算器各不相同其稳定性特征广泛。

另一个考虑因素是固定步长与自适应步长。固定步长使用常量向前移动模拟的时间步长，而自适应步长决定步长基于每个时间点解决方案的性质。自适应步长为所有常见的数字精度要求提供了理想的特性，效率和稳定性，如果可行，应始终予以考虑。然而，由于关键SpiNNaker体系结构上的大规模神经系统仿真方面从基本的固定时间步长行为移开的复杂性和开销系统作为一个整体使得自适应时间步长的全局使用不可行。本着因此，我们在本研究中使用了固定的时间步长行为（除来自用于比较的参考结果）。重要的是，当查看结果，因为它不可避免地会影响解决。

方法的顺序表示预测误差随着步长的增加而提高的速度减少，但顺序的确切含义因方法而异考虑过的。例如，较高的阶数通常会减少预测误差，但对于一些方法会增加稳定性，而其他方法则会降低稳定性。它需要承担请记住，在大多数情况下，这些结果背后的复杂数学分析假设精确的算术。这会使实际行为与正式行为大不相同预期，如以下各节所强调的。

屠夫(2003)提供了有用的最常用的数值算法分为以下三个维度：多步骤使用状态变量的先前值；多级计算一阶导数当前时间点和下一时间点之间的点；和多衍生用途的高级衍生品在当前时间点计算。与我们的要求相关的显式方法示例以下是（有关详细信息，请参阅Gear，1971;兰伯特，1973):

• •

  多步骤：线性多步骤方法-Adams-Bashforth in最简单的形式

• •

  多阶段：伦格-库塔

• •

  多重导数：泰勒级数和Parker-Sochacki

基本欧拉方法是所有三种类型中最简单的示例，因此，在某种意义上，定义此三维空间中的原点。此外，混合方法结合了一些或所有类型都是可能的。对这些潜在有趣内容的详细讨论混合动力车超出了本信的范围，但稍后会提到一些可能性。安外推法不适合这个方案；这个最著名的例子是Gragg-Bulirsch-Stoer（参见Lambert，1973; Press、Teukolsky、Vetterling和Flannery，1992; Stewart&Bair，2009). 这可能是非常准确的平滑解决方案时的首选行为系统是必需的，但当状态不连续时不建议使用变量是预期的，或者当导数计算非常重要时。因此，我们将初始考虑限制为以下明确的固定时间步长求解方法：线性多步法、Runge-Kutta、Taylor级数和Parker-Sochacki。

所有这些都是提供不同顺序的方法族。对于的Runge-Kutta方法阶数为2或更高时，每个阶数都有无限种算法，可通过定义要计算的显式方程的自由参数。这些考虑因素以屠夫画面的形式表达，对其进行了最全面的描述由Butcher提供(2003). 从计算上看Runge-Kutta方法需要对每个状态变量进行一次导数计算。这个泰勒级数方法在这方面类似于龙格库塔法，但在分析形式上是独特的对于每个订单，Parker-Sochacki也是如此。相关的显式线性多步方法（Adams-Bashforth）只需要对每个状态变量进行一次导数计算，与订单，使其潜在的效率非常高。

为了保持数值研究的可管理性，我们进一步修正了以下方法考虑过的。Parker-Sochacki方法使用截断的Picard迭代计算Maclaurin系列已经成功应用于Izhikevich模型，其中其他（关于这些方法如何结合在一起的细节，参见Stewart&Bair，2009). 它提供了许多额外的机会例如，通过自动调整以提供所需级别的幂级数、有理函数和超越函数的精确自然表达迭代多项式运算。这当然值得进一步研究，但首先对所需计算结果的分析表明，它不太可能具有竞争力就效率而言，这里和那里都没有被进一步考虑。

Adams-Bashforth方法在效率方面有很大潜力，但稳定性是差，尤其是对于更高的订单，并且在行为僵化的情况下，不太可能以实际的时间步长可靠地收敛。将其纳入预测-校正有时使用亚当斯-莫尔顿校正器和局部外推技术的算法称为Milne的方法显著增加了稳定性，但也增加了计算负荷（参见兰伯特，1991，了解此组合的详细信息，以及他的图4.1和4.2用于前后比较）。此外，它还需要额外的保留过去值的内存-需要的本地快速内存在SpiNNaker芯片并不能很好地处理状态变量中的不连续性没有额外的管理代码和相关开销。

因此，评估并相互比较的最终方法以及与使用Mathematica中提供的复杂ODE解算器计算的参考解环境（Wolfram Research，2014)是龙格库塔和泰勒级数.

与标准数值计算环境相比，其中浮点类型和操作可用，使用定点算法和各种约束随之而来的可能是一个挑战。历史上，定点算法通常使用整数算术的手工编码转换执行（例如在密切相关的上下文中，请参阅Jin、Furber和Woods，2008). 这种方法的优点是可以为生成自定义类型可以针对速度优化特定问题和关键算术运算，包括在适当的地方使用汇编语言。然而，这使得代码更加困难编写、阅读和维护，容易出现难以发现的错误。因此，它不是被认为是需要共享思想和代码的大型软件工程的理想选择群体之间，尤其是当这些群体可能包含更多神经科学家时，数学家和工程师，而不是计算机科学家。幸运的是，最近的草案标准TR18037 ISO/IEC(2008)提供定义明确的通过使用定点的类型和算术操作解决这个难题算法，同时允许C代码看起来与通常使用的代码非常相似浮点类型和操作，并执行所有必需的操作由编译器自动执行。本标准草案已在GCC中实施自4.7发布以来，开发工具链仅针对ARM目标。作为SpiNNaker架构基于ARM，我们可以透明地将此方法应用于我们的数值计算，因此，除非另有说明，否则本信其余部分均在此假设下工作否则。我们选择的默认类型在ISO草案中定义为累计-一个32位有符号类型，16个整数位和15个小数位（s16.15），允许ARM处理器上的高效算术操作。

考虑到这一点，其他重要问题是数字，尤其是，溢出、下溢和量化/舍入错误。如前所述第节中，ODE解算器的精度和稳定性的大多数标准结果都基于精确算术。即使是双精度浮点，我们离这个理想还有很长的路要走在我们使用s16.15类型的目标32位体系结构中，其他问题也会出现投入比赛。已开展工作，专门研究舍入误差对初始值ODE（单位：吉尔，1951; 亨利西，1962; 巴布什卡、普雷格和维塔塞克，1966; 卡汉，1965; 维塔塞克，1969; 海姆，1996)，但最近的大部分工作都是在浮点运算的背景下进行的，它具有不同于固定点的特性，特别是各种补偿建议用于浮点算法的求和方案（例如，参见海姆，1996)在定点中没有任何好处计算。

最早的三个参考文献是例外，吉尔(1951)为应用传统的四阶Runge-Kutta方法在不动点算法中的应用计算中的舍入错误。亨利西(1962)提供了四舍五入误差统计分布的分析和示例显式算法的数量。Babuška等人的第2章和第3章。(1966)提供定点和浮点的分析与ODE解决方案有关的错误及示例。尽管细节千差万别在他们的示例和方法之间，出现了一个特征V形图案随着步长的减小而减小，直到舍入/量化导致误差的点当步长进一步减小时再次上升。这种图案的斜率取决于细节以及使用的是浮点算法还是定点算法。

虽然有趣且与我们的一些结果相关这两种算术之间的区别超出了这封信的范围。一个然而，显而易见的是，舍入误差的敏感性是很难进行分析表征，因为可以实现相同的算法例如，在中期变量的选择和排序方面有所不同其本身可以对最终结果的准确性产生重大影响。

我们在每种情况下对最佳算法的搜索都是完全手动的，显然无法证明最优的。然而，在这种情况下，它是相对彻底的，着眼于大量数据变量和运算的可能方程分解和排序。即使有这种有限的方法，在如果按照代数定义考虑，则相同。我们希望至少展示了该方法的潜力。已通过以下方式完成此过程的自动化工作马特尔(2009)，Darulova等人。(2013)、高、拜利斯和君士坦丁尼德斯(2013)，这似乎是一个潜在的富有成效的方法继续。我们需要注意的是，定点类型需要非常多的坚定的先验知识，以便将过程转化为自动化方法，尤其是，严格保证算法中的变量不会溢出或下溢在任何可能的用例中都有代表性的范围，而且这些有时更难提供比最初看起来更合适的服务。

使用的关键数字细节累计s16.15 ARM上的类型目标是：

1. 绝对值0.000031变为零。

2. 绝对值包大约（饱和类型可用，但很多速度较慢）。

3. 精度是绝对的，而不是相对的带浮点。

4. 划分非常缓慢（部分原因是ARM指令集的详细信息），并且应该避免。

5. 不存在本机超越函数图书馆（尽管我们小组正在开发一个）。

正如所料，这些点对我们的算法有直接影响在这里进行讨论，在考虑实现和结果。

SpiNNaker的系统语言是ANSI C，因此这是用于执行。开发了一个API，它提供了一个通用但高效的接口可能感兴趣的神经模型API和ODE解算器之间。一个插件选择了能够测试新算法的体系结构，同时提高编译器优化、内存访问等问题的效率模式和寄存器的使用。经过一些考虑，我们找到了一个结构类似于Press等人第16章中使用的方法。(1992)而是关注我们的特殊需求。

API引用通用真实类型，可以在在我们的情况下，编译时间到任何类型的兴趣，累计，长的累计，浮动、和双重的。这允许头中只有一个类型定义的所有模型和算法的端到端比较文件。

ODE解算器API具有内部减小解算步长的能力与整个仿真时间步长无关。这是一种有用的增加效率线性成本的准确性（如需要）。解算器类型通常作为代码中的函数指针，以便易于更改，因此非常适合比较工作。然而，要为这种普遍性付出的代价是在每次调用时取消引用指针。对于生产代码，此开销可以减少，每一个连续的专门化都允许编译器减少所需的指令用于求解例程。如果事先知道有多少，这可以进一步专门化通过在解算器中对内部循环进行硬编码，状态变量将位于神经元定义中。效率提高的一些示例在第节中给出7.

另一种与专业化正交的速度改进可能性描述了静态内联函数的使用。使用此技术，代码在中定义头文件而不是源文件，类似于C++中常用的模板。尽管软件包现在在构建和交付方面不那么干净，但技术允许编译器进一步优化机会，特别是代码跨源文件的优化，这些源文件通常单独编译并链接到构建时间。该方法的效率优势也将在第节中讨论7.

在考虑实施时产生的一个重要想法如下。与任何在特定情况下，我们使用的是固定步长、显式解算器和已知ODE任何给定时间步长的计算都可以用代数来描述，而不是在算法上，采用常规的实现方法。然后可以对其进行操作，以便例如，简化并找到通用术语，这将减少计算开销并允许编译器查找更多优化。我们称之为显式求解器简化（ESR），事实证明，它不仅能够提高效率，而且选择约简可以帮助定点算法保持更高的精度并避免溢出和底流条件，因此也提高了精确度。我们在例如。我们称这些准闭式解为近似闭式解解决。其他软件工程的好处是使用静态内联函数不再需要，因为所有相关的编译器优化现在都可以在单个定义计算和任何给定方法/神经元组合的源文件交付预测试。

Stimberg、Goodman、Benichoux和布雷特(2014). 虽然我们已经独立地结论是，将ODE和求解器组合成代数有好处方程、动机和使用的方法不同。对于他们来说，重点是描述模型的通用机制、新思想的简单测试和自动代码世代；对我们来说，它被视为寻找最佳的代数方法准确度和效率之间的计算权衡，同时考虑考虑定点算法的特殊约束。

像大多数其他神经激发的ODE一样，Izhikevich模型是自主的，这意味着方程组没有明确提及右侧的时间变量；他们是神经参数与状态变量值的函数参考。这使得泰勒级数方法可以提前用代数方法实现使用Alsaker描述的优雅方法的任意数量的状态变量(2009). 这意味着这些解算器/神经元组合可以直接转换为准闭式解，并描述了所有优点之前。

以前的工作使用了不同的默认时间步长，通常为0.1 ms或1 ms。这个SpiNNaker上的默认值选择为1ms，以提高生物实时性能。作为一个结果，除非另有说明，大多数比较工作都使用这个时间步长。随着时间的推移，为了实现ODE解决方案的更精确性和具有短时间常数的突触，因此一些结果也将显示在s用于比较。

许多神经ODE的一个特殊特征是当一个神经元处于状态时出现不连续性变量击中截止点并触发。在Izhikevich模型的情况下此时变量的值会立即改变。这会产生两个数值解的要求：

1. 一个解决方案不受这些不连续性影响的方法

2. 一个阈值交叉点介于时间之间的校正机制步骤

第一个问题已经讨论过了。第二个我们在这里描述为时间量化（TQ）需要更多的思考。

任何全面的ODE解决方案代码都应处理TQ的正确处理，但纳ve代码将忽略此问题并处理截止点作为超过阈值后的下一个时间步。这并不难确保这将不可避免地导致国家进步的累积滞后变量。假设时间步之间交叉点的随机均匀分布这意味着这种方法会累积神经元每次激发的时间步长。改进是可能的，但是这需要更多的处理，因此必须高效地完成。最简单的方法是添加神经元按顺序启动后，状态变量的第一次进化需要一段时间让它弥补失去的时间。可以选择的明显值是预期值该损失的价值：时间步长，以便在触发事件接收时间值1时间步长。这消除了预期偏差，并提供了最大误差时间步长，引入到峰值时间约为0.29时间步长，平均绝对偏差为时间步长。这个简单的方案执行得很好，并且具有较低的间接费用。可以使用更复杂的方案，但最准确的方案需要一些迭代形式（例如，Stewart和Bair第3.2节中的Newton-Raphson，2009)因此，考虑到我们的具有挑战性的实时目标。

计算如图所示1，以及下面的代码片段，这里使用30 mV的截止值和1 ms的时间步长上一个点位于0，当前点位于1。对于绿线，括号显示一个=20和插入前三分之一（刚好）和红线一个=2和最后三分之一扣球：

在后面显示的比较中，该方案被描述为TQ3，而之前的更简单方案为TQ1。其他类型的方案可以作为本地的自动结果自适应时间步长，旨在使解算器在峰值之前更加稳定事件。这些已经过简单的研究，结果很有希望，并在中提到部分8.

应该注意的是，由于当前的时钟驱动模拟基础设施SpiNNaker，我们目前只能在模拟选择的时间网格上交付峰值已运行。因此，TQ只是试图避免滞后，否则会累积在状态变量轨迹，但无法使峰值时间本身增加虽然这在长时间的模拟中仍然是有益的。作为Hansel，马托、梅尼尔和内特纳(1998)指出除非有某种方法，否则高阶方法的通常优点将因该限制而降低可以实现插值峰值定时本身。这是一个就哪些神经模型和哪些解算器而言，值得进一步研究最易感；反复联系的人群可能特别敏感。正在研究的一种可能的解决方案是使用插值方法当一个尖峰信号被发送时，在其上贴上标签，从而恢复诸如使用较小的时间步长。

即使对所描述的方法进行了精简选择，也很容易生成难以总结的数据量。例如，选择的解算器方法如下Izhikevich模型中的四个参数允许行为范围广，有四种算术类型（32位：浮点浮动; 固定的指向累计，64位：浮点指向双重的; 固定的指向长期累积)至少两次例如，感兴趣的步骤，可以用不同的方式刺激神经元。为了避免在信息爆炸的情况下，已经做出了许多决定：

• 大多数人选择了一种参考神经元类型评估，其他两个示例的结果有限。Izhikevich描述作为常规峰值，并使用以下参数：⁠，⁠，毫伏，d日= 8.V（V）和U型在初始化毫伏和0，分别是。

• 每种情况下的参考结果都使用Mathematica中复杂的ODE解算器，自动使用自适应时间步长并透明地处理阈值交叉和刚度。我们已经执行了25个确保高精度的小数位数算法。

• 两个使用时间步长：大多数评估为1毫秒，有些评估为第条。

• 这两个32位算术类型的重点是：累计和浮动。如果64位类型(双重的和长堆积)展示特定行为明确地。

• 使用了两种不同的刺激。这个神经元最初不接受输入；则（1）在60 ms时，4.775 nA直流电流为交付并持续（直流输入）或（2）50 ms时，指数衰减输入时间常数为8ms（总充电量为80 pC）然后每50毫秒注射一次（突触输入）。一旦初始瞬时沉降，这些值提供约10 Hz的峰值频率，精确的脉冲间隔100 ms的间隔，然后可以很容易地与参考进行比较。直流输入为对ODE的精度性能进行更严格的测试解算器。

• 结果以图表的形式显示（1）薄膜电压V（V）相对于从50到350的参考毫秒和（2）n个^第个与前2个参考相比的峰值模拟的秒数（仅适用于较难的直流输入；所有方法都执行对于突触输入模拟来说非常好）。在膜电压图中重点应放在峰值时间和阈下行为的形状上与高度相反V（V）在扣球之前，因为后者主要是为生成图。

举个例子，RK-3（Heun）-ODE-TQ1-1000-累计指“Runge-Kutta 3rd带有Heun参数的订单作为算法ODE求解器实现，时间简单使用累积算术类型。”对于ODE解算器，使用静态内联函数的结果不会单独显示为它们应该只影响效率，而不是准确性。

选择用于调查的求解器方法如下：

• 欧拉伦格-库塔秩序泰勒级数顺序：这是最简单的显式解算器许多神经软件的可用性和默认选择包装。

• 伦格-库塔顺序：这是一个解算器家族，因此最常用的三个成员是被选中的。这一家族的名字在文献中经常被混淆，因此也会被混淆参见Lambert第4.3节(1973)，其中它们由参数化一₂如图所示在表中1我们将在粗体。

  Runge-Kutta三阶：另一类解算器，这次有两个自由参数。它们在文献中通常不会混淆，因此我们称之为他们的共同名字是：库塔、希恩、洛特金。正如这些不可避免的效率较低，我们没有像二阶模型那样彻底地研究它们方法。

Runge-Kutta二阶和三阶方法可以使用ESR或在标准ODE求解器框架。有很多理论与每个家族的哪个成员相关的方法应该是优越的，选择取决于这种优越性如何评估。例如，Ralston二阶方法和Lotkin三阶方法有它们的自由参数的选择是为了最小化局部截断误差。然而，在似乎存在实际算法的有限精度和其他实现细节这些理论结果不一定会转化为实际效益。

泰勒级数的二阶和三阶方法使用ESR和方程式实现来自阿尔萨斯的(2009)方法。

我们还简要调查了Chan和Tsai的方法(2010)将龙格-库塔方法与更高导数相结合，结果很有趣。二阶方法（我们称之为CT-2-ESR）产生与TS-2-ESR等式相同，因此是等效的。三阶方法是新的。

数字2和三给出了峰值滞后与第一个和选定的参考值的比较使用DC输入和两个32位类型的二阶方法(累计使用蓝线；浮动使用红线）。请注意，每个图都使用不同的年-缩放以提供最佳精度。

不出所料，RK1（=Euler=TS1）是表现最差的。在二阶方法中，RK-2（梯形）是一个明显的赢家累计和同样最好的浮动.明确使用浮动提高性能结束累计此处（以及所有后续结果中）。TS-2一般，但似乎更容易受到随机数字错误的影响围绕RK方法产生的线性滞后趋势的变化。在s、 改进了滞后性能（大约好45%）具有相同的近似顺序累计。使用浮动，除了RK-1没有那么好。

尽管RK-2（Ralston）被设计为最小化局部截断误差比其他两种RK-2算法更精确累计和相同其他三种类型。由于速度较慢，我们还没有显示这些结果，并将重点放在这里显示的ESR算法的ODE版本都是对于累计通常与另一个相同类型。

接下来的两组图形4和5、和6和7（每种刺激类型输入一个），显示来自膜电压（即状态变量V（V）)参照50至350毫秒。每组中的第一个以1ms时间步长比较不同的解算器，并与累计类型。第二个显示了改变时间步长和32位类型，带RK-2（梯形）-ESR解算器。

可以在膜电压中看到前面的图中测量的滞后，但超过使用的时间尺度短，方法之间的分离不明确。RK-1显示扣球项目后射门不足。这是因为欧拉很简单术语为线性外推，因此必然会超出状态变量轨迹包含一个实质性的二阶分量，其方向相反方向——峰值后的初始恢复将用于膜电压状态大多数集成和火灾模型中的变量。

减少时间步长和浮点的比较优势清楚地显示在图5.利益浮动在1ms是清晰的（橙色与红色），以及浮动在s很难与参考。有趣的比较是浮动在1 ms和累计在s（橙色对比绿色）。前者可能更快，但性能明显更好滞后曲线表明这一点是显而易见的。

突触输入是一种更容易的测试，滞后性能被视为非常相似方法之间。图中的峰值事件后，RK-1仍显示超调6各种方法在150和250 ms时遵循脉冲具有不同程度的保真度。

使用RK-2（梯形），较短时间步长和浮点的优点不太明显图中的突触输入7.

为了进行比较，尽管有额外的不可避免地会产生计算负载。（见图8和9.)

看到标准三阶方法在1 ms时间步长，即使是64位类型。在Runge-Kutta三阶方法中，只有RK-3（Heun）提供了可靠的直流输入解决方案，即使如此，它们也较少比二阶方法更精确。使用累积，仅ESR方法完全奏效了。以往使用ODE解算器的经验显示出了一个渐进的改进根据中的分析结果预测，阶数增加到五阶和六阶标准文本。然而，这种体验通常是通过自适应时间步长实现的算法和双精度浮点。

因此，需要进行一些诊断，以调查以下初步假设：下溢或溢出（或两者兼而有之），为此，制定了一个范围检查宏可以打印选定状态变量的值或内部的中间值算法。此宏适用于所有算术类型，可以在一个位置删除以进行速度测试。这表明高阶RK方法、步骤内预测和在（和特别是接近尾声时）在尖峰事件发生之前，时间步长可能会变得很大。高阶RK方法更容易受到这种现象的影响的一个原因是与低阶方法相比，在时间步长中进行预测或计算导数，这些影响在峰值附近迅速变得更加普遍。绘制高位在这些情况下，参考溶液的导数显示出非常大的值，这将溢出算术类型或违反评估局部错误时所做的假设存在精确算术的行为。本质上，膜电压表现为截止点附近的超指数行为，当相对较大解中使用了固定的时间步长，而高阶方法因其试图追踪这些非常大的高阶衍生品。使用的时间步长，没有发现任何问题算法。

当接近峰值时，本地切换到较小时间步长的方法进行了简要的调查，并显示了有希望的结果。它们明显提高了这里描述的三阶方法具有较小的速度损失，并提供了一个机会如果公式正确，则可以获得比TQ3更好的阈值交叉时间估计。它可能值得进一步探讨这些想法，因为它们对所有人都有一定的好处如果速度惩罚可以做得足够小的话。

虽然Runge-Kutta方法被广泛认为是一种可供选择的方法解中可能会遇到奇点（例如，参见Press等人。，1992)，一些人提到了在这种情况下应用高阶版本的危险（Gear第232页，1971以及恩赖特、海姆、奥伦和夏普的第21页，1995). 我们的结论是，结合固定时间步长的限制太大，无法有效跟踪接近峰值的超指数行为，在较小程度上，有限精度问题，是高阶方法在这里如此无效的原因。

以数字表示10至12，我们展示了一些膜电压与RK-1和RK-2（梯形）参考结果的比较求解器（后者具有不同的时间步长和算术类型）Izhikevich神经元使用直流输入。

所选类型包括：

1. 喋喋不休{⁠，⁠，⁠，d日=2｝

2. 快速峰值{⁠，⁠，⁠，d日=2}

对于ODE解算器来说，这两种神经元类型都将更难处理因为相对于固定时间步长，变化率更高。对于喋喋不休神经元，这些高变化率成簇出现；对于快速峰值神经元，对于给定输入。

使用喋喋不休很明显，要评估它要困难得多生成尖峰簇时尖峰行为的详细信息。图10放大90 ms到120 ms之间的时间以提高第一个集群的清晰度。只有RK-2（梯形）位于s的时间步长与四个峰值相匹配由参考解算器产生，尽管有一些滞后。1毫秒解算器均产生三个峰值，RK-2（梯形）浮动明显领先，RK-1累积明显滞后。图11在通常的时间轴上显示了一些模式存在于与参考神经元不同的滞后行为中。尽管最严重的滞后是RK-1累积，但RK-2（梯形）的时间步长并没有太好，似乎是随着时间的推移，情况越来越糟，这是一个惊喜，因为它与参考很好地配合神经元。此处还显示了RK-2（梯形）浮子领先其他浮子的趋势，但它似乎领先太多，至少在前350毫秒。也许令人惊讶的是在整体滞后方面表现最好的似乎是RK-2（梯形）-在1ms时间内累积步骤。

对于快速峰值图中的神经元12，故事简单了一点。所有1 ms时间步长结果均滞后于RK-1最差的参考结果，然后是RK-2（梯形）浮动和RK-2非常相似。RK-2（梯形）-累计s支架在第十二次峰值时，相隔距离远小于10毫秒。显然，时间步长较小更能处理较高的峰值率，这并不奇怪。

在所有测试的案例中，64位类型产生的结果与32位浮动类型，因此没有生成单独的图形。这个表明在这些情况下不会产生显著的效益，因此64位类型的计算和存储成本目前还不合理。

值得注意的是，对于突触输入（从表面上看是更多的现实案例），使用累计1 ms没有太多即使对于更简单的方法，也会受到相应的惩罚。然而，尽管真实网络可能正在接收来自（在某些情况下，在听觉系统中）a的输入单突触，更常见的情况是多达1000个或更多突触，而对于后者在这种情况下，这些输入在实际尖峰速率和1ms时间步长下的总和将在含有高斯噪声的直流输入下渐近收敛。这使得情况变得更糟在模拟实际行为方面挑战输入类型。缅因州和塞诺夫斯基(1995)还有Brette和Guigon(2003)表明DC周围的一些非周期变化在未来的工作中，组件可能是一个良好的确定性测试信号，模拟真实的伪随机输入噪声量并分析输出存在这种噪声输入的分布将更接近实际应用案例。

显然，要找到最佳的ODE解算器，需要使用一组现实的神经网络行为。从上一节中对其他两种常见类型的简要分析来看似乎一个能很好地处理一种神经元类型的解算器并不一定能传递这种信息性能优于其他类型。虽然作为另一个演示很有用性能特征，应该记住，其他两种神经元类型是更难的挑战，这些峰值频率在峰值神经网络中不太可能出现这与大部分大脑中的现实行为相匹配（尽管它们可以暂时存在或在某些区域）。无论如何，神经元类型、时间步长、，算术类型、解算器算法，甚至算术运算的精确顺序都使在这里很难达成明确的声明。然而，存在一些模式。这个在某些情况下，二阶方法通常会显著改进RK-1。更多完整的调查将需要在以便更彻底地搜索各种参数空间。

一些峰值滞后曲线的几乎线性性质可能会让人相信可以进行简单的校正，以便将此行恢复到x个-轴。然而，观测到的滞后是一个复杂的非线性函数输入动力学和神经元定义等等，解决这个问题在计算上比简单地使用更好的解算器更具挑战性。

表中的结果2显示一个SpiNNaker板以200 MHz的频率运行，在一个简单的循环中执行10000个神经元更新，除以10000。突触输入方案使用1ms时间步长。在初始值之后瞬态，峰值速率10赫兹，所以离散状态更新的相对比例正确。正时修正为消除生成输入、收集和计数峰值所需的代码开销，和其他行政职能。因此，每次更新的时间仅指每种算法特定的计算、内存访问和移动。编译器标志-使用Ofast，禁用缩略图代码以提高速度，而不是其他考虑因素允许编译器执行尽可能多的优化。²但是，不使用手动调整装配代码初步调查表明，这有可能提高性能例如，在不同定点类型的乘法可能传递提高了准确性。ARM968处理器上的浮点操作由执行软件库功能。

对于ODE解算器结果，我们对相同的方法进行了大量比较：（1） 实现为标头中的普通C源或静态内联函数文件（表示这些文件S.I.公司。)（2）最大使用函数指针的通用ODE驱动程序(发电机驱动器)而且是最多的一种专门的算法，它显式地对算法进行编码，并通过知道这一点来避免内部循环存在两个状态变量(规范驱动程序).CompErr公司意味着编译器错误不允许生成完成。

以RK-1为参考，目前仅考虑累计结果，对于ESR，二阶方法的额外成本在54%到91%之间，介于静态内联形式的ODE为59%和143%，使用ODE的ODE在191%和323%之间状态变量上的函数指针和通用循环。显然，ESR是有效的在这里以及在准确性方面。静态内联结果表明编译器有效利用跨源文件允许的优化否则是可能的。最快的二阶方法是RK-2（中点）-ESR。这两个最快的三阶方法成本分别为125%和168%。

对于浮动，相对成本有所不同。平均而言，解决方案所需时间延长6.6倍，但对于一般ODE，相对惩罚一致较小结果，可能是因为他们一开始就更糟糕了。更快的二阶与RK-1相比，方法通常稍慢，其中的顺序为大致相同，但TS-2-ESR现在是71%的最小处罚，尽管RK-2（中点）-ESR在75%时并不落后。

增加TQ3修正成本通常约为0.025至0.05s.这是可能比预期的要大，因为修正值仅在峰值出现时计算大约每100次更新发生一次，这表明每个TQ3的成本修正值约为5第条。完全消除TQ校正可以提高效率，原因有二：（1）由于固定的时间步长现在总是预先知道的，这简化了最佳ESR形式，导致更少的计算，（2）一个间接内存访问不再是必修的。例如，RK-2（MidPoint）-ESR的一个神经元更新成本降低到0.73s用于累计，提高了15%。随着算法变得越来越简单，达到了另一个阈值，这样的增益变得均匀更大，可能是因为编译器发现优化使用寄存器。例如，RK-1的剥离版本累计哪里删除所有与TQ相关的代码。每次更新的时间变为0.45 s、 比标准TQ1版本提高20%。

也许有了像SpiNNaker这样的创新技术结合了非常低的能耗和实时性能，如果需要做出一些权衡。与使用的神经求解器框架的比较双精度浮点算法，自适应时间步长，高阶/隐式求解器和复杂的阈值交叉算法将不可避免地显示一些状态变量轨迹和峰值时间方面的差异。然而，使用仔细选择求解器和TQ方法，甚至只使用浮动而不是累计每个核的神经元更少（或者，等价地，比实时慢的性能）将提供非常接近的性能这些耗能更高的系统，它们本身只提供了一小部分实时性能。其他可能性，如订单自适应性和步长内适应性在下文中进行了讨论。就像数学领域的几乎所有东西一样计算，需要在{准确性、速度}空间中做出明智的选择我们希望这项研究有助于阐明一些权衡。

Izhikevich神经元模型创造了超高的膜电压轨迹尖峰附近的指数具有较大的高阶导数，尖峰事件为不连续性。这些将始终对明确的固定时间步长提出挑战具有相对于变化率显著的时间步长的算法跟踪突然增加和不可避免的滞后。这是所示滞后的主要原因在精度结果上，高阶方法的失败在于1ms时间步长下的低阶。

尽管本研究已排除了全时间步长自适应性，但还有其他一些选择是可能的。Parker-Sochacki提供了一种订单自适应形式，添加了高阶项如有必要，本文在添加在我们的例子中，泰勒级数解的高阶项高达四阶。这很容易在这种情况下，由于状态变量的每个高阶更新可以很简单添加到上一个。目前的结果很有希望，但变化很大，因此还没有还没有详细报道。很可能决定何时添加高阶项这种方法成功的关键，这需要进一步研究。也有可能在存在数值的情况下，高阶项将更难计算错误。Runge-Kutta解算器也可以使用类似的方法，但这是不可避免的计算成本更高，或者至少复杂得多。

另一种可能性是有限的时间步长自适应性点，解算器将其内部时间步长除以一个整数d日和然后重复求解d日将更新的神经元返回到系统。这显然会创建一个“d日时间“开销，但仅限于点在需要的地方，衍生工具变得非常大，因此不太可能达到峰值为整个模拟产生了显著的开销。与订单自适应性一样关键是要对何时实施做出有效而稳健的决定。首字母实验显示以相对适中的成本对所有解算器都有有用的改进，使三阶方法更加稳定。

一种简单的机制，用于消除由于忽略提出并实现了阈值交叉。最简单的TQ1几乎没有开销并删除时间步长偏差。一个更复杂的TQ3补充道只需少量额外开销，并提供更精确的峰值后过渡平均绝对偏差为时间步长。允许时间步长在峰值附近自适应，以获得更好的准确性和如上所述，稳定性为测量基本上无需额外成本即可实现阈值交叉时间。然而，目前这些机制将在时间步长之间产生峰值。随着SpiNNaker系统的发展可能是某些基本规则发生了变化；例如，事件驱动而不是时钟驱动行为可能成为可能（参见Brette等人的第2.3节。，2007).

准确性和效率结果都对算法。显然，编译器根据优化重新排序指令方案。在这项工作中，人们普遍发现，一个在总共有80个数学运算的算法中，减法可以使显著差异，尤其是在精确度和32位类型方面。还有，最整洁的代数约简并不总是能够产生最佳精度（累计尤其是），但它们通常会产生最佳效率。显然，特定解决方案的优化也需要艺术和实验作为科学。（见马特尔，2009; Darulova等人。，2013; 和Gao等人。，2013，了解如何使该过程更加科学。）

其他神经元模型将提供不同的挑战和机遇。复杂的行为（例如，霍奇金-赫胥黎模型）几乎肯定需要更高阶的方法和可能是隐式解算器，对神经元更新性能和每个核可以实时支持的神经元。有效应用该车型的Parker-Sochacki如Stewart和Bair所示(2009)这种方法可能会引起更广泛的兴趣。

当然，如果持有一个，则上述所需的额外计算可能是不合理的在时间步长上固定状态变量，以解决其他引入的更多错误而不是精确的解决方案。求解耦合的二维积分系统两者的方程式V（V）和U型是真正的答案，但确实如此到目前为止，证明这是难以捉摸的，而且可能是不可能的。几乎可以肯定的是，如果这样的解决方案存在时，它需要使用超越函数和复数加强了对定点类型快速准确的此类操作库的需求可用。

如果ODE是自治的，并且可以使用显式解算器。我们已经证明，仔细实施的算法可以同时比同类优化程度最高的通用求解器更精确、更快使用静态内联函数优化。神经元解算器代码可以交付的预测试和标准编译库+表头形式也提供软件工程优势。

由于时间和空间的限制，我们只看了一小部分Izhikevich本研究中的参数（以及因此产生的神经元行为），因此不能将其视为详尽分析。我们希望我们选择的例子显示出可以推断的趋势得出更一般的结论，但当然也有可能是看不到的病理病例存在。即使是伊兹凯维奇指南针范围内其他两种神经元类型的小样本参数化显示了与参考神经元结果不同的模式。一个实现更好覆盖的方法是自动化测试机制。这部分是在SpiNNaker板上运行的主C事件循环中执行了研究，但为了全面实现这一点，在主机上编写测试机制脚本的能力将是理想。目前，使用定点类型并不简单，因为它们只有交叉编译到没有操作系统的ARM目标。如果这些类型可以在主机上支持，或者在可编写脚本的模拟器中支持尤其是在更大的示例空间上进行环境验证时更容易的。它还将有助于使用随机输入进行测试实际用例，导致基于使用例如Kullback-Leibler信息理论的峰值或峰值间隔测量，参考神经元提供参考输出分布。

存在多步骤、多阶段和多验证方法之间的混合，并且可以推理使用布彻第5章中所示的一般线性方法(2003)，他提到了一些更知名的示例：Rosenbrock和Obreshkov方法。这些混合方法可以在多大程度上在保持稳定性和保持在固定1ms时间步长和累计算术类型尚未得到充分的探索，但其潜力是不言而喻的结合多阶段和多衍生方法（通常称为Obreshkov或Turan方法）。Chan和Tsai中讨论的方法(2010)应该值得进一步研究，因为它们允许代数将泰勒级数和龙格库塔解结合起来，以获得更好的精度以最小的额外开销。初步测试非常有希望，CT-3-ESR正在进行中目前最精确的算法为1ms时间步长（除了累计类型，其中数值下溢破坏了迄今为止所有算法排序的准确性同时也是浮点类型最快的三阶算法。（有关工作，请参阅Shintani，1972; 三井，1982; 小野和吉田，2004.)

目前，根据迄今为止对Izhikevich神经元类型，速度和准确度之间的最佳权衡-使用累计类型和1ms时间步长，以便于实时操作是RK2-ESR求解方法之一；中点对于速度或梯形以确保准确性。

为了指示此表中状态变量和中间变量使用的典型值图中的直方图13显示2秒的样本使用1ms时间步长和直流输入和TQ3模拟参考神经元。

实现信中描述的几个ESR解算器核心的C代码是可在获取http://www.mitpressjournals.org/doi/suppl/10.1162/NECO_a_00772.

我们感谢Dave Lester和Jamie Knight讨论定点问题和静态内联功能优化；塞尔吉奥·戴维斯、路易斯·普莱纳、史蒂夫·坦普尔和西蒙·戴维森通过SpiNNaker体系结构和开发环境获得初步指导对于M.H。；吉姆·加赛德（Jim Garside）关注装配输出，以确定优化机会；Patrick Camilleri为LaTeX提供帮助；以及审阅者的有用意见和参考文献。

导致这些结果的研究得到了EPSRC（英国工程物理科学研究委员会授予EP/G015740/1和EP/G0115775/1），ARM有限公司欧洲联盟第七框架计划下的欧洲研究理事会（FP7/2007-2013）/ERC拨款协议320689和欧盟旗舰人脑项目（FP7-604102）。