Why do we need the `T` in Symmetric{T, S}?

General Usage
, ,
1

I am a bit curious about the question in the title. Since it is unanswered in the docstring, I list its definition as a reference.

# Symmetric and Hermitian matrices
struct Symmetric{T,S<:AbstractMatrix{<:T}} <: AbstractMatrix{T}
    data::S
    uplo::Char

    function Symmetric{T,S}(data, uplo::Char) where {T,S<:AbstractMatrix{<:T}}
        require_one_based_indexing(data)
        (uplo != 'U' && uplo != 'L') && throw_uplo()
        new{T,S}(data, uplo)
    end
end

Since Symmetric is merely a “shell”, I wonder what is the purpose of T (the design that singles it out). Consider the following example

julia> using LinearAlgebra

julia> function m() return rand(-9:.1:9, 2, 2) end
m (generic function with 1 method)

julia> A = m()
2Ă—2 Matrix{Float64}:
 1.9   3.7
 8.7  -1.8

julia> LinearAlgebra.symmetric_type(typeof(A[1, 1]))
Float64

julia> LinearAlgebra.symmetric_type(typeof(A))
Symmetric{Float64, Matrix{Float64}}

julia> B = [m() for i in 1:2, j in 1:2]
2Ă—2 Matrix{Matrix{Float64}}:
 [8.6 2.1; 0.8 8.9]   [0.7 -8.4; -2.3 -2.8]
 [6.4 4.1; -5.8 1.7]  [4.6 -3.0; 4.5 5.6]

julia> LinearAlgebra.symmetric_type(typeof(B[1, 1]))
Symmetric{Float64, Matrix{Float64}}

julia> LinearAlgebra.symmetric_type(typeof(B))
Symmetric{AbstractMatrix, Matrix{Matrix{Float64}}}

See the last line, AbstractMatrix is hardly instructive.
I mean, what about

image
image685Ă—463 12 KB

The related code is here

(and line 60 therein).

2

T here is the type of the elements inside of the Symmetric, indicated by <: AbstractMatrix{T} (where that same T designates the element type).

Are you perhaps looking for eltype?

1 Like
3

Yes. But then, isn’t this T redundant? making the type expression looks cumbersome. (I’ve revised the title of the topic). Is the T indispensable? We could have infer/dispatch conveniently, e.g. with the aid of eltype you suggested.

4

You need the T so you can write <: AbstractMatrix{T} in the definition. You can’t just do struct Symmetric{S} or struct Symmetric{S<:AbstractMatrix{T}} (since this last one needs T to be defined first).

3 Likes
5

Why not do this

julia> struct MySymmetric{S <: AbstractMatrix{T} where T <: Union{Number, AbstractMatrix}}
           data::S
       end

julia> using LinearAlgebra

julia> A = rand(2, 2); B = [rand(2, 2) for i in 1:2, j in 1:2];

julia> MySymmetric(A)
MySymmetric{Matrix{Float64}}([0.26880897363440703 0.6476423313174478; 0.6105024571886544 0.1442311885018166])

julia> LinearAlgebra.Symmetric(A)
2Ă—2 Symmetric{Float64, Matrix{Float64}}:
 0.268809  0.647642
 0.647642  0.144231

julia> MySymmetric(B)
MySymmetric{Matrix{Matrix{Float64}}}([[0.8054973493291133 0.7597395285524011; 0.14616471160829714 0.42034533671996954] [0.737336878028752 0.8408610455649348; 0.48620788006404236 0.15674775620016002]; [0.8856275147903988 0.7373043349918851; 0.3587458066422289 0.32907006075331047] [0.8815327773307925 0.8533663448400511; 0.15092667149027095 0.18201636320297676]])

julia> LinearAlgebra.Symmetric(B)
2Ă—2 Symmetric{AbstractMatrix, Matrix{Matrix{Float64}}}:
 [0.805497 0.75974; 0.75974 0.420345]    [0.737337 0.840861; 0.486208 0.156748]
 [0.737337 0.486208; 0.840861 0.156748]  [0.881533 0.853366; 0.853366 0.182016]
6 
julia> typeof(MySymmetric(A)) <: AbstractMatrix{Float64}
false

julia> typeof(Symmetric(A)) <: AbstractMatrix{Float64}
true

So you want to declare

struct MySymmetric{S <: AbstractMatrix{T}} <: AbstractMatrix{T} end

but that doesn’t work because T is not defined. You can fix it by adding the extra type parameter

struct MySymmetric{T, S <: AbstractMatrix{T}} <: AbstractMatrix{T} end
4 Likes
7

You can define that with a local T

struct MySymmetric{S<:AbstractMatrix} end
# or equivalently
struct MySymmetric{S<:AbstractMatrix{T where T}} end

but the main problem is then still that MySymmetric{Matrix{Float64} is never a subtype of AbstractMatrix{Float64} as one can see in the first typeof in the post above me. And that is something one wants for nice default dispatches for all eltypes.

3 Likes
8

Yes, your last line is good. I copy-paste it here

struct MySymmetric{T, S <: AbstractMatrix{T}} <: AbstractMatrix{T} end

But the LinearAlgebra.jl define it as

struct Symmetric{T,   S <: AbstractMatrix{<:T}} <: AbstractMatrix{T} end

We see that the latter uses <: T, what is the notivation? so that

Looks a bit inconsistent.

I now understand this point, thank you.

9

I’m unsure if this is intentional, probably just an oversight that has not been corrected since.

10

There is a discussion in the relevant PR but it did not clarify the issue for me, I don’t understand the possibility they mention.

1 Like
11

My understanding is that it’s because of matrices of matrices (transpose of a number has the same type, but transpose of a matrix is of type Transpose).

Example behavior:

julia> A = [1 2;
            3 4]
2Ă—2 Matrix{Int64}:
 1  2
 3  4

julia> B = [[A] [A];
            [A] [A]]
2Ă—2 Matrix{Matrix{Int64}}:
 [1 2; 3 4]  [1 2; 3 4]
 [1 2; 3 4]  [1 2; 3 4]

julia> C = Symmetric(B)
2Ă—2 Symmetric{AbstractMatrix, Matrix{Matrix{Int64}}}:
 [1 2; 2 4]  [1 2; 3 4]
 [1 3; 2 4]  [1 2; 2 4]

julia> C[1, 1]
2Ă—2 Symmetric{Int64, Matrix{Int64}}:
 1  2
 2  4

julia> C[1, 2]
2Ă—2 Matrix{Int64}:
 1  2
 3  4

julia> C[2, 1]
2Ă—2 transpose(::Matrix{Int64}) with eltype Int64:
 1  3
 2  4

(three different eltypes, all of which are <: AbstractMatrix)

3 Likes
12

Do make sense. But in this case, T is not some “eltype” but a “least common eltype”. And in your example, T is AbstractMatrix, indicating that the dispatch advantage is lost. (But I guess this usage is not practical)

13

Wait a minute. There is still a question.
The LinearAlgebra.jl defines

struct Symmetric{T, S <: AbstractMatrix{<:T}} <: AbstractMatrix{T} end

And you call C = Symmetric(B), where B isa Matrix{Matrix{Int64}}.
At the time of construction (which is “earlier”), T can surely be Matrix{Int64}.
It’s after some manipulations (at a “later” time) that julia find AbstractMatrix can serve as a “least common supertype” for all the eltypes.
The weird thing is: Did I tell julia that what role I want T to play? Did I stipulate that T should be the “least common supertype” at that “later” time? Seems I hadn’t.
Then why can’t T also be Matrix{Int64}?

14

You did not tell it explicitly. T is derived automatically by the code (not magically by Julia, by the code LinearAlgebra.jl/src/symmetric.jl at ae208d6aae535461ba29812bbdbf6c6dcda9f6ce · JuliaLang/LinearAlgebra.jl · GitHub).

If T was Matrix{Int64}, then C[i, j] should always return a Matrix{Int64}. That’s certainly possible, but it would involve copying when C[i, j] should be a transposed version of B[i, j]. Currently, we get just a lazy wrapper (of type Transpose{Int64, Matrix{Int64}}) – to avoid copies, I suppose.

In summary, C[i, j] can return Transpose, Symmetric, or Matrix, depending on i, j, and so AbstractMatrix is the narrowest common type (Transpose <: AbstractMatrix, too).

2 Likes
15

The assumption that construction happens earlier is wrong. For any value one wants to construct, the type must exist before the construction is eventually over (new returns). I don’t really understand the questions in your post, but one thing’s certain: the type system checks the constraints on the type parameters during type application, so before it creates the type, so before the value can be constructed. For example:

julia> using LinearAlgebra
 
julia> Symmetric{Int, Int}
ERROR: TypeError: in Symmetric, in S, expected S<:(AbstractMatrix{<:Int64}), got Type{Int64}
Stacktrace:
 [1] top-level scope
   @ REPL[2]:1
1 Like
16

Yes, it is by the code in LinearAlgebra.jl. But the link you gave is not precise, since that S is very general therefore with no chance to be called. The actually code being called should be

I propose to modify this function as my_symmetric_type as follows

julia> using LinearAlgebra

julia> @isdefined symmetric_type
false

julia> symmetric_type = LinearAlgebra.symmetric_type
symmetric_type (generic function with 4 methods)

julia> @isdefined promote_op
false

julia> promote_op = Base.promote_op
promote_op (generic function with 1 method)

julia> function my_symmetric_type(::Type{T}) where {S<:AbstractMatrix, T<:AbstractMatrix{S}}
         return Symmetric{typejoin(S, promote_op(transpose, S), symmetric_type(S)), T}
       end
my_symmetric_type (generic function with 1 method)

julia> A = rand(-9:9, 2, 2);

julia> M = [A for i in 1:2, j in 1:2]; t = typeof(M)
Matrix{Matrix{Int64}} (alias for Array{Array{Int64, 2}, 2})

julia> my_symmetric_type(t) # me proposed
Symmetric{AbstractMatrix{Int64}, Matrix{Matrix{Int64}}}

julia> LinearAlgebra.symmetric_type(t) # compared with the existing
Symmetric{AbstractMatrix, Matrix{Matrix{Int64}}}

Is it good?
also @jishnub

17

Thanks, it’s instructive.

Therefore my current guess about Symmetric{T, S} is:

  1. T is the “least common supertype” of eltype of all the entries (finally after the symmetric matrix object has been constructed). The aim of defining this is to facilitate dispatch.
  2. S is simply the type of the data stored.