Avoid copy in `chol_lower`

[st--]

Fixes #143

[st--]

Side note: chol_lower(::CholTypeSparse) seems to be working fine, but I'm not sure I understand enough about the implementation details of SuiteSparse to guarantee this.

[st--]

I agree that CholType is redundant but please make that change in a separate PR.

Done, see #145

[st--]

I think you can just make this chol_lower(a::Matrix) = cholesky(Hermitian(a, :L)).L.

Writing it in terms of chol_lower(cholesky(a)) seems cleaner to me, and it should compile to the same code; what benefits do you see in writing it as this ?

[andreasnoack]

and it should compile to the same code

No. The default is upper

julia> cholesky(A).uplo
'U': ASCII/Unicode U+0055 (category Lu: Letter, uppercase)

julia> cholesky(Hermitian(A, :L)).uplo
'L': ASCII/Unicode U+004C (category Lu: Letter, uppercase)

[st--]

Ok, not literally the same code, but chol_lower(cholesky(A)) which is equivalent to cholesky(A).U' is no less efficient (more efficient, actually) than introducing the Hermitian:

julia> using LinearAlgebra, PDMats

julia> A=rand(100,100); C=A'A;

julia> @allocated cholesky(Hermitian(C, :L)).L;  # compile

julia> @allocated cholesky(Hermitian(C, :L)).L
80208

julia> @allocated cholesky(C).U';  # compile

julia> @allocated cholesky(C).U'
80192

julia> @allocated PDMats.chol_lower(cholesky(C));  # compile

julia> @allocated PDMats.chol_lower(cholesky(C))
80192

[andreasnoack]

It's just an artifact of the not wrapping it in a function

julia> chol_lower_an(A::Matrix) = cholesky(Hermitian(A, :L)).L;

julia> @allocated chol_lower_an(A);

julia> @allocated chol_lower_an(A)
80144

It's also best to avoid wrapper types when possible.

[st--]

It's just an artifact of the not wrapping it in a function

Ah, thanks for pointing that out.

It's also best to avoid wrapper types when possible.

Could you expand on what you mean ?

[andreasnoack]

Could you expand on what you mean ?

Sure. The more types layers you have the less likely it is that a specific method has been written for it and you risk hitting a less optimized fallback. The compiler can also sometimes completely remove these wrapper types but I believe that the more layers you have, the harder it is for the compiler to do so.

[st--]

Ok, so do I understand correctly, you would prefer

-chol_lower(a::Matrix) = chol_lower(cholesky(a))
+chol_lower(a::Matrix) = cholesky(Hermitian(a, :L)).L

because the former returns a LowerTriangular{Adjoint{Matrix}}, whereas the latter returns a LowerTriangular{Matrix} directly?

[andreasnoack]

That is correct

[st--]

Thanks, that makes sense. I've changed it accordingly:)

[codecov-commenter]

Codecov Report

Merging #144 (1ce765d) into master (a0e7191) will not change coverage. The diff coverage is 100.00%.

@@           Coverage Diff           @@
##           master     #144   +/-   ##
=======================================
  Coverage   88.62%   88.62%           
=======================================
  Files           8        8           
  Lines         378      378           
=======================================
  Hits          335      335           
  Misses         43       43           

Impacted Files	Coverage Δ
src/chol.jl	100.00% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update a0e7191...1ce765d. Read the comment docs.

[st--]

Downside of merging it in with Hermitian() is that this breaks AD :S

[andreasnoack]

Oh no. Do you have an example?

[st--]

Yes, it turns out to be a Zygote issue, I've made a PR to fix it here: https://github.com/FluxML/Zygote.jl/pull/1114

[st--]

@andreasnoack looks like the Zygote issue might be a bit more complex to resolve; how about (at least temporarily) changing it from Hermitian to Symmetric?

[andreasnoack]

Given that we have this restriction https://github.com/JuliaStats/PDMats.jl/blob/635a957abc71aa9b63d276f96a185db1bed9de74/src/PDMats.jl#L36 it should be fine