# Find $X$ such that $XX'=AA'-BB'$ where $A$ and $B$ are known real matrices

Suppose that $A$ and $B$ are known real matrices (not necessarily triangular or square), and that it is known that $AA'-BB'$ is positive semi-definite (where $'$ denotes transpose).

What is the most numerically stable way to find $X$ (also not necessarily triangular or square) such that $XX'=AA'-BB'$?

If $A$ were triangular, one could perform a Cholesky down-date, but what about the case when $A$ is not triangular?

One always has the option of calculating $AA'-BB'$ and then taking a factorisation, but it seems likely that this looses precision.

(E.g. to factorize a symmetric matrix $Z$, take $Z=UDU'$, where $D$ is diagonal and $U$ is unitary, from the Schur decomposition, then if $Z$ is known to be p.s.d. we can safely zero any negative elements of the diagonal, giving $D\approx\text{diag}[d,0]$ where $d$ is a strictly-positive vector, so $Z\approx \left( U_{\cdot 1} \text{diag}\sqrt{d} \right)\left( U_{\cdot 1} \text{diag}\sqrt{d} \right)'$, where $U$ is conformably partitioned with $D$, and square-root acts element-wise on vectors.)

• If $A$ is not triangular, then you could apply a Cholesky factorization to $AA'$ to get a triangular replacement for $A$. – Omnomnomnom Aug 5 '16 at 19:48
• Doesn't help. Cholesky factoring $AA'$ throws away as much numerical precision as just factoring $AA'-BB'$. $AA'$ is already factored and we ought to be able to use that information. – cfp Aug 5 '16 at 20:42
• Are you claiming that a Cholesky down-date doesn't throw away numerical precision? How does that work if $B$ is something other than a rank $1$ matrix? – Omnomnomnom Aug 5 '16 at 20:58
• No, I'm not claiming that down dates don't throw away numerical precision. But they're likely more accurate than factorising ex post. There are Cholesky down-date algorithms even when B has rank > 1. But in my circumstance there's no Cholesky decomposition to start with, so a Cholesky down date is not what is needed. – cfp Aug 6 '16 at 13:30
• Well one can take $A'=RQ$ by the RQ decomposition, which is almost in the form you suggest. In particular, if you define $L=[R' 0]$, then $L$ is lower triangular and one has $(QX)(QX)'=LL'-(QB)(QB)'$. Alternatively, you could have defined $A'=QR$ by the QR decomposition, so that with $L=[R' 0]$, $AA'=R'R=LL'$. (Not clear which approach would be more accurate.) But with either approach, one still can't apply the standard Cholesky downdate, since $XX'$ is only known to be p.s.d., not positive-definite. (Indeed $XX'$ must have some zero eigenvalues when $A$ is not square.) – cfp Aug 7 '16 at 6:56