Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Avoid allocation when accessing cholesky factors (.L, .U) #42920

Open
wants to merge 9 commits into
base: master
Choose a base branch
from

Conversation

st--
Copy link
Contributor

@st-- st-- commented Nov 3, 2021

In master, accessing fields L or U of Cholesky or CholeskyPivoted calls copy() on the transpose. Here I add tests thayt check that the accessors don't allocate more memory than is needed for e.g. the Adjoint() construction, and remove the unnecessary copy(). Generally julia seems to be trying hard to avoid unnecessary allocations, so I would consider this to be a bugfix, not a breaking change.

There is a related discourse thread: https://discourse.julialang.org/t/implementation-for-accessing-cholesky-factors/56219

Side note 1: there is also a copy() in the definition of _chol!(A::AbstractMatrix, ::Type{UpperTriangular}) - it's a bit less clear to me in what cases this will get run and whether it's okay to remove the copy(), so let me know if you think it'd be worth spending the time to go deeper into that also.

Side note 2: some of the code tests uplo using Cuplo === char_uplo(d), some of it using sym_uplo(Cuplo) == d, is there any particular reason for this, or would it be worth unifying it?

@st--
Copy link
Contributor Author

st-- commented Nov 3, 2021

I opened this PR as suggested by @dlfivefifty in the linked discourse thread (which I found following discussion with @KristofferC) and tagging the people who also reviewed #38389: @Keno, @dkarrasch.

@andreasnoack
Copy link
Member

I believe this makes the return type a Union and that the original motivation for the copy was to have a concrete return type.

@st--
Copy link
Contributor Author

st-- commented Nov 3, 2021

I believe this makes the return type a Union and that the original motivation for the copy was to have a concrete return type.

If you think this change should not get merged, because having a concrete return type is the most important consideration, then I can change it to add a comment that explains the reasons for explicit copy, plus making this more explicit in the docstrings (which currently give no indication at all that accessing a field would involve allocations, which I would argue is very counter-intuitive behaviour - property getters shouldn't add any overhead!).

A (higher effort, probably breaking, perhaps cleaner overall) change to get both performance (minimal allocations) and a concrete return type would be to turn uplo from a field into a type parametrization.

@andreasnoack
Copy link
Member

I think this can mostly be fixed with documentation. You can specify if you want the L or the U version when you compute the factorization and then you just need to use the same one when extracting the factor.

@st--
Copy link
Contributor Author

st-- commented Nov 3, 2021

You can specify if you want the L or the U version when you compute the factorization and then you just need to use the same one when extracting the factor.

That's assuming that you're both in control of the code constructing the factorization, and in the code using it. Which is not necessarily the case in practice. E.g. I have code that gets passed an MvNormal (from Distributions.jl) whose covariance is a PDMat (from PDMats.jl) which contains a Cholesky, and I don't actually get to say whether that Cholesky uses L- or U-representation.
It's possible of course to write my own function that checks the uplo field and returns either L or U' depending on which is more efficient - as is done in https://github.com/JuliaStats/PDMats.jl/blob/master/src/chol.jl - but I also just identified a bug in that code (JuliaStats/PDMats.jl#143). Would be much nicer if the stdlib made it easy doing the right thing!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants