如何高效求解1e6阶Hilbert矩阵大规模线性方程组

你现在用的实现思路从根上就不适合这个场景，跑不出结果是必然的：

• numpy.linalg.solve()是给通用稠密矩阵设计的高斯消元实现，时间复杂度O(n³)，空间复杂度O(n²)。n=1e6规模下，光存储双精度格式的完整稠密Hilbert矩阵就需要8TB内存，你碰到的内存占满、十几小时无结果，本质是系统在疯狂做磁盘交换页操作，别说18小时，跑几十天也出不了有效结果。之前n=1e5规模能16分钟跑完，纯粹是因为该规模下稠密矩阵占80GB内存，你的设备内存刚好能扛住大部分负载，交换页开销没到完全卡死的程度。
• Hilbert矩阵是出了名的极度病态结构化矩阵，n=12左右用双精度求解就会出现明显数值误差，n>30之后双精度算出来的结果已经完全没有参考价值，n=1e6时矩阵条件数是远超浮点数表示上限的天文数字，就算你有无限内存、算力拉满，用标准双精度算出来的结果和理论值的偏差会大到根本没法用。

• 如果是完成课业答题：题目已经明确给出理论解是长度100万的全1向量，直接跑x = np.ones(10**6)就是标准答案，根本不需要跑数值求解流程。
• 如果要做数值求解流程验证，立刻抛弃稠密矩阵存储+通用高斯消元的思路：
  • 用好Hilbert矩阵的特殊结构：Hilbert矩阵属于Hankel矩阵，元素值只和行列索引之和i+j有关，这类结构化矩阵有专门的快速求解算法，时间复杂度能压到O(n log²n)，而且不需要存储完整矩阵，只需要存第一行、最后一列总共2n-1个元素就行，n=1e6规模下内存占用才十几MB，完全不会爆内存。
  • 别用双精度浮点数计算，换用足够字长的多精度算术库，不然阶数过100之后计算结果就没有数值意义。要提前说明：n=1e6规模下就算用快速算法，多精度算术带来的计算开销依然极高，没有实际工程计算价值，只适合小阶数场景做原理验证。
• 小阶数场景下还可以直接用Hilbert矩阵逆的闭式解析公式计算逆元素，再和右端向量相乘得到解，效率比通用消元法高几个量级。

内容的提问来源于stack exchange，提问作者ortunoa