# 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

## 2 Answers

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.

Computing the SVD is, of course, very expensive. In this case avoid divide and conquer SVD algorithm implementations since those, although often faster to compute, tend to be less accurate for the very small singular values.

Without special properties of matrixes you can not have a smart decompositions or think of using special methods of resolution for the linear system (for example iterative methods with preconditioning techniques to mitigate the ill-conditioning, always with in mind Wolfgang Bangerth's comment).

Considering Toepliz’s matrix you meet sometimes. From a more linked to numerical libraries perspective for what concerns lapack certainly until early 2013, it was not available a special function for Toepliz’s matrix, see lapack archives. Also now I don't found dedicate functions. About mkl I don't know special function for Toepliz’s matrix.

Perhaps it may be useful these functions of Netlib (so we remain in an environment similar to lapack/mkl) for Toepliz’s matrix.