Cholesky decomposition vs LDL decomposition

Question

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?

Answer 1

As the commenters have already said, for your specific problem, QR factorization might be a better approach not because of speed but because of numerical stability. It's still a good question to ask in general.

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.

Answer 2

Not a full answer, but I will just mention two subtle details:

1. OP is applying LDLT to matrices that are positive semidefinite in exact arithmetic; hence one would expect that, barring catastrophic cancellation errors, LDLT always uses 1x1 pivots rather than 2x2 pivots. Hence the benchmark results for those matrices may differ from generic ones (based on random matrices, I presume).

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