Cholesky type stability

[willtebbutt]

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]

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

[willtebbutt]

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]

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]

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--]

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?