JuliaStats / PDMats.jl

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

Avoid Disrupting Array Types #177

Closed ParadaCarleton closed 1 year ago

ParadaCarleton commented 1 year ago

As an example:

julia> using DimensionalData, Distributions

julia> z = randn(3, 3) |> x->x' * x |> x->DimArray(x, (ax1=[:x, :y, :z], ax2=[:x, :y, :z]))
3×3 DimArray{Float64,2} with dimensions: 
  Dim{:ax1} Categorical{Symbol} Symbol[:x, :y, :z] ForwardOrdered,
  Dim{:ax2} Categorical{Symbol} Symbol[:x, :y, :z] ForwardOrdered
        :x         :y         :z
  :x   3.2005    -0.550321  -0.538308
  :y  -0.550321   2.14419   -2.98709
  :z  -0.538308  -2.98709    5.92229

julia> typeof(MvNormal(z))
ZeroMeanFullNormal{Tuple{OneTo{Int64}}} (alias for MvNormal{Float64, PDMats.PDMat{Float64, Array{Float64, 2}}, FillArrays.Zeros{Float64, 1, Tuple{Base.OneTo{Int64}}}})

julia> rand(MvNormal(z))
3-element Vector{Float64}:
 -0.14355464178021315
  1.076531030969517
 -0.666861987973966

This makes it impossible to use DimensionalData.jl or other array wrappers (such as CuArrays) together with PDMats.jl.

devmotion commented 1 year ago

Distributions is not needed for the example, a bit more minimal and clearer:

julia> PDMat(z)
3×3 PDMat{Float64, Matrix{Float64}}:
  1.87205    3.20828  -0.993003
  3.20828    5.97114  -1.88158
 -0.993003  -1.88158   0.597444

I'm not sure if it's an actual issue with PDMats though, the main problem seems to be

julia> using LinearAlgebra

julia> cholesky(z)
Cholesky{Float64, Matrix{Float64}}
U factor:
3×3 UpperTriangular{Float64, Matrix{Float64}}:
 1.36823  2.34484   -0.725759
  ⋅       0.687651  -0.261456
  ⋅        ⋅         0.0485701
devmotion commented 1 year ago

Seems to be https://github.com/rafaqz/DimensionalData.jl/issues/514/.