use*_*579 3 linear-algebra numerical-methods julia
我广泛使用 julia's linear equation solver res = X\b。由于参数变化，我必须在我的程序中使用它数百万次。这工作正常，因为我使用的是小尺寸（最多30）。现在我想分析更大的系统，直到1000，线性求解器不再有效。
我认为可以有一个解决办法。但是我必须说，有时我的 X 矩阵是密集的，有时是稀疏的，所以我需要对这两种情况都能正常工作的东西。
该b向量是一个全为零的向量，除了一个条目，它总是1（实际上它总是最后一个条目）。此外，我不需要所有的res向量，只需要它的第一个条目。
通常当人们谈论加速线性求解器时res = X \ b，它是针对多个bs 的。但既然你b没有改变，你只是一直在改变X，因此这些技巧都不适用。
从数学的角度来看，加速这一过程的唯一方法似乎是确保 Julia 为 选择最快的求解器X \ b，即，如果您知道X是正定的，请使用 Cholesky 等。Matlab 的流程图了解它如何选择求解器用于X \ b，用于密集和稀疏X，是可用的——很可能 Julia 也实现了一些接近这些流程图的东西，但同样，也许你可以找到一些简化或快捷方式的方法。
所有与编程相关的加速（多线程——虽然每个单独的求解器可能已经是多线程的，但当每个求解器使用的线程少于内核数时，并行运行多个求解器可能是值得的；@simd如果您愿意深入研究求解器本身； OpenCL/CUDA 库；等等）然后可以应用。
如果你的问题是形式的(A - µI)x = b，这里µ是一个可变的参数，并且A，b是固定的，你可能与角化工作。
令A = PDP°whereP°表示 的倒数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。