# Cholesky decomposition vs LDL decomposition

In different books and on Wikipedia, you can see mentions of Cholesky decomposition and only sometimes of LDL decomposition.

As far as I understand, LDL decomposition can be applied to a broader range of matrices (we don't need a matrix to be positive-definite).

I am specially looking to solve hundreds of thousands last squares problems like bellow:

$$\mathbf{X} \cdot \mathbf{\beta} = z$$ $$\mathbf{\hat{\beta}} = \left(\mathbf{X}^T \mathbf{X}\right)^{-1} \mathbf{X}^T z$$

Where sometimes $$\mathbf{X}^T \mathbf{X}$$ due colinearity is not very well conditioned. I choosed to not use QR factorization or SVD due: 1. being much slower than cholesky, 2. I can afford discard some $$\beta$$s if performance gain is worth it.

Forgot to mention but I am also considering those algorithms on a GPU performance optimization perspective. Considering, when possible (SVD is not good on GPU), I can use a batched version to launch many thousands problems at once.

Is not LDL in general better than pure Cholesky? Is it much slower?

• How large are your problems? QR is not significantly slower than LDL/Cholesky when $X$ has very imbalanced dimensions. – Federico Poloni Feb 29 at 7:54
• You should be concerned about why the $X^{T}X$ matrix is sometimes ill-conditioned and what that means for the corresponding estimate of $\beta$. – Brian Borchers Feb 29 at 16:44
• @Federico Poloni $X$ in average is around 300 rows by 30 columns for some around hundred thousands of least squares. – eusoubrasileiro Feb 29 at 18:02
• @Brian Borchers $X^TX$ is ill-conditioned sometimes due colienearity on $X$ columns. $X$ comes from a statistical hypothesis test of unit root (augmented dickey fuller test). That happens when for some specific data-sets there is "similarities between lagged dickey fuller tests". I will get the link of another question here soon. – eusoubrasileiro Feb 29 at 18:10
• I agree with Federico, QR factorization seems to be the right tool for this problem. In general it's a slower algorithm than cholesky or LDL on square inputs, but in least squares (tall/skinny QR) you don't need $\mathbf X^T \mathbf X$, so you avoid the cost to form it (using syrk or gemm) and the squaring of condition number. – rchilton1980 Feb 29 at 18:47

One of the advantages you cite is that $$LDL^*$$ can be used for indefinite matrices, which is definitely a point in its favor. The linear algebra library Eigen, which I highly recommend, has some benchmarks about this which seem to show that $$LL^*$$ is much faster for large matrices (> 1000 x 1000). This table on their website is also interesting. It seems to say that they haven't implemented blocking optimizations to improve cache utilization for $$LDL^*$$, which might account for why it's slower on large matrices at present. On the other hand, that same table also says that the accuracy of the $$LL^*$$ form degrades more with the condition number than the $$LDL^*$$ form, which should be true whether or not they've implemented cache blocking. Which form is better then becomes a situational question depending on the size and condition number of your system.
Naively, I would expect that an optimized implementation of $$LDL^*$$ would be faster because it uses fewer square root and division operations, which require something like 4x more cycles than floating point add or multiply on modern CPUs. (This 4x factor might be different on GPUs.) It could be that I was wrong because of something fundamental I'm missing or because one of the implementations has been more optimized than the other. In either case, it just illustrates the importance of making your own measurements.
2. Moreover, if LDLT were to take a 2x2 pivot even once on one of those matrices, then it means that the computational errors have perturbed $$X^T X$$ by more than its distance to singularity. So, arguably, in this case the results given by LDLT and the computed value of $$\beta$$ are likely to be garbage. (The results given by Cholesky on the same example are likely to be garbage, too, so this one is not a point in favor of Cholesky. It just means that the 2x2 code path in LDLT does not really matter here.)