求解 PSD 线性系统
Solve PSD linear system
我正在尝试使用 BLAS/LAPACK 矩阵 $A - A B^{-1} A$ 进行计算，其中 A 是对称的，B 是 PSD。
最聪明的方法是什么？
我正在考虑计算 Cholesky 分解 $B = LL^T$，然后求解 $L^T X = A$。
但是最后一步怎么做呢？是否可以利用 $A$ 是对称的这一事实？
我正在使用 Cython 和提供的 BLAS/LAPACK 接口 (https://docs.scipy.org/doc/scipy/reference/linalg.cython_blas.html)。
感谢您的帮助！
假设 A 是对称的 A^T=A 并且 B 是对称正定的。
通过 LAPACK ?potrf 找到 B=LL^T 的 Cholesky 分解。
这会将L作为下三角矩阵保存到B中（注意在?potrf中设置第一个参数uplo='L'）。
作为下一步，您可以使用 LAPACK ?trtrs 为 X 求解 LX=A。
确保再次设置 uplo='L'。
使用以下计算
A B^{-1} A
= A (LL^T)^{-1} A
= A L^{-T} L^{-1} A
= (L^{-1} A)^T L^{-1} A
= X^T X
很明显只需要乘X^T X。
这可以通过 BLAS ?syrk 来完成。
以下代码计算 ABinvA := X^T X
call syrk(X, ABinvA, trans='T')
最终结果是一个简单的减法运算
res = A - ABinvA
我正在尝试使用 BLAS/LAPACK 矩阵 $A - A B^{-1} A$ 进行计算，其中 A 是对称的，B 是 PSD。 最聪明的方法是什么？ 我正在考虑计算 Cholesky 分解 $B = LL^T$，然后求解 $L^T X = A$。 但是最后一步怎么做呢？是否可以利用 $A$ 是对称的这一事实？
我正在使用 Cython 和提供的 BLAS/LAPACK 接口 (https://docs.scipy.org/doc/scipy/reference/linalg.cython_blas.html)。
感谢您的帮助！
假设 A 是对称的 A^T=A 并且 B 是对称正定的。
通过 LAPACK ?potrf 找到 B=LL^T 的 Cholesky 分解。
这会将L作为下三角矩阵保存到B中（注意在?potrf中设置第一个参数uplo='L'）。
作为下一步，您可以使用 LAPACK ?trtrs 为 X 求解 LX=A。
确保再次设置 uplo='L'。
使用以下计算
A B^{-1} A
= A (LL^T)^{-1} A
= A L^{-T} L^{-1} A
= (L^{-1} A)^T L^{-1} A
= X^T X
很明显只需要乘X^T X。
这可以通过 BLAS ?syrk 来完成。
以下代码计算 ABinvA := X^T X
call syrk(X, ABinvA, trans='T')
最终结果是一个简单的减法运算
res = A - ABinvA