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

Question

7

我正在尝试编写Python版本的Wilson谱密度因式分解算法[1]。该算法将一个[QxQ]矩阵函数迭代地分解为其平方根（类似于牛顿-拉夫逊求解谱密度矩阵的扩展）。

问题是我的实现只能收敛于45x45或更小的矩阵大小。因此，在20次迭代后，矩阵之间的平方差之和约为2.45e-13。但是，如果输入大小为46x46，则需要大约100次迭代才能收敛。对于47x47或更大的矩阵，矩阵永远不会收敛；误差在约100次迭代中波动在100到1000之间，然后开始迅速增长。

您会如何尝试调试这样的问题？似乎没有任何特定的点使它变得疯狂，而且矩阵太大了，我无法手动进行计算。是否有人有关于发现此类奇怪数字错误的提示/教程/启发？

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

谢谢， - 丹

[1] G.T.威尔逊。 "Matricial Spectral Densities的分解"。 SIAM J. Appl. Math（Vol 23，No.4，Dec.1972）

- Dan

“only converges for matrixes of size 45x45” 是什么意思？那么，小于45x45的矩阵也会失败吗？ - badp

不好意思，我会编辑这篇文章。它在大小为45x45及更小的情况下成功收敛。 - Dan

2个回答

2

区间算术可以帮助，但我不确定在您感兴趣的矩阵大小上性能是否足够支持有意义的调试（您必须考虑到几个数量级的减速，因为需要将高度硬件辅助的“标量”浮点运算替换为软件密集型的“区间”运算，并添加检查哪些区间过于宽阔、何时、何地和为什么）。

- Alex Martelli

那么...你的意思是将区间矩阵插入到SciPy中吗？我甚至不确定是否可以在不重写区间数学中的linpack的情况下完成这个任务，我能吗？ - Dan

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- robince · Accepted Answer

我建议你在scipy-user邮件列表上提出这个问题，最好附上你的代码示例。通常，列表上的人似乎对数值计算非常有经验，并且非常乐于助人，仅仅关注列表本身就是一种学习。

否则，恐怕我没有任何想法...如果你认为这是一个数值精度/浮点舍入问题，你可以尝试将所有的数据类型升级到float128，看看是否有任何区别。