Julia: transpose for Symmetric is not recursive

[tkoolen]

From https://github.com/JuliaArrays/StaticArrays.jl/pull/358#issuecomment-358525234:

julia> a = Symmetric([[row * col * i for i = 1 : 2] for row = 1 : 3, col = 1 : 3])
3×3 Symmetric{Array{Int64,1},Array{Array{Int64,1},2}}:
 [1, 2]  [2, 4]   [3, 6] 
 [2, 4]  [4, 8]   [6, 12]
 [3, 6]  [6, 12]  [9, 18]

julia> transpose(a)[1]
2-element Array{Int64,1}:
 1
 2

julia> transpose(Array(a))[1]
1×2 Transpose{Int64,Array{Int64,1}}:
 1  2

[martinholters]

Unfortunate titlecasing here (EDIT: Hehe, funny timing) because Transpose is, transpose isn't:

julia> transpose(a)[1]
2-element Array{Int64,1}:
 1
 2

julia> Transpose(a)[1]
1×2 Transpose{Int64,Array{Int64,1}}:
 1  2

But maybe for a Symmetric matrix we'd want a[i,j] == transpose(a[j,i])? While this looks weird a first glance, for adjoint and Hermitian, we have the same situation, and maybe the condition for Hermitian should be a[i,j] == adjoint(a[j,i]) instead of a[i,j] == conj(a[j,i])? Or, put it another way, shouldn't a == Transpose(a) for a::Symmetric?

[fredrikekre]

https://github.com/JuliaLang/julia/blob/c4b21ca160e40ec61d735487c5e03f6c8b93a53b/base/linalg/symmetric.jl#L248-L249

Since we have that^ I think one shouldn't wrap things in Symmetric that isn't invariant under (recursive) transposition.

[martinholters]

Or Symmetric would need to take care of the recursive transpose. E.g. for the a above, maybe we actually do want this:

julia> a[1,2]
2-element Array{Int64,1}:
 2
 4

julia> a[2,1]
1×2 Transpose{Int64,Array{Int64,1}}:
 2  4

Problems: What should eltype(a) be and what to do about the main diagonal?

[tkoolen]

Problems: What should eltype(a) be and what to do about the main diagonal?

Yes, this seems like a real problem, as well as the getindex type instability arising from a[i,j] == transpose(a[j,i]).

I think one shouldn't wrap things in Symmetric that isn't invariant under (recursive) transposition.

That seems hard to check for automatically, correctly and efficiently. Also, even though I personally haven't used Symmetric to wrap types that aren't invariant under transposition, I suspect there may be legitimate use cases for e.g. Symmetric{T, Matrix{SMatrix{2, 2, T, 4}}}.

I think that simply making it so that Array(transpose(a)) == transpose(Array(a)) for a::Symmetric (i.e., by changing the transpose(A::Symmetric) method) would be the least surprising solution for users.

[tkoolen]

On, the other hand,

shouldn't a == Transpose(a) for a::Symmetric

is a good point. Maybe the eltype of a should still just be the eltype of its parent matrix, but getindex should always convert to eltype(a)? Then the case in the issue description would error (and maybe that would be fine), but Symmetric{T, Matrix{SMatrix{2, 2, T, 4}}} would adhere to a == Transpose(a), and be pretty efficient and type stable.