eltype of Symmetric or Hermitian matrix of matrices is AbstractMatrix

[araujoms]

MWE:

using LinearAlgebra
m = Symmetric(fill(ones(2,2), 2, 2))
eltype(m)
n = Hermitian(fill(ones(2,2), 2, 2))
eltype(n)

On the other hand, if we do this with UpperTriangular or LowerTriangular the eltype is Matrix{Float64}, as expected.

It's not a burning issue, but prevents the faster kronecker product introduced in JuliaLang/julia#53186 from working on Symmetric matrices of matrices, and requires the less elegant syntax used in JuliaLang/julia#54413.

I've tried to understand the problem, but it seems to come directly from the constructor so I'm stuck.

[jishnub]

Does the result improve if the eltype is a Union of three types? That's generally the best that this may be narrowed down to.

[araujoms]

I don't know, it depends on whether promote_op gives a sensible answer when given this Union.

The issue is that I need to allocate a matrix to store the result of the kronecker product, which I do with Matrix{promote_op(*, T, S)}. Now if T,S == Matrix{Float64} this gives me a Matrix{Matrix{Float64}}, which is exactly what I need. But if T,S == AbstractMatrix this gives me a Matrix{Any} which ruins everything.

EDIT: I just tested, and it does work. Let T = Union{Symmetric{Float64, Matrix{Float64}}, Transpose{Float64, Matrix{Float64}}, Matrix{Float64}}. Then promote_op(*, T, T) == Matrix{Float64} and everything is fine. Ditto if we have Hermitian and Adjoint instead of Symmetric and Transpose.

[dkarrasch]

The reason is that that Symmetric wrappers only look at one, say the upper, triangle, and as for elements of the lower triangle, we pick elements from the upper and transpose them. Similarly for Hermitian. For Number elements, that transpose/adjoint doesn't make a type difference, but for matrix elements, the obtained result is a Transpose/Adjoint, which combined with the original matrices returns the AbstractMatrix. The function to look at is

function symmetric_type(::Type{T}) where {S<:AbstractMatrix, T<:AbstractMatrix{S}}
    return Symmetric{AbstractMatrix, T}
end

so it seems nobody had a better answer so far.

[araujoms]

Thanks for the information. AbstractMatrix still seems like a bad choice, as it erases information. Union{Matrix{T},Transpose{Matrix{T}},Symmetric{Matrix{T}}} would still make it possible to allocate the necessary matrix.

[araujoms]

On a second thought, the easiest solution would be for promote_op(*,AbstractMatrix,AbstractMatrix) to return AbstractMatrix. Does anything depend on it returning Any?

[dkarrasch]

On a second thought, the easiest solution would be for promote_op(*,AbstractMatrix,AbstractMatrix) to return AbstractMatrix.

Good question why that is not the case.

[oscardssmith]

not really, promote_op of almost any abstract type will return Any because the compiler can't prove that no one will introduce a new method.

[araujoms]

Do you mean that someone could define a new Prank subtype of AbstractMatrix, and define a method for * that takes two Pranks to a String?

I suppose that's possible, but I don't see why promote_type should indulge this behaviour. Isn't type inference supposed to be just a heuristic anyway?

[oscardssmith]

promote_op uses type inference, and while type inference is a heuristic, it is required to be a conservative heuristic. Specifically, type inference is always allowed to return a non-concrete supertype of the correct result, but it would be a very serious bug if type inference ever returned an incorrect result.

[araujoms]

Is anyone willing to narrow the eltype from AbstractMatrix to the union of three types? I'm afraid that requires touching the C code, and that's not something I do without getting paid.