Implementation for accessing Cholesky factors

Hi everyone!

I’ve been a Julia user for several years, but this is my first time posting here, so apologies if this has been asked somewhere before. I was curious about Julia’s implementation of Cholesky factorization — specifically accessing and using the factors of the Cholesky data structure.

Say I have some (positive definite) matrix A = L L^T, I perform the Cholesky factorization as
C = cholesky(A), and I want to use L for some calculations (e.g., a matrix-vector multiplication). I realized recently that, given some vector v,


is quite a bit slower than


because C.L actually forms a copy of the underlying factors matrix every time it’s called (and because storage of the factor in the upper part is the default behavior). The relevant code in cholesky.jl is:

if d === :U
    return UpperTriangular(Cuplo === char_uplo(d) ? Cfactors : copy(Cfactors'))
elseif d === :L
    return LowerTriangular(Cuplo === char_uplo(d) ? Cfactors : copy(Cfactors'))
elseif d === :UL
    return (Cuplo === 'U' ? UpperTriangular(Cfactors) : LowerTriangular(Cfactors))

I could imagine that many a user might naively call C.L without realizing the significant performance hit, so I’m curious: why not return something like LowerTriangular(Cfactors') instead of LowerTriangular(copy(Cfactors')) (and similarly for C.U if the factor is stored in the lower part)? Of course, there’s a difference in meaning by returning an adjoint, and maybe something would break (I can’t think of anything at the moment) for an AbstractMatrix data structure that’s not just a simple array of floats. But for most use cases (e.g., multiplication, solves, changing coefficients, accessing coefficients, etc.), I would imagine this is advantageous. Are there some deep reasons why this choice of copying was made?

Thanks :slight_smile:!


Please submit a PR. I did a similar change for QR

Thanks, I’ll do that.