Julia: Symmetric(A::Symmetric) should be a no-op?

[fredrikekre]

I find it surprising that this is the case for Symmetric (and Hermitian):

julia> typeof(Symmetric(Symmetric(rand(2, 2))))
Symmetric{Float64,Symmetric{Float64,Array{Float64,2}}}

instead of just Symmetric{Float64,Array{Float64,2}}. Is this behaviour ever useful?

The simplest solution would be to define these:

Symmetric(A::Symmetric, uplo::Symbol=:U) = A
Hermitian(A::Hermitian, uplo::Symbol=:U) = A

and just ignore the uplo. A more general fix would be the following, but that will sometimes create a new Symmetric/Hermitian instance, depending on the input uplo:

function Symmetric(A::Symmetric{T,S}, uplo::Symbol=:U) where {T,S}
    if A.uplo == char_uplo(uplo)
        return A
    else
        transpose!(A.data, copy(A.data)) # to keep correct aliasing
        return Symmetric{T,S}(A.data, char_uplo(uplo))
    end
end

[andreasnoack]

This is not intentional so a PR would great. I think the second version is the right direction. We'll have to make a decision about aliasing though. Your version will alias the output depending on the value of uplo and usually (but not always) we don't want aliasing to be value dependent. Usually, we also don't want contructors to modify their argument and we don't want to make unnecessary copies. We could consider throwing if the uplos don't match. If the wrapped matrix can be Symmetric already then I think it is unlikely that the user would care about the setting the value of uplo.

[fredrikekre]

Not sure I understood you correctly, but will the second version not always alias the .data field?

julia> A = [1. 2; 3 4]; ASym = Symmetric(A);

julia> pointer(Symmetric(ASym, :U).data)
Ptr{Float64} @0x00007f73837f6890

julia> pointer(Symmetric(ASym, :L).data)
Ptr{Float64} @0x00007f73837f6890

Usually, we also don't want constructors to modify their argument and we don't want to make unnecessary copies.

I agree, so my new proposal would be the following:

Symmetric(A::Symmetric) = A
Symmetric(A::Symmetric, uplo::Symbol) = A.uplo == char_uplo(uplo) ? A : throw(ArgumentError("something"))

Will submit a PR soonish.

[andreasnoack]

You are right. I misread it. Btw, you can use the unexported LinAlg.copytri! for this. Great with a PR.

[fredrikekre]

Yea, that would have been cleaner, but would not work with for example Symmetric{T, SparseMatrixCSC}, or user defined matrices. So I think throwing if the uplos don't match is the way to go here.

[andyferris]

I've wondered for a while whether uplo should be a type var. it would make matching uplos into a dispatch problem, here. And uplo is usually set by programmer intdnt, not dynamic data, so it should be safe AFAICT.

[TotalVerb]

@andyferris That has the unfortunate consequence of compiling two copies of each function for Symmetric and Hermitian.

[andreasnoack]

I don't think it is worth it. The example here is very much a corner case. I don't think we have seen examples of more common use cases where it would be useful to dispatch on uplo

[andyferris]

It could help slightly for StaticArrays, I suppose. However, in that case a symmetric matrix that is more compact (only storing the necessary numbers) would be a bit better, if we are worried about this level of optimization.