Slow lower triangular solves when compared to full Cholesky

I hear you. In general I don’t need a specific factor. But for some applications of Hutchinson-type stochastic trace estimation, having a symmetric factor can be enormously beneficial. I shouldn’t let the discussion in this (solved) thread drift too much, but I often need to estimate traces like \text{tr}(A^{-1} B), where A = L L^T. And as it turns out, estimating the trace with

\sum_{j=1}^M \textbf{u}_j L^{-1} B L^{-T} \textbf{u}_j

instead of

\sum_{j=1}^M \textbf{u}_j A^{-1} B \textbf{u}_j

(where the \textbf{u}_j vectors have iid Rademacher components) has substantial benefits in variance reduction. For the specific application in Gaussian processes that I work on, there is a brief discussion that the symmetrized version is at least as good here, although the theorem itself pertains to the non-symmetrized estimator.

In any case, I’ll just work with U instead and everything is fine. At this point I’m just complaining because I personally would prefer the behavior of getfield(M, :L) to give me the lazy getfield(M, :U)'. But life goes on.

1 Like