3
$\begingroup$

Background

I am decomposing (and then operating on) a symmetric positive definite matrix (a covariance matrix) as part of a larger system. Dimensions range from O(10) to O(300). The covariance matrices have reasonably dimensioned diagonal elements (i.e. do not differ by too many orders of magnitude), but are very strongly correlated and hence have poor condition number, sometimes as large as $10^{13}$ and thus is close to positive semi-definite.

I've noted that Eigen offers both standard Cholesky $LL^T$ decomposition and pivoted $LDL^T$ modified Cholesky decomposition. Notably, it does not seem to offer un-pivoted modified Cholesky decomposition (e.g. Algorithm 4.1.1 in Golub).

Golub claims that "if A is symmetric and positive definite, then Algorithm 4.1.1 not only runs to completion, but is extremely stable." Later, in s.4.2.6, Golub discusses the stability of the standard Cholesky decomposition $LL^T$, calling "the eigenvalue $\lambda_{min}$ is the 'distance to trouble' in the Cholesky setting." Golub then discusses pivoting strategies for $LDL^T$ without discussing the stability of the un-pivoted $LDL^T$ decomposition itself.

For reasons of fitting into the larger problem, I would like to avoid going down the pivoting path.

Questions

  1. What, if any, numerical advantages are there to the un-pivoted $LDL^T$ compared to standard Cholesky $LL^T$?
  2. Are there any criteria or "rules of thumb" when choosing one over the other?
  3. Much of the higher-level algorithms involve rank-1 updates of the $LL^T$ or $LDL^T$ decomposition. Is one variant considered to be substantially more stable than the other in this scenario?
  4. At what point does a pivoting strategy become unavoidable?

References

Golub, G. and van Loan, C., "Matrix Computations", 4th Ed. 2013.

$\endgroup$
3
  • $\begingroup$ You'd like to avoid pivoting because finding the pivot is expensive? $\endgroup$ Jul 17, 2023 at 22:46
  • $\begingroup$ The $LDL^T$ version can avoid square roots of numbers that are negative due to numerical precision. Note that your $LL^T$ can be written as $(LD^{1/2})(LD^{1/2})^T$. For very large matrices it could happen that you get a negative sign due to numerical error where it should have been positive. Taking a square root of that will ruin things. $LDL^T$ circumvents that. I can't say much about the pivoted variants though because I haven't implemented those. $\endgroup$
    – lightxbulb
    Jul 18, 2023 at 11:28
  • 1
    $\begingroup$ Just read a paper ... Stewart [1] claims in context of Kalman Filtering that LDL^T is bad news since the conditioning of D is the same as the original covariance (with cond(L) = 1), whereas for LLT, the conditioning of L will have half the exponent. "Therefore using the LDLT Kalman does not possess the wordlength advantages of Cholesky, and furthermore will suffer from overflow/underllow." Could this also be said for other problems? [1] Stewart & Chapman "Fast stable Kalman filtering utilizing the square root", IEEE ASSP 1990. ieeexplore.ieee.org/document/115844 $\endgroup$
    – Damien
    Nov 6, 2023 at 20:46

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.