Making AbstractPDMat inherit from AbstractMatrix

[theogf]

Following issue #40
Implemented getindex for all types for a nicer display of the PDMats

[theogf]

So following https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array-1
I realized I should implement size and setindex!.
But what do you think setindex! should do? :

• Should it just return an error (we assume all AbstractPDMat are entirely immutable) ?
• Should it allow to change diagonal elements (recomputing Cholesky afterwards) ?
• Should it allow to change off-diagonal elements and automatically symmetrize the matrix accordingly (recomputing Cholesky afterwards) ?

[jlumpe]

I thin you can leave setindex! undefined for non-mutable arrays. size should be implemented for AbstractArrays.

[theogf]

I added size and left out setindex!, any more details I should take care?
Note that the advantages of inheriting from AbstractMatrix also involves having fallback algebraic operations (i.e Base.:-(PDMat,Matrix) etc), and a nicer display.

[theogf]

Bump

[theogf]

Bump

[andreasnoack]

@theogf Why were these / methods added? They seem wrong. E.g.

julia> randn(3)/randn(3,3)
ERROR: DimensionMismatch("Both inputs should have the same number of columns")
Stacktrace:
 [1] /(A::Vector{Float64}, B::Matrix{Float64})
   @ LinearAlgebra ~/julia/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/generic.jl:1143
 [2] top-level scope
   @ REPL[2]:1

[theogf]

Why were these / methods added?

Because they were missing :)

But you're right, the current definition for PDiagMat is wrong
To be fair it's also wrong for \ :

julia> x = rand(3)
3-element Array{Float64,1}:
 0.7725835504446747
 0.8135553481482751
 0.2409921677507807

julia> x' isa AbstractVecOrMat
true

julia> rand(3, 3) \ x'
ERROR: DimensionMismatch("B has leading dimension 1, but needs 3")
Stacktrace:
 [1] getrs!(::Char, ::Array{Float64,2}, ::Array{Int64,1}, ::Array{Float64,2}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/lapack.jl:996
 [2] ldiv! at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/lu.jl:391 [inlined]
 [3] \(::LU{Float64,Array{Float64,2}}, ::Adjoint{Float64,Array{Float64,1}}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/factorization.jl:101
 [4] \(::Array{Float64,2}, ::Adjoint{Float64,Array{Float64,1}}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/generic.jl:1116
 [5] top-level scope at REPL[18]:1

julia> 

[theogf]

Should it be changed to:

Suggested change

	/(x::AbstractVecOrMat, a::PDiagMat) = a.inv_diag .* x
	/(x::AbstractMatrix, a::PDiagMat) = x .* a.inv_diag'

[andreasnoack]

Ah. We allow the 1x1 case. A bit annoying that we allow that but it's too late to change.

[oschulz]

Thanks for this @theogf !