Inverse of `Symmetric{ , Diagonal}` returns a dense matrix

[vandenman]

using LinearAlgebra

julia> A = I(2)
2×2 Diagonal{Bool, Vector{Bool}}:
 1  ⋅
 ⋅  1

julia> inv(A) # correct, diagonal matrix
2×2 Diagonal{Float64, Vector{Float64}}:
 1.0   ⋅ 
  ⋅   1.0

julia> inv(Symmetric(A)) # dense matrix?
2×2 Matrix{Float64}:
 1.0  0.0
 0.0  1.0

possible fix:

inv(A::Symmetric{<:Any, <:Diagonal}) = inv(A.data)

I'd be happy to open a PR if this seems reasonable.

[araujoms]

There are few methods that hit such doubly-wrapped matrices. But I think adding this one is fine. You should just do instead

inv(A::Symmetric{<:Any,<:Diagonal}) = Symmetric(inv(A.data), sym_uplo(A.uplo))
inv(A::Hermitian{<:Any,<:Diagonal}) = Hermitian(inv(A.data), sym_uplo(A.uplo))

in order to preserve the wrapper, and get the analogous method for Hermitian.

[stevengj]

inv(A::Hermitian{<:Any,<:Diagonal}) = Hermitian(inv(A.data), sym_uplo(A.uplo))

This isn't correct, because Hermitian(::Diagonal) ignores the imaginary parts:

julia> D = Diagonal(randn(ComplexF64, 3))
3×3 Diagonal{ComplexF64, Vector{ComplexF64}}:
 0.535383-0.0897475im          ⋅                    ⋅    
          ⋅            0.90297+0.710358im           ⋅    
          ⋅                    ⋅           -1.33897+0.465176im

julia> A = Hermitian(D)
3×3 Hermitian{ComplexF64, Diagonal{ComplexF64, Vector{ComplexF64}}}:
 0.535383+0.0im          ⋅               ⋅    
          ⋅      0.90297+0.0im           ⋅    
          ⋅              ⋅      -1.33897+0.0im

julia> inv(A) ≈ Hermitian(inv(D))
false

julia> inv(A) == Hermitian(inv(real(D)))
true

A correct method (that avoids unnecessary allocations) would be:

inv(A::Hermitian{<:Number,<:Diagonal}) = Hermitian(Diagonal(inv.(real.(A.data.diag))), sym_uplo(A.uplo))

To make it work for <:Any I think you need to replace real. with hermitianpart., but hermitianpart is currently not defined for scalars, annoyingly (#1438).

[araujoms]

That's true, but you want instead

inv(A::Hermitian{<:Any,<:Diagonal}) = Hermitian(inv(Diagonal(hermitianpart.(A.data.diag))), sym_uplo(A.uplo))

in order to hit the inv(::Diagonal) method.

[stevengj]

Hitting the inv(::Diagonal) method requires an extra allocation, so isn't it better to fuse the two steps, eliminating the inv(::Diagonal) call? On the other hand, maybe we don't care too much here as long as it is linear time, and the inv(::Diagonal) method has some additional checks that we wouldn't want to have to replicate here.

[vandenman]

sounds reasonable, but the Diagonal method already throws a SingularException when an element is zero, and then we'd need to duplicate the logic?

Alternatively, we could define inv!(D::Diagonal) and call that?

I'd be happy to provide a pr either way.

[stevengj]

In either case, you currently need to define two methods due to #1438:

inv(A::Hermitian{<:Number,<:Diagonal}) = Hermitian(inv(Diagonal(real(A.data.diag))), sym_uplo(A.uplo))
inv(A::Hermitian{<:Any,<:Diagonal}) = Hermitian(inv(Diagonal(hermitianpart.(A.data.diag))), sym_uplo(A.uplo))

[stevengj]

Actually, hermitianpart doesn't work here either because Hermitian(::Diagonal) doesn't actually call hermitianpart on the diagonal entries, it just calls Hermitian on them. So you'd want something like:

inv(A::Hermitian{<:AbstractMatrix,<:Diagonal}) = Hermitian(inv(Diagonal(Hermitian.(A.data.diag, sym_uplo(A.uplo)))), sym_uplo(A.uplo))

though this then produces a Hermitian wrapper around a Diagonal wrapper around a Hermitian wrapper, ugh, which is probably not quite what we want either.

[araujoms]

On a second thought, taking the inverse of a symmetric/Hermitian matrix of matrices is a rare enough operation, that currently Julia does not support, so I think you can leave it out with a clear conscience. Just add the methods

inv(A::Symmetric{<:Number,<:Diagonal}) = Symmetric(inv(A.data), sym_uplo(A.uplo))
inv(A::Hermitian{<:Number,<:Diagonal}) = Hermitian(inv(real(A.data)), sym_uplo(A.uplo))

and be happy.