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

Cholesky type stability #1117

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

Conversation

willtebbutt
Copy link
Member

@willtebbutt willtebbutt commented Nov 10, 2021

The methods for literal_getproperty on the Cholesky factorisation are type-unstable. This instability stems from the value of uplo not being known at compile time as it's dynamic.

The only way that I can see around this is to ensure that we arrive at the same type for the factors field regardless whether uplo is :L or :U.

I'm not 100% happy with using collect -- it feels like it'll be a bad option for things like CuArrays, to which this rule does in principle apply I believe. @oxinabox @mcabbott @DhairyaLGandhi @devmotion any thoughts on another way of achieving the same thing?

Note: this should not prevent us moving forward with #1114 as this doesn't interfere with the rrule for cholesky itself.

@DhairyaLGandhi
Copy link
Member

Collecting on CuArrays is a bad idea generally. Is there an mwe to show the value in fixing the instability?

test/gradcheck.jl Outdated Show resolved Hide resolved
@willtebbutt
Copy link
Member Author

Collecting on CuArrays is a bad idea generally.

Indeed -- hence my desire to find a better way to do this that doesn't involve collecting, but which does solve the instability.

Is there an mwe to show the value in fixing the instability?

In my application, I have a call to cholesky(A).U buried somewhere deep inside my code, and the instability is causing most of the reverse-pass to be unstable. I don't know how I would demonstrate the value in fixing it beyond fixing instabilities generally being a good thing though. Did you have something particular in mind?

@mcabbott
Copy link
Member

Is this branch needed at all? What's returned here goes back into the adjoint of cholesky, presumably. Does it read the wrong triangle at all? Maybe something like (uplo=nothing, info=nothing, factors=Hermitian(Δ)) would always be fine?

@devmotion
Copy link
Collaborator

What's the type of Δ in these pullbacks? Could you use triu(Δ) and tril(Δ) instead of triu!(collect(Δ)) and tril!(collect(Δ))?
I was a bit worried if tril and triu would return the same type but it seems to be the case at least for e.g. Matrix and triangular matrices:

julia> using LinearAlgebra

julia> A = rand(2, 2);

julia> typeof(triu(A)) === typeof(tril(A)) === Matrix{Float64}
true

julia> typeof(triu(UpperTriangular(A))) === typeof(tril(UpperTriangular(A))) === UpperTriangular{Float64,Matrix{Float64}}
true

julia> typeof(triu(LowerTriangular(A))) === typeof(tril(LowerTriangular(A))) === LowerTriangular{Float64,Matrix{Float64}}
true

@st--
Copy link
Contributor

st-- commented Nov 16, 2021

Related: JuliaLang/julia#42920 - question on whether we could have uplo be a type parametrization instead of a field on Cholesky?

@willtebbutt

In my application, I have a call to cholesky(A).U buried somewhere deep inside my code

are you using PDMats.chol_lower?

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

Successfully merging this pull request may close these issues.

5 participants