Search code examples
cythonlapackblas

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!

Solution

  • 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