问题标签 [numerical-stability]

我正在尝试为 Python 编写 Wilson 的谱密度分解算法 [1] 的实现。该算法迭代地将 [QxQ] 矩阵函数分解为其平方根（它是用于谱密度矩阵的 Newton-Raphson 平方根查找器的一种扩展）。

问题是我的实现只收敛于大小为 45x45 或更小的矩阵。所以经过 20 次迭代后，矩阵之间的平方和差约为 2.45e-13。但是，如果我输入大小为 46x46 的输入，它直到第 100 次左右迭代才会收敛。对于 47x47 或更大，矩阵永远不会收敛；误差在 100 到 1000 之间波动大约 100 次迭代，然后开始快速增长。

你将如何尝试调试这样的东西？似乎没有任何具体的疯狂点，而且矩阵太大，我无法真正尝试手动进行计算。有没有人有提示/教程/启发式来找到像这样奇怪的数字错误？

我以前从未处理过这样的事情，但我希望你们中的一些人...

谢谢， - 丹

[1] GT威尔逊。“矩阵谱密度的分解”。SIAM J. 应用程序。数学（第 23 卷，第 4 期，1972 年 12 月）

我正在尝试对共现矩阵 (C) 进行因子分析，该矩阵是从术语文档矩阵 (TD) 计算得出的，如下所示： C=TD*TD'

理论上 C 应该是半正定的，但事实并非如此，因此因子分析算法无法使用它。由于速度原因，我无法更改算法）。

我查了一下，这可能是一个数值稳定性问题： 一种用于生成半正定矩阵的简单算法 - 答案 2。

在这里进行的好方法是什么？

许多数值算法倾向于在 32/64 位浮点上运行。

但是，如果您可以使用精度较低（且功耗较低）的协处理器会怎样？那么如何在数值算法中使用呢？

有谁知道解决这些问题的好书/文章？

谢谢！

我有一个用 Python 编写的贝叶斯分类器，问题是当我将特征概率相乘时，我得到非常小的浮点值，例如 2.5e-320 或类似的值，然后突然变成 0.0。0.0 显然对我没有用，因为我必须根据哪个类返回 MAX 值（更大的值）来找到“最佳”类。

处理这个问题的最佳方法是什么？我考虑找到数字的指数部分（-320），如果它太低，则将该值乘以 1e20 或类似的值。但也许有更好的方法？

我正在寻找一种算法来计算 ln(1-x)。x 通常很小 (<0.01)，但有时它可能会更大。算法需要准确，而且不能太慢。我宁愿不使用 ln(x) 的库，因为我可能会失去准确性。

我正在尝试在 Java 中实现模糊 C-Means 算法的一个版本，并且我试图通过只计算一次可以计算一次的所有内容来进行一些优化。

这是一个迭代算法，关于矩阵的更新，像素 x 簇成员矩阵U（一行中的值之和必须为 1.0），这是我要优化的更新规则：

其中 x 是矩阵的元素X（像素 x 特征），v 属于矩阵V （簇 x 特征）。是一个参数，m范围从1.1到，是簇的数量。使用的距离是欧几里得范数。infinityc

如果我不得不以一种平庸的方式实现这个公式，我会这样做：

通过这种方式，一些优化已经完成，我预先计算了 和 之间所有可能的平方距离，X并将V它们存储在一个矩阵D中，但这还不够，因为我循环遍历元素V两次，导致两个嵌套循环。查看公式，分数的分子与总和无关，因此我可以独立计算分子和分母，并且每个像素只需计算一次分母。所以我找到了这样的解决方案：

在计算距离时，我将分母预先计算到矩阵的另一列中D：

通过这样做，我遇到了“平庸”计算和第二个计算之间的数值差异，导致不同的数值U（不是在第一次迭代中，但很快）。我想问题在于将非常小的数字取幂到高值（ 的元素U可以从 0.0 到 1.0 并且exp，对于m = 1.1，是10）导致非常小的值，而通过将分子和分母相除然后取幂，结果似乎在数值上更好。问题是它涉及更多的操作。

更新

我得到的一些价值观ITERATION 0：

D这是未优化的矩阵的第一行：

384.6632 44482.727 17379.088 1245.4205

这是矩阵的第一行D优化方式（注意最后一个值是预先计算的分母）：

384.6657 44482.7215 17379.0847 1245.4225 1.4098E-26

这是U未优化的第一行：

0.99999213 2.3382613E-21 2.8218658E-17 7.900302E-6

这是U优化后的第一行：

0.9999921 2.338395E-21 2.822035E-17 7.900674E-6

ITERATION 1：

D这是未优化的矩阵的第一行：

414.3861 44469.39 17300.092 1197.7633

这是矩阵的第一行D优化方式（注意最后一个值是预先计算的分母）：

414.3880 44469.38 17300.090 1197.7657 2.0796E-26

这是U未优化的第一行：

0.99997544 4.9366603E-21 6.216704E-17 2.4565863E-5

这是U优化后的第一行：

0.3220644 1.5900239E-21 2.0023086E-17 7.912171E-6

最后一组值表明，由于传播错误（我仍然希望我犯了一些错误），它们非常不同，甚至违反了这些值的总和必须为 1.0 的约束。

难道我做错了什么？是否有可能的解决方案来优化代码和数值稳定？任何建议或批评将不胜感激。

我正在为算法课程做一些线性编程练习，在此过程中，我正在手动解决许多带有分数的运算。在这样做的过程中，我意识到人类不会遭受数字不稳定的困扰：我们只是将值保留为分数表示，然后我们最终评估（可能通过使用计算器）表达式的值。

是否有任何技术可以自动执行此操作？

我正在考虑实现某种符号计算，在内部简化数字并最终仅在表达式评估期间产生值的东西。

我有一些数学（在 C++ 中）似乎生成了一些非常小的、接近于零的数字（我怀疑 trig 函数调用是我真正的问题），但我想检测这些情况以便我可以研究它们更多详情。

我目前正在尝试以下方法，是否正确？

其次，数学的本质是三角函数（也就是使用大量的弧度/度数转换和//sin调用等），我可以做什么样的转换来避免数学错误？costan

显然，对于乘法，我可以使用对数变换- 还有什么？

在我正在使用 C 语言的数值求解器中，我需要反转一个 2x2 矩阵，然后在右侧乘以另一个矩阵：

我一直在使用倒置 2x2 矩阵的以下定义：

在我的求解器的前几次迭代中，这似乎给出了正确的答案，但是，经过几步之后，事情开始增长并最终爆炸。

现在，与使用 SciPy 的实现相比，我发现同样的数学不会爆炸。我能找到的唯一区别是 SciPy 代码使用scipy.linalg.inv()，它在内部使用 LAPACK 来执行反转。

当我用上述计算替换调用时inv()，Python 版本确实会爆炸，所以我很确定这就是问题所在。计算中的微小差异正在蔓延，这让我相信这是一个数值问题——对于反演运算来说并不完全令人惊讶。

我正在使用双精度浮点数（64 位），希望数值问题不会成为问题，但显然情况并非如此。

但是：我想在我的 C 代码中解决这个问题，而不需要调用像 LAPACK 这样的库，因为将它移植到纯 C 的全部原因是让它在目标系统上运行。此外，我想了解这个问题，而不仅仅是呼叫一个黑匣子。如果可能的话，最终我也希望它以单精度运行。

所以，我的问题是，对于这么小的矩阵，是否有一种数值上更稳定的方法来计算 A 的逆？

谢谢。

编辑：目前试图弄清楚我是否可以通过求解来避免反转C。

给定二维贝塞尔曲线（P0、P1、P2、P3）的点，我想找到给定 x 坐标的 y 坐标。由于以下限制，该问题得到了很好的定义：

• P0 = (0,0), P3 = (1,1)
• P1 = (t, 1-t) 对于 0, 1 之间的 t
• P2 = 1 - P1（x 和 y）

我有以下函数来计算答案，将上述所有限制都放入CubicBezier.html的贝塞尔曲线公式中。我正在使用 Newton-Raphson 来计算我想要的点的参数，我知道这是可行的，因为我不会让循环完成，直到它完成（在定义的公差内）。

我正在使用此功能将对比度滤镜应用于图像。因为这个 0.5 会返回相同的图像，0.0 会最大程度地降低对比度，而 1.0 会最大程度地增加。

编辑以下功能已更正，现在可以完美运行。

如果我们设置 p = 0.5，我们应该得到一条直线，但是当我对 linspace 执行此操作并绘制点时，我会得到 0.5 和 1.0 之间的弯曲。谁能明白为什么会发生这种情况，如果有什么我能做的吗？