# Lapack symmetric update $B^{-1}AB^{-T}$

Does Lapack have a routine that, given symmetric $$A=A^T$$ and $$B$$, computes the symmetric matrix $$B^{-1}AB^{-T}$$ (while preserving symmetry exactly)?

It would be enough to have this routine for triangular $$B$$, because of course the first step here will be computing a factorization of $$B$$ (LU or LDL, according to whether $$B$$ is symmetric, too).

Motivation: I think it would be useful in this question, especially if one factorizes the expression as suggested in the second part of my answer.

Unfortunately I don't think there's a great way to do this, at least not without some effort.

In the event that $$\mathbf A$$ is nonsingular, it might be useful to point out that the desired matrix $$\mathbf C = \mathbf B^{-1} \mathbf A \mathbf B^{-T}$$ is the inverse of the symmetric Schur complement $$\mathbf S = \mathbf B^T \mathbf A^{-1} \mathbf B$$. I find that BLAS/LAPACK has much better support for forming a matrix like $$\mathbf S$$ than one like $$\mathbf C$$, so maybe it's profitable to compute $$\mathbf S$$ first, then apply a symmetry-preserving inversion routine? The utility of this idea depends on the conditioning of $$\mathbf A$$ and $$\mathbf B$$.

If you're lucky enough that $$\mathbf A$$ is positive definite in addition to symmetric, you can use [potrf] to Cholesky factor/overwrite $$\mathbf A = \mathbf L \mathbf L^T$$. Then $$\mathbf S = \mathbf B^T \left(\mathbf L \mathbf L^T \right)^{-1} \mathbf B = \left( \mathbf L^{-1} \mathbf B \right) ^T \left( \mathbf L^{-1} \mathbf B \right)$$. Using [trsm], you can overwrite $$\tilde{\mathbf B} = \mathbf L^{-1} \mathbf B$$, then use [syrk] to compute $$\mathbf S = \tilde{\mathbf B}^T \tilde{\mathbf B}$$ into a temporary (or just overwrite $$\mathbf A$$). If $$\mathbf B$$ has full column rank, then $$\mathbf S$$ is also positive definite, and the desired output $$\mathbf C = \mathbf S^{-1}$$ can be computed using [potrf] followed by [potri]. The symmetry of $$\mathbf A$$, $$\mathbf S$$ and $$\mathbf C$$ is enforced explicitly at the API level, because all of [syrk], [potrf] and [potri] operate only upon a single triangle (with the opposing triangle assumed to match).

When $$\mathbf A$$ is invertible but indefinite, you can follow the same basic idea but computing $$\mathbf S$$ is more involved due to pivoting considerations. If you'd like to see that procedure too, just leave a comment and I will add it as an edit.

• Thanks! Even if $A$ is invertible, I would be worried about its conditioning. However, I think the idea you suggested works also without the inversion part: factor $A=LL^T$ (or $LDL^T$ if indefinite), compute $C=B^{-1}L$, then the sought quantity is $CC^T$. This might even be cheaper than the other alternatives in terms of flops. Thanks for the good idea! – Federico Poloni Oct 10 '19 at 16:46
• It's probably pretty close to a wash but only experimentation can show for certain. The answer proposes 2 [potrf]'s, 1 [trsm], 1 [syrk], 1 [potri]. The comment proposes 1 [getrf], 2 [trsm's], 1 [syrk]. Recalling that [getrf] is twice as expensive as [potrf], in the end it just boils down to 1 [potri] vs 1 [trsm]. It's worth noting that the [potri] might actually not be essential, because perhaps leaving the output $\mathbf C$ in factored form, as a solve-by-$\mathbf S$ operator, is sufficient? Might depend on context. – rchilton1980 Oct 10 '19 at 17:38
• I will re-emphasize, though, that the indefinite case is messier. Unlike the positive case [potrf], the analogous indefinite factorization [sytrf] does not render $\mathbf L$ in a form that is suitable for other BLAS routines like [trsm], instead $\mathbf L$ is all tangled up with $\mathbf D$. Furthermore, the desired output is not $\mathbf C \mathbf C^T$ anymore, it's $\mathbf C \mathbf D \mathbf C^T$. This quantity cannot be computed readily using [syrk]/BLAS. – rchilton1980 Oct 10 '19 at 17:46
• I agree with you that the indefinite case is messier. However, I am not sure about your complexity estimates for the posdef case. The method suggested in my comment needs a potrf for $LL^T$, a trsm for $B^{-1}L$ (and maybe there is still something to gain here because the RHS is triangular), and a syrk for $CC^T$. That's strictly less than the answer, am I missing anything? – Federico Poloni Oct 10 '19 at 17:48
• I've assumed there's two trsms for $\mathbf B$ / for applying $\mathbf B^{-1}$, as it is factored into LU. Unless $\mathbf B$ itself is also symmetric? – rchilton1980 Oct 10 '19 at 17:49