调试数值稳定性问题的策略？
我正在尝试为Python编写Wilson的谱密度分解算法的实现。这个算法是通过反复分解一个[QxQ]的矩阵函数，找到它的平方根（可以看作是牛顿-拉夫森平方根查找方法在谱密度矩阵上的一种扩展）。
问题是，我的实现只对45x45及以下的矩阵收敛。经过20次迭代后，矩阵之间的平方和差大约是2.45e-13。但是，如果输入的矩阵是46x46，它直到第100次左右才收敛。对于47x47或更大的矩阵，根本无法收敛；误差在100到1000之间波动大约100次迭代，然后开始迅速增长。
你会怎么去调试这样的情况呢？似乎没有特定的点让它变得不正常，而且这些矩阵太大了，我根本无法手动计算。有没有人能提供一些建议、教程或经验法则，帮助找到这种奇怪的数值错误呢？
我之前从未遇到过这样的情况，但我希望你们中的一些人有经验...
谢谢，
- Dan
[1] G. T. Wilson. "矩阵谱密度的分解". SIAM应用数学杂志（第23卷，第4期，1972年12月）
2 个回答
2
区间算术可能会有所帮助，但我不确定它的性能是否足够好，能让你在你关注的矩阵大小上进行有效的调试。你得考虑到，使用这种方法可能会慢上几个数量级，因为你要把依赖硬件加速的“标量”浮点运算换成计算量大的“区间”运算，还要增加一些检查，看看哪些区间变得太宽了，什么时候、在哪里以及为什么会这样。
4
我建议你可以在scipy-user这个邮件列表上问问这个问题，最好能附上你的代码示例。一般来说，那里的人对数值计算非常有经验，而且很乐于助人，光是关注这个列表就能学到很多东西。
另外，我恐怕没有其他的想法……如果你觉得这是数值精度或者浮点数舍入的问题，你可以先试试把所有的数据类型都改成float128，看看有没有什么变化。