1
$\begingroup$

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.

|cite|improve this question
$\endgroup$
2
$\begingroup$

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.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ 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! $\endgroup$ – Federico Poloni Oct 10 at 16:46
  • $\begingroup$ 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. $\endgroup$ – rchilton1980 Oct 10 at 17:38
  • $\begingroup$ 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. $\endgroup$ – rchilton1980 Oct 10 at 17:46
  • $\begingroup$ 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? $\endgroup$ – Federico Poloni Oct 10 at 17:48
  • $\begingroup$ 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? $\endgroup$ – rchilton1980 Oct 10 at 17:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.