问题标签 [numerical-stability]

我正在通过此维基百科链接在 Matlab 中实现 Baum-Welch 算法：Baum-Welch algorithm。这是我在金融时间序列中预测波动性的一部分。
我有两个问题：

1：在维基百科页面的最后一个更新步骤中，有人告诉我们“这些步骤现在迭代地重复，直到达到所需的收敛水平。”。
那么如何声明一个条件来停止循环呢？此外，应该评估哪些变量以查看它们是否可以接受？

2：如果您注意维基百科的公式kesi：

alpha分子 ( and beta) 和分母 (just) 中有两个因子被缩放alpha。所以他们不会取消彼此的效果。我已经实现了方程（例如，在一个重复 20 次的循环中）并且我已经完成了缩放过程。但是“转移概率矩阵”和“初始分布”和“排放矩阵”得到了NaN值！！！

我已经阅读了这个问题的答案Baum-Welch 许多可能的观察结果。我已经根据那里提到的教程进行了缩放。
这是我的代码：

那么我现在应该为这个缩放问题做些什么呢？这个NaN值是因为缩放问题还是其他原因？
谢谢你的时间。

我编写了一个 Python 3 包，可以集成常微分方程系统。在写这篇文章时，我遇到了一个我似乎找不到任何信息的问题。

问题：

基本上我有一个自适应算法，可以根据局部误差调整步长（实际方法是使用 Cash-Karp 系数的 RKF 算法的实现），并且在某个点之后，如果步长足够大，该方法将超过结束时间。

我的解决方案：

我已经设置好了，如果当前步长会导致积分超过结束时间，那么我将步长设置为剩余时间值的一半。方程形式：dt = (t_end - t_current)/2

但这对自适应算法中使用的相对误差有一定的要求，会导致时间步长变得微不足道。大约在 1e-15 的数量级上，这会导致数值不稳定，因为由于机器精度的限制，t_current 永远不会增加。

我的问题：

我想知道是否有更有效的方法来实现积分算法，以便结束时间相对于预期结束时间精确到小数点后 15 位，而不会产生这种数值不稳定性。

我目前正在根据一篇科学论文实施一种用于 3D 点云过滤的算法。

在计算特定值的旋转矩阵时，我遇到了一些问题。目标是将点旋转到由法线向量（Z 轴）的方向定义的坐标系中。由于以下查询在 X、Y 轴上旋转对称，因此这些轴的方向无关紧要。R 定义如下：Rotationmatrix

n 被归一化。当 n_z 变得非常小或为零时，就会出现问题。因此，我考虑在计算第 2 行的叉积之前对第 1 行进行归一化。

然而，行列式变为-1。旋转矩阵会导致正确的结果吗？R 是正交的，但 det|R| 不 +1

感谢您的任何建议

我必须对具有以下形式的 ODE 系统进行数值求解：

其中u_j(t)和v_j(t)是时间 的复值标量函数t，f_i和g_i是给定的函数， 和j = -N,..N。这是一个初值问题，任务是在某个时间找到解决方案T。

如果，则可以独立求解g_i(t) = h_i(t) = 0不同 值的方程。j在这种情况下，我借助四阶 Runge-Kutta 方法获得了稳定且准确的解。然而，一旦我打开耦合，结果就时间网格步长和函数的显式形式而言变得非常不稳定g_i，h_i。

我想尝试使用隐式 Runge-Kutta 方案是合理的，在这种情况下它可能是稳定的，但如果我这样做，我将不得不评估一个巨大的 size 矩阵的逆矩阵4*N*c，其中c取决于顺序每个步骤的方法（例如c = 3，对于 Gauss-Legendre 方法）。当然，该矩阵将主要包含零并具有块三对角形式，但它似乎仍然非常耗时。

所以我有两个问题：

1. g_i有没有一种稳定的显式方法，即使在耦合函数h_i（非常）大的情况下也能工作？

2. 如果隐式方法确实是一个很好的解决方案，那么最快的块三对角矩阵求逆方法是什么？目前我只是执行一个简单的高斯方法，避免由于矩阵的特定结构而出现的冗余操作。

可能对我们有帮助的其他信息和详细信息：

• 我使用 Fortran 95。

• 我目前考虑g_1(t) = h_1(t) = g_2(t) = h_2(t) = -iAF(t)sin(omega*t)，其中i是虚数单位，A并且omega给定常数，并且F(t)是一个平滑的包络线，首先从 0 到 1，然后从 1 到 0，所以F(0) = F(T) = 0。

• 最初u_j = v_j = 0除非j = 0. 绝对值很大的函数u_j和对所有 来说都非常小，因此初始峰值不会到达“边界”。v_jjt

我如何评估 xe^x/(e^x-1) 与零附近的数值稳定性以及当 x 非常积极或消极？我可以访问 和 中的所有常用数学numpy函数scipy。

对于接近 0 的 x 值，有没有办法在 python 中正确计算 log(1+x)/x 的值？当我通常使用 np.log1p(x)/x 执行此操作时，我得到 1。当我使用 np.log(x) 时，我似乎得到了正确的值。log1p 不应该更稳定吗？

IEEE 754（IEC 559）浮点标准中的乘法是否x + x可以互换加法，或者更一般地说，是否有任何保证并且总是给出完全相同的结果？2 * xcase_addcase_mul

Is it guaranteed by the C++ standard that angle == std::acos(std::cos(angle)) if angle is in the range [0, Pi], or in other words is it possible to restore the exact original value of angle from the result of std::cos using std::acos given the mentioned range limit?

The marginal cases when angle is infinity or NaN are omitted.

我正在尝试Ax = b在 MATLAB 中求解系统，其中A是一个 30x30 三角矩阵，（非零）值范围从1e-14到0.7，并且b是一个 30 元素列向量，值范围从1e-3到1e3。当我输入时x = A\b，我得到一个答案并且没有警告消息，但答案是不合理的（看起来只是向量底部的随机数）。我认为这是由于数字错误。

此页面上的消息 5 建议分解/缩放矩阵以避免数值错误，但我无法弄清楚如何计算缩放矩阵。

所以问题是：这确实是数值不稳定性的一个例子，如果是这样，我怎样才能重新调整我的矩阵 A，或者改变 MATLAB 执行计算的方式，以避免它？

这是产生问题的矩阵和向量：

可以在此处找到具有全精度变量的 .mat 文件。

这是我在我的机器上得到的结果（OS X 10.10.5 上的 Matlab R2013a）：

虽然norm(A*x-b)返回一个 1e-13 的值，但考虑到我要解决的问题，结果在物理上并不合理（x 中的值应该是单调递减的，并且没有一个应该是负数）。例如，这是一个类似的数据集，它返回具有相同矩阵 A 的正确（看起来）解决方案：

这是另一个问题/答案的分支

我想要一个与此等效的功能：

这是稳定和健壮的，即它不会溢出小的值tau，也不会溢出大的值x。由于这将用于计算概率，因此输出的总和应为 1。

换句话说，我传入了一些值（和一个温度），并且我想要一个与输入和 tau 一起“缩放”的概率数组作为输出。

例子：

因此，当 tau 趋向于 0 时，输出中的最高概率位于最高元素的位置。随着 tau 变大，概率变得越来越接近。

（可选）有关链接答案的问题。在那里，尼尔给出了以下选择：

然而，这个输出总和不等于 1，解释是函数返回一个只有N-1自由参数的分类分布，最后一个是1 - sum(others)。但是在运行时，我注意到对于长度为 3 的向量，它返回一个长度为 3 的向量。那么缺少的向量在哪里？我可以使它等同于上面的例子吗？

为什么这个答案是稳定的？一个人如何从这个简单的公式得到softmax这个？

可能相关的问题： General softmax but without temperature