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

Question

我正在尝试为Python编写一个威尔逊谱密度因式分解算法1的实现。该算法迭代地将QxQ矩阵函数分解为其平方根(它是谱密度矩阵的牛顿-拉夫森平方根搜索器的一种扩展)。

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

你将如何尝试调试这样的东西呢？似乎没有任何特定的点它会变得疯狂，而且矩阵太大了，我无法实际尝试手动计算。有没有人有小贴士/教程/启发式方法来发现像这样的奇怪的数字bug？

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

谢谢，-丹

1G.T.威尔逊。“矩阵谱密度的因子分解”。SIAM J. Appl.数学(第23卷，第4期，1972年12月)

Answer 1

我建议在scipy-user邮件列表上问这个问题，也许可以用您的代码示例。一般来说，名单上的人似乎在数值计算方面经验丰富，并且真的很有帮助，仅仅遵循名单本身就是一种教育。

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

Answer 2

Interval arithmetic可以提供帮助，但我不确定性能是否足以在您感兴趣的矩阵大小下进行有意义的调试(您必须计算出几个数量级的减慢，将高度HW帮助的“标量”浮点操作替换为SW繁重的“间隔”操作，并添加检查哪些间隔变得过宽、何时、何处和原因)。