4

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

现在我注意到 GSL 的 BLAS 部分具有执行此操作的功能，但前提是 A 是三角形的。这有什么具体原因吗？另外，进行此计算的最快方法是什么？我应该使用 LU 分解来反转 A，还是有更好的方法？

FWIW，A 具有 P' G P 的形式，其中 P 是正规矩阵，P' 是它的逆矩阵，G 是对角矩阵。

非常感谢:)

4

2 回答 2

4

我相信 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()

我希望它有帮助！

于 2010-09-05T10:44:25.710 回答
3

简而言之：

仅支持三角矩阵的事实仅仅是因为此操作对于三角矩阵非常容易执行（称为反向替换），并且 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）。

希望这可以帮助。

一个。

于 2010-08-23T16:06:02.480 回答