JuliaStats / PDMats.jl

Uniform Interface for positive definite matrices of various structures
Other
104 stars 43 forks source link

Implement fallbacks for (un)whiten(!) and (inv)quad(!) to AbstractArray #162

Closed oxinabox closed 2 years ago

oxinabox commented 2 years ago

This is the partner to https://github.com/JuliaStats/Distributions.jl/pull/1552#issuecomment-1130483202 which will allow resolving https://github.com/JuliaStats/Distributions.jl/issues/1219

It doesn't go and implement everything falling back to AbstractArray, just the operations required for Multivariate Normals, and their family (so I can test them)

In some cases that implementation is just making the PDMat one generic, (since cholesky(x::PDMat)=x.chol is defined, etc.

codecov-commenter commented 2 years ago

Codecov Report

Merging #162 (334ca80) into master (92a5178) will decrease coverage by 0.14%. The diff coverage is 92.00%.

@@            Coverage Diff             @@
##           master     #162      +/-   ##
==========================================
- Coverage   90.33%   90.19%   -0.15%     
==========================================
  Files           8        8              
  Lines         414      418       +4     
==========================================
+ Hits          374      377       +3     
- Misses         40       41       +1     
Impacted Files Coverage Δ
src/pdmat.jl 97.87% <ø> (-0.38%) :arrow_down:
src/utils.jl 91.66% <80.00%> (ø)
src/generics.jl 87.50% <95.00%> (+2.88%) :arrow_up:

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 92a5178...334ca80. Read the comment docs.

oxinabox commented 2 years ago

Looks good to me in principle, but I would like to get someone else's opinion (@andreasnoack maybe?).

When commenting on the function signatures in the docstrings, I noticed that the general definitions for AbstractMatrix could cause surprising slow-downs for AbstractPDMats with non-StridedArrays if their specializations for (inv)quad! etc. are only defined for StridedArrays. So maybe it would be good to check if such restrictions still exist and if we can remove them.

AFAIK those don't slow down. They were MethodErrors before and now they are AmbiguityErrors

devmotion commented 2 years ago

Are you sure? To me it seems e.g. that https://github.com/JuliaStats/PDMats.jl/blob/92a51785c69a303a6da5142433b7556e1d68d66d/src/pdiagmat.jl#L113 defines an optimization of quad! for a::PDiagMat and x::StridedMatrix but by adding the generic fallback in this PR a call of e.g. quad!(r::AbstractArray, a::PDiagMat, x::SubArray{Float64,2}) will not use this specialized implementation but the less efficient fallback.

oxinabox commented 2 years ago

Looks fine to me:

With this PR:

julia> @which invquad!([2,2], PDiagMat([1,1]), [3 3; 3 3])
invquad!(r::AbstractArray, a::PDiagMat, x::StridedMatrix{T} where T) in PDMats at /home/oxinabox/JuliaEnvs/Distributions.jl/dev/PDMats/src/pdiagmat.jl:127

julia> @which invquad!([2,2], PDiagMat([1,1]), @view([3 3; 3 3][:,:]))
invquad!(r::AbstractArray, a::PDiagMat, x::StridedMatrix{T} where T) in PDMats at /home/oxinabox/JuliaEnvs/Distributions.jl/dev/PDMats/src/pdiagmat.jl:127

julia> @which invquad!([2,2], PDiagMat([1,1]), @view([3 3; 3 3][[2 4; 1 3]]))
invquad!(r::AbstractArray, a::AbstractMatrix, x::AbstractMatrix) in PDMats at /home/oxinabox/JuliaEnvs/Distributions.jl/dev/PDMats/src/generics.jl:140

Without this PR:

julia> @which invquad!([2,2], PDiagMat([1,1]), [3 3; 3 3])
invquad!(r::AbstractArray, a::PDiagMat, x::StridedMatrix{T} where T) in PDMats at /home/oxinabox/.julia/packages/PDMats/pIpaj/src/pdiagmat.jl:127

julia> @which invquad!([2,2], PDiagMat([1,1]), @view([3 3; 3 3][:,:]))
invquad!(r::AbstractArray, a::PDiagMat, x::StridedMatrix{T} where T) in PDMats at /home/oxinabox/.julia/packages/PDMats/pIpaj/src/pdiagmat.jl:127

julia> @which invquad!([2,2], PDiagMat([1,1]), @view([3 3; 3 3][[2 4; 1 3]]))
ERROR: no unique matching method found for the specified argument types

The only difference I see for this case is that the last one is now longer an error

devmotion commented 2 years ago

The only difference I see for this case is that the last one is now longer an error

Yes, exactly, that's what I meant - it's not an error anymore but it uses the generic fallback instead of the method for diagonal matrices.

oxinabox commented 2 years ago

ah right, there is a lot of code here that I think only works for StridedVectors, or rather code that only really works for things indexed by 1:length(x). I think easiest would be to go through and generalize every single method in the package. Which I am happy to make a follow up PR to do, but i think it would make this PR to large. I will start that PR now, branching off this one.

oxinabox commented 2 years ago

Given #163 is open, how do you feel about merging this? The only behavourally change is to things that errored before. And we have PRs open that will also make some of those fast.

matbesancon commented 2 years ago

Any idea where the nightly tests error comes from? https://github.com/JuliaStats/PDMats.jl/runs/6593682307?check_suite_focus=true#step:6:185

BandedMatrices: Error During Test at /home/runner/work/PDMats.jl/PDMats.jl/test/specialarrays.jl:50
  Got exception outside of a @test
  UndefVarError: liblapack not defined
  Stacktrace:
    [1] pbtrf!(uplo::Char, m::Int64, kd::Int64, A::Matrix{Float64})
      @ BandedMatrices ~/.julia/packages/BandedMatrices/dt29n/src/lapack.jl:210
    [2] banded_chol!
      @ ~/.julia/packages/BandedMatrices/dt29n/src/symbanded/BandedCholesky.jl:4 [inlined]
devmotion commented 2 years ago

Seems to be caused by the tests with BandedMatrices, so it might not be related to (or rather caused by) PDMats.

devmotion commented 2 years ago

The same test errors we've seen for a while on the master branch (eg https://github.com/JuliaStats/PDMats.jl/runs/6388016866?check_suite_focus=true), so it does not seem related to this PR.

oxinabox commented 2 years ago

can this please be merged? this is blocking some work i want to wrap up

matbesancon commented 2 years ago

done, thanks!