我正在使用 GNU GSL 进行一些矩阵计算。我正在尝试将矩阵 B 与矩阵 A 的逆相乘。
现在我注意到 GSL 的 BLAS 部分具有执行此操作的功能，但前提是 A 是三角形的。这有什么具体原因吗？另外，进行此计算的最快方法是什么？我应该使用 LU 分解来反转 A，还是有更好的方法？
FWIW，A 的形式为 P'GP，其中 P 是正规矩阵，P' 是它的逆矩阵，G 是对角矩阵。
非常感谢:)
【问题讨论】：
标签： math matrix linear-algebra gsl
我相信 Adrien 是正确的，因为 BLAS 没有方阵的反函数。这取决于您用来优化其逆演算的矩阵。
一般情况下，您可以使用 LU 分解，它适用于任何方阵。即，类似于：
gsl_linalg_LU_decomp(A, p, signum);
gsl_linalg_LU_invert(A, p, invA);
其中 A 是您想要其逆的方阵，p 是 gsl_permutation 对象（对置换矩阵进行编码的置换对象），signum 是置换的符号，invA 是 A 的逆。
由于您声明A = P' G P 是P 正常和G 对角线，因此A 可能是正常矩阵。我已经有一段时间没有使用它们了，但是它们必须有一个因式分解定理，您可以在 Chapter 14 of the GSL reference manual 中找到它，或者在线性代数书中找到更好的方法。
举个例子，如果您有对称的正定矩阵并且想要找到它的逆矩阵，您可以使用 Cholesky 因式分解，它针对这种矩阵进行了优化。然后你可以在 GSL 中使用gsl_linalg_cholesky_decomp() 和gsl_linalg_cholesky_invert() 函数来提高效率。
希望对你有帮助！
【讨论】：
简而言之：
仅支持三角矩阵的事实仅仅是因为此操作对于三角矩阵非常容易执行（称为反向替换），并且 BLAS 仅提供低级函数的例程。任何更高级别的函数通常都可以在 LAPACK 中找到，它在内部使用 BLAS 进行块操作。
在您正在处理的特定情况下：A:=P'*G*P，如果P 是一个正规矩阵，那么A 的逆矩阵可以写成inv(A) := P'*inv(G)*P。然后你有两个选择：
inv(A) 并将其与B 相乘。inv(A) 的不同因子相乘：将 B 乘以 P（在左侧），然后用 G 的对角元素的倒数重新缩放结果的每一行，然后相乘再次与P'（再次离开）。 根据具体情况，您可以选择您的方法。无论如何，假设双精度，您将需要dgemm（矩阵-矩阵乘法，Lvl3 BLAS）和dscal（线的缩放，Lvl 1 BLAS）。
希望这会有所帮助。
一个。
【讨论】：