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?

  • 3
    $\begingroup$ How large are your problems? QR is not significantly slower than LDL/Cholesky when $X$ has very imbalanced dimensions. $\endgroup$ – Federico Poloni Feb 29 '20 at 7:54
  • 1
    $\begingroup$ 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$. $\endgroup$ – Brian Borchers Feb 29 '20 at 16:44
  • $\begingroup$ @Federico Poloni $X$ in average is around 300 rows by 30 columns for some around hundred thousands of least squares. $\endgroup$ – eusoubrasileiro Feb 29 '20 at 18:02
  • $\begingroup$ @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. $\endgroup$ – eusoubrasileiro Feb 29 '20 at 18:10
  • 3
    $\begingroup$ 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. $\endgroup$ – rchilton1980 Feb 29 '20 at 18:47

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.


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


Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.