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!
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