问题标签 [numerical-stability]

这是一个指数时间差分法的算法，使用牛津的原始 Matlab 代码

Python代码

Python 代码运行在某种程度上表明 NaN 有一些价值，但对于 matlab 却没有。完全相同的代码

所以我想知道为什么以及发生了什么，你可以自己运行代码。

更让我觉得奇怪的是，第一次迭代，两者的搭配非常好！

背景：我正在模拟来自Cox 比例风险模型的数据。Bender et 提供了一种方法。人（2005 年）。在这篇文章中，它在 R 中得到了很好的解释和实现。Cox 模型的复杂性对这个问题并不重要。

问题：我需要评估 exp(-1 * X * beta)。由于我正在处理的问题的一些细节，X 或 beta 需要相对较大，介于 50 到 100 之间。我相信这会导致准确性问题。

1.0000000 0.9048374 0.8187308 0.7408182 0.6703200 0.6065307 0.5488116 0.4965853 0.4493290 0.4065697 0.3678794

1.000000E+00 6.737947E-03 4.539993E-05 3.059023E-07 2.061154E-09 1.388794E-11 9.357623E-14 6.305117E-16 [9] 24 8.756511e-27 5.900091e-29 3.975450e-31 2.678637e-33 [17] 1.804851e-35 1.216099e-37 8.194013e-40 5.521082e-42 3.720076e-44

我不相信最后的价值观是正确的。R 真的将 exp(x) 计算为 44 位数字吗？最终，我将这些值插入一个函数中，我想我得到了无意义的答案。我是对 R 缺乏信心还是我是正确的。如果这确实是一个问题，我能做些什么吗？

我考虑过只使用 e^x 的泰勒级数表示的前 2-3 项。我可以毫无问题地做到这一点。我不知道这是否能解决我的问题，因为我认为 R 在幕后做了类似的事情，

归结为我的问题，即使参数很大，我可以从 exp() 函数中获得精确的值吗？

假设我用 c 写了一个表达式，例如

这里 a,b,c,d,e,f 都是双精度变量。例如，当我使用带有一些优化设置的 gcc 编译器编译代码时，编译器是否尊重我编写的表达式的特定形式，还是修改表达式以使代码运行得更快？例如，带有 -O2 优化设置的 gcc 是否可以将上述表达式编译为

或者它会保持表达式不变？我担心编译器会改变我的方程式的形式，这可能会影响我的表达式的数值稳定性。

为了学习深度学习神经网络的更精细细节，我用自制的所有东西（优化器、层、激活、成本函数）编写了自己的库。

在 MNIST 数据集上进行基准测试并仅使用 sigmoid 激活函数时，它似乎工作正常。

不幸的是，在用relus替换这些时我似乎遇到了问题。

这是我在约 500 个示例的训练数据集上 50 个 epoch 的学习曲线：

前 8 个时期一切都很好，然后我在虚拟分类器的分数上完全崩溃（~0.1 准确度）。我检查了relu的代码，看起来还不错。这是我的向前和向后传球：

罪魁祸首似乎在于 relu 的数值稳定性。我为相同的结果尝试了不同的学习率和许多参数初始化器。Tanh并sigmoid正常工作。这是一个已知的问题？relu它是函数非连续导数的结果吗？

我有一个功能：

对于那些没有看过的人来说，这比后者更可取，x0 + (x1-x0) * alpha因为后者并不能保证lerp(1.0f, x0, x1) == x1.

现在，我希望我的lerp函数有一个额外的属性：我想要lerp(alpha, x0, x1) == lerp(1-alpha, x1, x0). （至于为什么：这是一个更复杂功能的玩具示例。）我想出的似乎可行的解决方案是

这种双重减法具有接近零和接近一的舍入效果，因此如果alpha = std::nextafter(0)(1.4012985e-45), then 1 - alpha == 1and so 1 - (1-alpha) == 0。据我所知，这始终是正确的1.0f - x == 1.0f - (1.0f - (1.0f - x))。似乎也有这样的效果w0 + w1 == 1.0f。

问题：

1. 这是一个合理的方法吗？
2. 我可以相信我的编译器会做我想做的事吗？特别是，我知道在 Windows 上它有时对部分结果使用更高的精度，而且我知道编译器可以做一些代数；显然 1-(1-x)==x 代数。

这是在 C++11 中使用 Clang、VisualStudio 和 gcc。

我想计算, , 和mean使用std一次性skewness算法。我发现的最简单和最快的一种方法是由Berkeley Research Group 的 Stuart McCrary发表的。例如，一个人可能会使用：kurtosiscovariancestd

我读到这种方法不够好，因为它在数值上不稳定。不幸的是，我对数值稳定性没有深入的了解，但据我所知，这是一些问题，这是由于浮点运算的精度有限而发生的。

就我而言，我将只处理范围内的整数10^1-10^6。

我可以在我的情况下使用这种方法并且不关心数值稳定性吗？

我正在使用 CVXPY（1.0 版）求解二次规划（QP），我经常遇到这个异常：

SolverError：求解器“xxx”失败。尝试另一个求解器。

这使我的程序非常脆弱。我尝试过不同的求解器，包括 CVXOPT、OSQP、ECOS、ECOS_BB、SCS。他们都有或多或少相同的问题。我注意到，当我使求解器的停止标准更严格（例如，降低绝对误差容限）时，我得到SolverError的频率更高，而当我不那么严格时，SolverError问题会减弱甚至消失。我还发现 CVXPY 抛出的方式SolverError是随机的：如果我多次运行同一个程序，有些运行会得到最佳结果，SolverError而另一些会得到最佳结果。

虽然我可以通过尝试更多次并降低停止标准来避免 SolverError，但我真的很想了解异常背后的真正具体原因

SolverError：求解器“xxx”失败。尝试另一个求解器。

这个错误并没有真正提供信息，我不知道如何提高解决问题的稳健性。其原因是否特定于求解器？是否会针对一组明确定义的情况引发此异常？或者它只是一种说“由于未知原因出现问题”的方式？这些可能是什么原因？

当使用 IEEE 754 浮点数 a 和 b 进行评估时，就 (a - b) + b 之和的 a 和 b 的大小而言，最坏情况下的误差是多少？我可以期望它有多接近？

我想以数值方式评估线性生死过程的转移概率

其中是二项式系数和

对于大多数参数组合，我能够以可接受的数值误差（使用对数和 Kahan-Neumaier 求和算法）对其进行评估。

当加数在符号上交替出现并且数值误差占总和时会出现问题（在这种情况下，条件数趋于无穷大）。这发生在

例如，我在评估p(1000, 2158, 72.78045, 0.02, 0.01). 它应该是 0 但我得到的值非常大log(p) ≈ 99.05811，这对于概率来说是不可能的。

我尝试以许多不同的方式重构总和，并使用各种“精确”求和算法，例如Zhu-Hayes。我总是得到大致相同的错误值，这让我认为问题不在于我对数字求和的方式，而是每个加数的内部表示。

由于二项式系数，值很容易溢出。我尝试进行线性变换，以使每个（绝对）元素保持在最低正态数和 1 之间的总和中。它没有帮助，我认为这是因为许多类似量级的代数运算。

我现在处于死胡同，不知道如何进行。我可以使用任意精度的算术库，但是对于我的马尔可夫链蒙特卡罗应用程序来说，计算成本太高了。

当我们无法在 IEEE-754 双精度中以足够好的精度存储部分总和时，是否有适当的方法或技巧来评估这些总和？

这是一个基本的工作示例，我仅按最大值重新调整值并使用 Kahan 求和算法求和。显然，大多数值最终都是 Float64 的次正规。

我正在尝试使用 Python3 用欧拉法求解这个微分方程：

根据 Wolfram Alpha 的说法，这是正确方程的图。

同样，根据 Wolfram Alpha，在这种情况下，经典的欧拉方法不应该是稳定的，正如您在区间结束时所看到的：

但是，在我的实现中，欧拉方法提供了一个稳定的结果，这很奇怪。我想知道我的实现由于某种原因是错误的。尽管如此，我找不到错误。

我生成了一些点和一个图，比较了我的近似值和函数的分析输出。蓝色为对照组的分析结果。红色，我的实现的输出：

那是我的代码：

谢谢您的帮助。

==================================================== ==========

更新

在@SourabhBhat 的帮助下，我能够看到我的实现实际上是正确的。确实，它产生了不稳定。除了增加步长之外，我还需要进行一些放大才能看到它的发生。

下图不言自明（步长为 0.22）：