Improve chol_lower and chol_upper

[devmotion]

This PR improves chol_lower and chol_upper by using the field factors directly:

julia> using PDMats, LinearAlgebra, BenchmarkTools

julia> C = cholesky(Symmetric(Matrix(I, 250, 250), :L));

julia> @btime PDMats.chol_lower($C);
  # master: 11.995 ns (1 allocation: 16 bytes)
  # this PR: 8.097 ns (1 allocation: 16 bytes)

julia> @btime PDMats.chol_upper($C);
  # master: 15.072 ns (2 allocations: 32 bytes)
  # this PR: 8.106 ns (1 allocation: 16 bytes)

julia> C = cholesky(Symmetric(Matrix(I, 250, 250), :U));

julia> @btime PDMats.chol_lower($C);
  # master: 15.248 ns (2 allocations: 32 bytes)
  # this PR: 8.116 ns (1 allocation: 16 bytes)

julia> @btime PDMats.chol_upper($C);
  # master: 12.136 ns (1 allocation: 16 bytes)
  # this PR: 8.100 ns (1 allocation: 16 bytes)

[andreasnoack]

Is the type instability a concern?

[devmotion]

If you mean if it caused me problems in any downstream application, then no, it didn't. I noticed it when I played around with FluxML/Zygote.jl#1117 and this PR seemed like a simple improvement. But I'm not too invested in the PR and therefore didn't push it, I just remembered it when I went through my list of open PRs yesterday.

[andreasnoack]

I was considering the type instability in the new version a.uplo === 'U' ? UpperTriangular(a.factors) : UpperTriangular(a.factors'). The alternative would be something like a.uplo === 'U' ? UpperTriangular(a.factors) : UpperTriangular(copy(a.factors')). It would make it type stable at the cost of the copy.

[devmotion]

Ah, sorry, I misunderstood. The PR doesn't change this type instability, the return types are exactly the same as on the master branch. For completeness, I benchmarked the version on the master branch, the version in this PR, and the version in this PR + copy:

julia> using BenchmarkTools, LinearAlgebra

julia> chol_lower(a::Cholesky) = a.uplo === 'L' ? a.L : a.U'
chol_lower (generic function with 1 method)

julia> chol_upper(a::Cholesky) = a.uplo === 'U' ? a.U : a.L'
chol_upper (generic function with 1 method)

julia> function chol_lower_pr(a::Cholesky)
           return a.uplo === 'L' ? LowerTriangular(a.factors) : LowerTriangular(copy(a.factors'))
       end
chol_lower_pr (generic function with 1 method)

julia> function chol_lower_pr(a::Cholesky)
           return a.uplo === 'L' ? LowerTriangular(a.factors) : LowerTriangular(a.factors')
       end
chol_lower_pr (generic function with 1 method)

julia> function chol_lower_pr_copy(a::Cholesky)
           return a.uplo === 'L' ? LowerTriangular(a.factors) : LowerTriangular(copy(a.factors'))
       end
chol_lower_pr_copy (generic function with 1 method)

julia> function chol_upper_pr_copy(a::Cholesky)
           return a.uplo === 'U' ? UpperTriangular(a.factors) : UpperTriangular(copy(a.factors'))
       end
chol_upper_pr_copy (generic function with 1 method)

julia> C = cholesky(Symmetric(Matrix(I, 250, 250), :L));

julia> typeof(chol_lower(C))
LowerTriangular{Float64, Matrix{Float64}}

julia> typeof(chol_lower_pr(C))
LowerTriangular{Float64, Matrix{Float64}}

julia> typeof(chol_lower_pr_copy(C))
LowerTriangular{Float64, Matrix{Float64}}

julia> @btime chol_lower($C);
  12.179 ns (1 allocation: 16 bytes)

julia> @btime chol_lower_pr($C);
  8.839 ns (1 allocation: 16 bytes)

julia> @btime chol_lower_pr_copy($C);
  4.598 ns (0 allocations: 0 bytes)

julia> typeof(chol_upper(C))
UpperTriangular{Float64, Adjoint{Float64, Matrix{Float64}}}

julia> typeof(chol_upper_pr(C))
UpperTriangular{Float64, Adjoint{Float64, Matrix{Float64}}}

julia> typeof(chol_upper_pr_copy(C))
UpperTriangular{Float64, Matrix{Float64}}

julia> @btime chol_upper($C);
  16.554 ns (2 allocations: 32 bytes)

julia> @btime chol_upper_pr($C);
  8.933 ns (1 allocation: 16 bytes)

julia> @btime chol_upper_pr_copy($C);
  69.869 μs (2 allocations: 488.36 KiB)

julia> C = cholesky(Symmetric(Matrix(I, 250, 250), :U));

julia> typeof(chol_lower(C))
LowerTriangular{Float64, Adjoint{Float64, Matrix{Float64}}}

julia> typeof(chol_lower_pr(C))
LowerTriangular{Float64, Adjoint{Float64, Matrix{Float64}}}

julia> typeof(chol_lower_pr_copy(C))
LowerTriangular{Float64, Matrix{Float64}}

julia> @btime chol_lower($C);
  16.602 ns (2 allocations: 32 bytes)

julia> @btime chol_lower_pr($C);
  8.655 ns (1 allocation: 16 bytes)

julia> @btime chol_lower_pr_copy($C);
  69.314 μs (2 allocations: 488.36 KiB)

julia> typeof(chol_upper(C))
UpperTriangular{Float64, Matrix{Float64}}

julia> typeof(chol_upper_pr(C))
UpperTriangular{Float64, Matrix{Float64}}

julia> typeof(chol_upper_pr_copy(C))
UpperTriangular{Float64, Matrix{Float64}}

julia> @btime chol_upper($C);
  12.153 ns (1 allocation: 16 bytes)

julia> @btime chol_upper_pr($C);
  8.617 ns (1 allocation: 16 bytes)

julia> @btime chol_upper_pr_copy($C);
  4.588 ns (0 allocations: 0 bytes)

As expected the version with copy fixes the type instability, and thus gets rid of the remaining allocation if the factor is not copied. This leads also to better performance (for the example above 4.5ns vs 9ns with this PR). However, if the factors have to be copied performances is much much worse: in the example here it increases from around 9ns with this PR to around 70μs (!).

So I don't think it's worth fixing the type instability.

[codecov-commenter]

Codecov Report

Merging #149 (e9bea7f) into master (8727724) will increase coverage by 0.47%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #149      +/-   ##
==========================================
+ Coverage   89.67%   90.14%   +0.47%     
==========================================
  Files           8        8              
  Lines         426      406      -20     
==========================================
- Hits          382      366      -16     
+ Misses         44       40       -4     

Impacted Files	Coverage Δ
src/chol.jl	100.00% <100.00%> (ø)
src/pdsparsemat.jl	94.11% <0.00%> (-0.23%)	⬇️
src/pdiagmat.jl	98.95% <0.00%> (-0.07%)	⬇️
src/utils.jl	91.66% <0.00%> (+3.29%)	⬆️

---

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 8727724...e9bea7f. Read the comment docs.