在 Julia 中有效地求解一个特定的线性系统答案

【问题标题】：Solve a particular linear system efficiently in julia在 Julia 中有效地求解一个特定的线性系统
【发布时间】：2016-11-24 00:16:03
【问题描述】：

我广泛使用 julia 的线性方程求解器res = X\b。由于参数变化，我必须在我的程序中使用它数百万次。这工作正常，因为我使用的是小尺寸（高达30）。现在我想分析更大的系统，直到1000，线性求解器不再有效。

我认为可以解决这个问题。但是我必须说，有时我的 X 矩阵是密集的，有时是稀疏的，所以我需要在这两种情况下都能正常工作的东西。

b 向量是一个全为零的向量，除了一个始终为 1 的条目（实际上它始终是最后一个条目）。此外，我不需要所有的res 向量，只需要它的第一个条目。

【问题讨论】：

这可能不是您想要的，但我已经非常迷恋 GPU 的强大功能。如果您可以访问该硬件，它可能对您的情况非常有用。
你必须解决它数百万次。保存 LU 分解呢？如果您不更改 X 每次迭代，这将起作用。只需执行X=lufact(X) 然后X\b 应该会更快。但这是如果你改变 b （这发生在很多 PDE 求解器中）。
也许它会起作用，但X 确实会有所不同。它以X = x*eye(N,N) - T 的形式编写，其中T 是一个固定矩阵。 x是我必须解决系统res = X\b的参数。
减去与单位矩阵成比例的东西会以难以预测的方式改变分解，所以我认为分解在这里没有帮助。 You may want to ask here 如果有人有数学技巧可以更快地解决这个问题（使用 x*eye(N,N) - T 的因式分解来获得不同 x 的因式分解）。有人可能有一篇研究文章，因为这可能与进行相关的特征值计算有多么相关，但我认为做好这件事在数学上比预期的要困难。
如果 GPU 是一个选项，您可能想尝试cublas<t>getrfBatched() 和cublas<t>getrsBatched()，它们针对一次解决一批小型线性系统进行了优化。

标签： julia linear-algebra numerical-methods

【解决方案1】：

提高效率的最佳方法是使用：JuliaMath/IterativeSolvers.jl。对于A * x = b 的问题，我会推荐x = lsmr(A, b)。

第二个最佳选择是向编译器提供更多信息：如果 Cholesky 分解适合您，则不要使用 x = inv(A'A) * A' * b，而是使用 x = inv(cholfact(A'A)) A' * b。否则，你可以试试U, S, Vt = svd(A) 和x = Vt' * diagm(sqrt.(S)) * U' * b。

不确定x = pinv(A) * b 是否经过优化，但可能比x = A \ b 更有效。

【讨论】：

【解决方案2】：

如果您的问题是(A - µI)x = b 的形式，其中µ 是一个可变参数，而A、b 是固定的，您可以使用对角化。

设A = PDP° 其中P° 表示P 的倒数。那么(PDP° - µI)x = b可以转化为

(D - µI)P°x = P°b, 
P°x = P°b / (D - µI), 
x = P(P°b / (D - µI)).

（/ 操作表示各个向量元素除以标量Dr - µ。）

在对 A 进行对角化后，计算任何 µ 的解会减少到两个矩阵/向量乘积，或者如果您还可以预先计算 P°b，则减少为单个乘积。

数值不稳定性将出现在A的特征值附近。

【讨论】：

【解决方案3】：

通常当人们谈论加速线性求解器res = X \ b 时，它是针对多个bs。但是由于你的b 没有改变，而你只是不断地改变X，这些技巧都不适用。

从数学角度来看，加快这一速度的唯一方法似乎是确保 Julia 为 X \ b 选择最快的求解器，即，如果您知道 X 是正定的，请使用 Cholesky 等.Matlab’s flowcharts 了解它如何选择用于X \ b 的求解器，用于密集和稀疏X，可用——Julia 很可能也实现了与这些流程图相似的东西，但同样，也许你可以找到一些简化的方法或快捷方式。

所有与编程相关的加速（多线程 - 虽然每个单独的求解器可能已经是多线程的，但当每个求解器使用的线程数少于内核时，可能值得并行运行多个求解器；@simd 如果您愿意深入了解求解器本身；OpenCL/CUDA 库等）然后可以应用。

【讨论】：

不出所料，我错关了数学加速的大门！见@Yves Daoust’s answer。