Problem with eltype of a non-numeric Symmetric matrix

**lhatzis** · April 16, 2019, 5:12pm

Hi everyone,

I have defined a custom type to wrap a Python library similarly to what SymPy.jl does. So I will focus on a custom defined for SymPy but essentially the problem is the same for my custom type too.

If I define a symmetric matrix as:

using LinearAlgebra, SymPy

x = Matrix{Sym}(undef, 2, 2)
Symmetric(x)

the type of the output is 2×2 Symmetric{Any,Array{Sym,2}}, while if I defined as

Symmetric{eltype(x),typeof(x)}(x, 'L')

the return type is 2×2 Symmetric{Sym,Array{Sym,2}}.

That means that the eltype in the definition of Symmetric returns type Any. Why is this happening and how can I resolve that? Thanks!

**mbauman** · April 16, 2019, 5:10pm

You’re looking at a really old version of Julia there. You’re probably running 1.1 or somesuch — the function that’s now called is this:

github.com

JuliaLang/julia/blob/release-1.1/stdlib/LinearAlgebra/src/symmetric.jl#L48-L51


function Symmetric(A::AbstractMatrix, uplo::Symbol=:U)
    checksquare(A)
    return symmetric_type(typeof(A))(A, char_uplo(uplo))
end

Symmetric is now recursive so it needs to reason about the transpose of things — and if it can’t then it just gives up and says Any. I’m not sure what all you need to do for a custom type.

**lhatzis** · April 16, 2019, 7:10pm

I haven’t noticed that I’m looking into an older version. Thanks for pointing it out!

**tkoolen** · April 16, 2019, 8:13pm

I think the main problem is that transpose(::Sym) is not type-stable, so that Base.promote_op(transpose, Sym) == Any. promote_op is used here:

github.com

JuliaLang/julia/blob/7132c635ab8da723c0177a9d6617bf01f4e7460a/stdlib/LinearAlgebra/src/symmetric.jl#L75


symmetric(A::Number, ::Symbol) = A


"""
    symmetric_type(T::Type)


The type of the object returned by `symmetric(::T, ::Symbol)`. For matrices, this is an
appropriately typed `Symmetric`, for `Number`s, it is the original type. If `symmetric` is
implemented for a custom type, so should be `symmetric_type`, and vice versa.
"""
function symmetric_type(::Type{T}) where {S, T<:AbstractMatrix{S}}
    return Symmetric{Union{S, promote_op(transpose, S), symmetric_type(S)}, T}
end
symmetric_type(::Type{T}) where {T<:Number} = T


struct Hermitian{T,S<:AbstractMatrix{<:T}} <: AbstractMatrix{T}
    data::S
    uplo::Char


    function Hermitian{T,S}(data, uplo) where {T,S<:AbstractMatrix{<:T}}
        @assert !has_offset_axes(data)
        new{T,S}(data, uplo)

and promote_op uses Core.Compiler.return_type in its implementation.

This one case is easy enough to fix with a PR, but I’m thinking this is likely a pretty pervasive problem throughout SymPy.

**j_verzani** · April 16, 2019, 10:08pm

Thanks for pointing this out. I have a PR to fix this particular issue, but how widely it needs doing remains to be seen.