# numerically stable routines to compute $M = B A^{-1} B$

Rather than

1. gesv -> solve $AX = B$
2. gemm -> compute $M = BX$,

somehow I feel there are better ways to compute $M$ with lapack/mkl?

• Do $A$ and $B$ have some special properties/structure, or are just generic dense, non-symmetric, square matrices? – Stefano M Nov 7 '16 at 21:18
• In particular, is $A$ symmetric? Because in that case, you can factorize $A=L^TL$. If furthermore $B$ is also symmetric, then you have $M=(L^{-1}B)^T(L^{-1}B)$. – Wolfgang Bangerth Nov 7 '16 at 21:28
• @StefanoM dense in general, sometimes toepliz like,,, – lorniper Nov 7 '16 at 21:29
• @WolfgangBangerth In general not, I am just concerned about ill-conditionning in general case if we do gesv then gemm... – lorniper Nov 7 '16 at 21:31
• @lorniper -- if your matrix is ill-conditioned, then there is nothing you can do about it. But I would be surprised if you have a matrix that is ill-conditioned (but not singular) and still small enough that you can consider forming $M$ explicitly. – Wolfgang Bangerth Nov 7 '16 at 22:41

If numerical stability and robustness is your only concern, you might try computing the SVD of $A$, and using that to compute the product. Specifically, let the SVD of $A$ be: $$A = U \Sigma V^T.$$ Then $$B A^{-1}B = (BV)\Sigma^{-1}(U^TB).$$ This may even work in cases where $A$ is ill-conditioned, but $B A^{-1} B$ is not, due to the spaces on which $B$ and $B^T$ span, and their interaction with the dominant subspaces of $A$. And if it fails, you can look and see how exactly it is failing based on which modes that are being amplified, then decide what to do from there.