GSL/BLAS：矩阵与逆矩阵相乘

Question

我在使用GNU GSL进行一些矩阵计算。我正在尝试将矩阵B乘以矩阵A的逆。

现在我注意到GSL的BLAS部分有一个函数可以做到这一点，但仅当A为三角形时才能使用该函数。这是有特定原因吗？此外，最快的计算方法是什么？我应该使用LU分解来求A的逆，还是有更好的方法？

顺便说一下，A的形式为P'GP，其中P是一个正常矩阵，P'是它的逆，而G是一个对角矩阵。

非常感谢 :)

- Tom

Answer 1

我认为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是一个正常的矩阵。我有一段时间没有使用它们了，但是必须有一个针对它们的分解定理，您可以在GSL参考手册第14章或更好的线性代数书中找到它们。

只是一个例子，如果您有对称正定矩阵并希望找到其逆，则可以使用Cholesky分解，这是针对该类型矩阵进行优化的。然后，您可以使用GSL中的gsl_linalg_cholesky_decomp()和gsl_linalg_cholesky_invert()函数使其高效。

希望对您有所帮助！

- Alejandro Cámara

Answer 2

简而言之：仅支持三角矩阵是因为三角矩阵的这种操作非常容易执行（称为回代），而BLAS仅提供低级函数的例程。任何更高级别的功能通常在LAPACK中找到，它内部使用BLAS进行块操作。

在您处理的特定情况下：A := P'*G*P，如果P是一个正常矩阵，则A的逆可以写成inv(A) := P'*inv(G)*P。然后您有两个选择：

1. 构建inv(A)并将其与B相乘。
2. 按顺序将B与inv(A)的不同因子相乘：先将B与P相乘（左侧），然后重新调整结果的每一行，使其与G的对角线元素的倒数成比例，然后再次与P'相乘（再次在左侧）。

根据具体情况，您可以选择自己的方法。无论如何，您都需要dgemm（矩阵-矩阵乘法，Lvl3 BLAS）和dscal（行缩放，Lvl 1 BLAS），假设是双精度。

希望这可以帮助您。

A.

- Adrien

---

谢谢您的回复。我不确定排列来自何处。您是否可能混淆了什么，因为我的矩阵之一被称为“P”？您说得对，本质上我可以将问题重新表述为(P')GP*X = B，但我不明白您如何进一步推导出您提出的内容？您能详细说明一下吗？此外，请注意，虽然我的矩阵A由P'GP组成，但这并不是某种特征值分解，在您考虑这些方面时请注意。 - Tom

我的错误，我之前正在使用置换矩阵，我会编辑帖子进行更正。 - Adrien