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.)

  • $\begingroup$ If $A$ is not triangular, then you could apply a Cholesky factorization to $AA'$ to get a triangular replacement for $A$. $\endgroup$ – Omnomnomnom Aug 5 '16 at 19:48
  • $\begingroup$ 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. $\endgroup$ – cfp Aug 5 '16 at 20:42
  • $\begingroup$ 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? $\endgroup$ – Omnomnomnom Aug 5 '16 at 20:58
  • $\begingroup$ 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. $\endgroup$ – cfp Aug 6 '16 at 13:30
  • 1
    $\begingroup$ 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.) $\endgroup$ – cfp Aug 7 '16 at 6:56

Your Answer

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

Browse other questions tagged or ask your own question.