search

JuliaStats / PDMats.jl

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

Fix PDiagMat when inv(v) isn't a subtype of v #141

Closed Wynand closed 2 years ago

Wynand commented 2 years ago

Currently for any AbstractVector v where typeof(inv(v)) isn't a subtype of typeof(v)) creating a PDiagMat will fail

For example, if we use a range such a 0.1:0.1:0.5 to construct a PDiagMat we get this error:

julia> using PDMats
[ Info: Precompiling PDMats [90014a1f-27ba-587c-ab20-58faa44d9150]

julia> PDiagMat(0.1:0.1:0.5)
ERROR: MethodError: Cannot `convert` an object of type Vector{Float64} to an object of type StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}
Closest candidates are:
  convert(::Type{T}, ::AbstractRange) where T<:AbstractRange at range.jl:148
  convert(::Type{T}, ::LinearAlgebra.Factorization) where T<:AbstractArray at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/factorization.jl:58
  convert(::Type{T}, ::T) where T<:AbstractArray at abstractarray.jl:14
  ...
Stacktrace:
 [1] PDiagMat{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}(d::Int64, v::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, inv_v::Vector{Float64})
   @ PDMats ~/git-repos/open_source/julia_packages/PDMats.jl/src/pdiagmat.jl:9
 [2] PDiagMat(v::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, inv_v::Vector{Float64})
   @ PDMats ~/git-repos/open_source/julia_packages/PDMats.jl/src/pdiagmat.jl:15
 [3] PDiagMat(v::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}})
   @ PDMats ~/git-repos/open_source/julia_packages/PDMats.jl/src/pdiagmat.jl:18
 [4] top-level scope
   @ REPL[2]:1

julia>

To fix this I've changed this function to create a PDiagMat with a Union of the two AbstractVector types as the types for the two vectors

devmotion commented 2 years ago

An alternative would be to remove the inv_diag field from PDDiagMat (and similarly inv_value from ScalMat). From a quick look it seems that it should not be needed if multiplications with inv_diag are replaced with divisions by diag. This would also avoid the issue in this PR and unwanted promotions: Eg, the fix suggested by @andreasnoack would cause integer-valued diagonal matrices to be promoted to floating point numbers which would affect eg multiplication which should not depend on the type of the inverse (IMO).

devmotion commented 2 years ago

Fixed by https://github.com/JuliaStats/PDMats.jl/pull/146.