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

Question

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

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

您将如何尝试调试此类内容？似乎没有任何特定的点,它变得疯狂,并且矩阵太大,我实际上不能手动进行计算.有没有人有寻找奇怪的数字错误的提示/教程/启发式？

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

谢谢, - 丹

[1] GT Wilson."基质谱密度的因式分解".SIAM J. Appl.数学(Vol 23,No.4,1972年12月)

Answer 1

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

否则，恐怕我没有任何想法......如果您认为这是数值精度/浮点舍入问题，您可以尝试的第一件事就是将所有数据类型提高到最大float128，看看是否有任何区别。