CODEWITHSUNDEEP
cythonlapackblas
Sap
Sap

Reputation: 23

Solve PSD linear system

I am trying to compute with BLAS/LAPACK the matrix $A - A B^{-1} A$ where A is symmetric and B is PSD. What would be the smartest way to do this? I was thinking to compute a Cholesky decomposition $B = LL^T$ and then solve $L^T X = A$. But how to do this last step? Is it possible to exploit the fact that $A$ is symmetric?

I am using Cython and the provided BLAS/LAPACK interface (https://docs.scipy.org/doc/scipy/reference/linalg.cython_blas.html).

Thanks for your help!

Upvotes: 2

Views: 204

Answers (1)

jack
jack

Reputation: 1790

Assuming A is symmetric A^T=A and B is symmetric positive definit.

Find the Cholesky decomposition of B=LL^T by LAPACK ?potrf. This will save L into B as a lower triangular matrix (note to set the first argument uplo='L' in ?potrf).

As a next step you can solve LX=A for X by using LAPACK ?trtrs. Make sure to set uplo='L' again.

Using the following computation

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

it is clear that you only need to multiply X^T X. That can be done by BLAS ?syrk. The following code computes ABinvA := X^T X

call syrk(X, ABinvA, trans='T')

The final result is a simple subtraction operation

res = A - ABinvA

Upvotes: 1

Related Questions