调试数值稳定性问题的策略？

Question

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

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

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

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

谢谢， - 丹

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

Answer 1

我建议在scipy-user 邮件列表中提出这个问题，也许可以提供您的代码示例。一般来说，名单上的人似乎在数值计算方面经验丰富，并且非常乐于助人，仅仅关注名单本身就是一种教育。

否则，恐怕我没有任何想法......如果您认为这是一个数值精度/浮点舍入问题，您可以尝试的第一件事是将所有 dtypes 提高到 float128 和看看是否有什么不同。

Answer 2

Interval arithmetic 可以提供帮助，但我不确定性能是否足以真正允许在您感兴趣的矩阵大小上进行有意义的调试（您必须考虑几个数量级的减速，在替换高度硬件帮助的“标量”浮点运算与 SW 密集的“间隔”运算，并添加了关于哪些间隔增长过宽、何时、何地以及为什么增长的检查）。