make Symmetric/Hermitian recursive #25688

martinholters · 2018-01-22T16:40:59Z

Testing the waters here. I have no clue what a Symmetric(::Matrix{Matrix}) might be used for in practice, but I'd say the additional tests look like we'd want them to pass. The alternative would be to make a::Symmetric obey a == permutedims(a) and a::Hermitian obey a == conj(permutedims(a)), where especially the latter looks extremely artificial.

(I admit to having used promote_op here but that's what Transpose/Adjoint do and I'm afraid I have to be consistent with those...)

Fixes JuliaLang/LinearAlgebra.jl#497.

~~Note: This is on top of #25687 because it triggered the bug fixed therein.~~

andreasnoack · 2018-01-22T21:31:50Z

stdlib/LinearAlgebra/src/symmetric.jl

@@ -1,7 +1,7 @@
 # This file is a part of Julia. License is MIT: https://julialang.org/license

 # Symmetric and Hermitian matrices
-struct Symmetric{T,S<:AbstractMatrix{T}} <: AbstractMatrix{T}
+struct Symmetric{T,S<:AbstractMatrix{<:T}} <: AbstractMatrix{T}


Could you explain this change here?

The type of transpose(A.data[i, j]) may not be the same as that of A.data[i, j], so in order to not lie about the return type of getindex(A, i, j), eltype(A) may need to be widened with respect to eltype(A.data).

I’m not sure this is quite right. The Adjoint stuff might not even match this signature. IIRC, for Transpose and Adjoint we leave this free and check it later (in constructors, getindex, etc).

I'm confused, maybe I'm misunderstanding you. If A is the argument to the Symmetric outer constructor, then T == Union{eltype(A),promote_op(transpose, eltype(A))} should always be at least eltype(A) == eltype(S) (so that eltype(S) <: T), right?

Do you mean that it would be possible to call the inner constructor with bad parameters? That I agree with. Could make the current outer constructor an inner constructor.

I second @tkoolen's explanation and share his confusion.

This definition now clicks on this end. Thanks for the explanations! :)

tkoolen · 2018-01-23T03:13:53Z

I was just thinking through JuliaLang/LinearAlgebra.jl#497 some more, and I think the approach in this PR is better.

Case 1: `A = Symmetric([rand(T, n, n) for i = 1:m, j = 1:m])`

	This PR	JuliaLang/LinearAlgebra.jl#497
`eltype(A)`	`Union{Matrix{T}, Transpose{T, Matrix{T}}}`	`Matrix{T}`
`eltype(A)` concrete	❌	✅
`getindex(A, i, j)` allocations	O(1)	O(n^2)

Here eltype(A) is a non-isbits union, but I think this is outweighed by the getindex allocations, and perhaps the allocations for this PR can even go away entirely in the future.

Case 2: `A = Symmetric([rand(SMatrix{n, n, T}) for i = 1 : m, j = 1 : m])`

	This PR	JuliaLang/LinearAlgebra.jl#497
`eltype(A)`	`SMatrix{n, n, T, n^2}`	`SMatrix{n, n, T, n^2}`
`eltype(A)` concrete	✅	✅
`getindex(A, i, j)` allocations	0	0

At least right now, StaticArrays defines a non-lazy transpose, so transposed elements of A.data have the same type as non-transposed elements (my mistake in JuliaLang/LinearAlgebra.jl#497), and both approaches result in the same behavior.

Case 3: same as case 2, but hypothetically, `transpose(::SMatrix}` uses `Transpose`

	This PR	JuliaLang/LinearAlgebra.jl#497
`eltype(A)`	`Union{Transpose{SMatrix{n, n, T, n^2}}, SMatrix{n, n, T, n^2}}`	`SMatrix{n, n, T, n^2}`
`eltype(A)` concrete	❌	✅
`getindex(A, i, j)` allocations	0	0

I believe there would not be any allocations for the approach in this PR given the union-of-bitstypes optimization.

Case 4: `A = Symmetric([rand(SVector{n, T}) for i = 1 : m, j = 1 : m])`

	This PR	JuliaLang/LinearAlgebra.jl#497
`eltype(A)`	`Union{SMatrix{n, 1, T, n}, SVector{n, T}}`	N/A
`eltype(A)` concrete	❌	N/A
`getindex(A, i, j)` allocations	0	N/A

JuliaLang/LinearAlgebra.jl#497 doesn't do this case, because you can't convert an SMatrix{n, 1, T, n} (the transpose element type) to an SVector{n, T}.

martinholters · 2018-01-23T08:00:28Z

stdlib/LinearAlgebra/src/symmetric.jl

    end
 end
 @inline function getindex(A::Hermitian, i::Integer, j::Integer)
    @boundscheck checkbounds(A, i, j)
    @inbounds if (A.uplo == 'U') == (i < j)
        return A.data[i, j]
    elseif i == j
-        return eltype(A)(real(A.data[i, j]))
+        return convert(typeof(A.data[i, j]), real(A.data[i, j]))


Ah, this is probably wrong if A.data[i, j] is a matrix: should probably rather be Hermitian(...) then..

martinholters · 2018-01-23T08:06:05Z

Just noticed that

julia> Symmetric(hcat([[1 2; 3 4]]))
1×1 Symmetric{Union{Array{Int64,2}, Transpose{Int64,Array{Int64,2}}},Array{Array{Int64,2},2}}:
 [1 3; 2 4]

isn't quite symmetric and

julia> Hermitian(hcat([[1 2im; -2im 3]]))
1×1 Hermitian{Union{Array{Complex{Int64},2}, Adjoint{Complex{Int64},Array{Complex{Int64},2}}},Array{Array{Complex{Int64},2},2}}:
 [1+0im 0+0im; 0+0im 3+0im]

looks pretty wrong. Probably Symmetric and Hermitian should actually recurse on the diagonals (funny I only realized this now considering this PRs title...)

martinholters · 2018-01-24T08:28:10Z

Updated to fix these cases. I added (non-exported) functions symmetric/hermitian which mimic (lowercase) transpose / adjoint in that they eagerly evaluate on numbers and create lazy wrappers on matrices and use thos for recursing on the diagonal elements.

tkoolen · 2018-01-24T13:43:11Z

Not sure if we should care, but with 9ab815d,

julia> data = [[Complex(row * i, col * i + 1) for i = 1 : 2] for row = 1 : 3, col = 1 : 3]
3×3 Array{Array{Complex{Int64},1},2}:
 [1+2im, 2+3im]  [1+3im, 2+5im]  [1+4im, 2+7im]
 [2+2im, 4+3im]  [2+3im, 4+5im]  [2+4im, 4+7im]
 [3+2im, 6+3im]  [3+3im, 6+5im]  [3+4im, 6+7im]

julia> Symmetric(data)
ERROR: MethodError: no method matching symmetric_type(::Type{Array{Complex{Int64},1}})
Closest candidates are:
  symmetric_type(::Type{T<:AbstractArray{S,2}}) where {S, T<:AbstractArray{S,2}} at /Users/twan/code/julia/julia/usr/share/julia/site/v0.7/LinearAlgebra/src/symmetric.jl:51
  symmetric_type(::Type{T<:Number}) where T<:Number at /Users/twan/code/julia/julia/usr/share/julia/site/v0.7/LinearAlgebra/src/symmetric.jl:53
Stacktrace:
 [1] symmetric_type(::Type{Array{Array{Complex{Int64},1},2}}) at /Users/twan/code/julia/julia/usr/share/julia/site/v0.7/LinearAlgebra/src/symmetric.jl:51
 [2] Type at /Users/twan/code/julia/julia/usr/share/julia/site/v0.7/LinearAlgebra/src/symmetric.jl:45 [inlined]
 [3] Symmetric(::Array{Array{Complex{Int64},1},2}) at /Users/twan/code/julia/julia/usr/share/julia/site/v0.7/LinearAlgebra/src/symmetric.jl:44
 [4] top-level scope

julia> Hermitian(data)
ERROR: MethodError: no method matching hermitian_type(::Type{Array{Complex{Int64},1}})
Closest candidates are:
  hermitian_type(::Type{T<:AbstractArray{S,2}}) where {S, T<:AbstractArray{S,2}} at /Users/twan/code/julia/julia/usr/share/julia/site/v0.7/LinearAlgebra/src/symmetric.jl:102
  hermitian_type(::Type{T<:Number}) where T<:Number at /Users/twan/code/julia/julia/usr/share/julia/site/v0.7/LinearAlgebra/src/symmetric.jl:104
Stacktrace:
 [1] hermitian_type(::Type{Array{Array{Complex{Int64},1},2}}) at /Users/twan/code/julia/julia/usr/share/julia/site/v0.7/LinearAlgebra/src/symmetric.jl:102
 [2] Type at /Users/twan/code/julia/julia/usr/share/julia/site/v0.7/LinearAlgebra/src/symmetric.jl:96 [inlined]
 [3] Hermitian(::Array{Array{Complex{Int64},1},2}) at /Users/twan/code/julia/julia/usr/share/julia/site/v0.7/LinearAlgebra/src/symmetric.jl:95
 [4] top-level scope

andreasnoack · 2018-01-24T14:03:52Z

Does it make sense to talk about a symmetric matrix of vectors? You also used that example in JuliaLang/LinearAlgebra.jl#497 so is it a real use case?

martinholters · 2018-01-24T14:44:16Z

What should be on the diagonal in that case? What would be a reasonable transformation to turn [1+2im, 2+3im] into some x such that x == transpose(x) holds?

tkoolen · 2018-01-24T16:13:08Z

I don't have a use case for that, no. It's not completely inconceivable that someone would want it, but I'm fine with not supporting it, as long as it's thought through and documented.

Here's a more realistic case though: as it stands, wrapping a matrix of e.g. JuMP.AffineExpr (which is not <:Number) in a Symmetric wouldn't work out of the box. Side note: I personally believe that JuMP.AbstractScalar should be a Number and I am aware of JuMP's own handling of symmetric matrix variables, but still believe that there are legitimate use cases for this when setting up an optimization problem by calling generic code.

Should that be handled by fallbacks for symmetric_type and symmetric in LinearAlgebra, or should JuMP provide specializations for JuMP.AbstractScalar?

In the former case, (basically widening the ::Number methods of symmetric_type and symmetric to ::Any), the data::AbstractMatrix{<:AbstractVector} case would no longer result in a MethodError straight away, but then @martinholters' question about what should be on the diagonal still stands. Perhaps construction should fail with a descriptive error message, or maybe Symmetric's transpose/adjoint should only recurse for data::AbstractMatrix{<:AbstractMatrix}?

Sacha0 · 2018-01-24T19:19:07Z

stdlib/LinearAlgebra/src/symmetric.jl

@@ -162,11 +181,17 @@ end
 similar(A::Union{Symmetric,Hermitian}, ::Type{T}, dims::Dims{N}) where {T,N} = similar(parent(A), T, dims)

 # Conversion
-Matrix(A::Symmetric) = copytri!(convert(Matrix, copy(A.data)), A.uplo)
+function Matrix(A::Symmetric)
+    B = copytri!(convert(Matrix, copy(A.data)), A.uplo)


Can convert(Matrix, copy(A.data)) be simplified to Matrix(A.data)?

That sounds plausible, but it was like that before, so I don't know what considerations led to the copy+convert construct.

Sacha0

Technically (i.e. not speaking to broader design questions) this looks great! Thanks @martinholters! :)

andreasnoack · 2018-01-25T08:49:30Z

I think it is fair to require that if a new type doesn't want to be a subtype of, say, Number, then it would have to define more methods. The alternative is to define method Any fallback but they'd be generally wrong until the user overwrites the behavior. I think an error is better than silently returning the wrong result in some cases.

martinholters · 2018-01-25T11:02:15Z

symmetric(A::T, ::Symbol) where {T} = convert(T, 0.5*(A + transpose(A)))
symmetric_type(::Type{T}) where {T} = T

might make a reasonable fallback. (Similar for hermitian.)

But more generally, I'm doubtful whether we should do this. What we get is that instead of

S::Symmetric == permutedims(S)
H::Hermitian == conj(permutedims(H))

we get

S::Symmetric == transpose(S)
H::Hermitian == adjoint(H)

Now at first glance, it looks like a no-brainer we want this. But it comes at the expense of additional complexity. Further, if we go with this change, we have nothing to express the former behavior. So before I try to wrap this up by adding docs---are we sure we prefer the latter behavior? Is that more useful in practice? Is it worth it?

tkoolen · 2018-01-25T20:54:24Z

I don't have a strong opinion either way. The original issue didn't arise from a practical use case on my end, but from a perceived discrepancy between Matrix and Symmetric that showed up when I was coding up JuliaArrays/StaticArrays.jl#358. Maybe @andyferris or @stevengj have stronger opinions?

andreasnoack · 2018-01-25T21:00:28Z

@martinholters I'm pretty sure I prefer the latter version.

andyferris · 2018-01-29T01:48:43Z

Yes, the latter looks correct to me.

martinholters · 2018-01-29T08:59:46Z

Ok, I've added docstrings which I hope are enough documentation here. The need to create a Symmetric view of something where the element type is neither AbstractMatrix nor Number should be rare. If you try to do anyway, you get a method error. Looking at the docs of that method point you in the right direction. That said, ideas for improvement of the docs are more than welcome!

martinholters · 2018-01-30T14:19:01Z

Should this be nodded off by someone else or is good to go?

andreasnoack approved these changes Jan 22, 2018

View reviewed changes

martinholters commented Jan 23, 2018

View reviewed changes

martinholters force-pushed the mh/symmetric branch from e77b3a5 to 9ab815d Compare January 24, 2018 08:25

Sacha0 reviewed Jan 24, 2018

View reviewed changes

Sacha0 reviewed Jan 25, 2018

View reviewed changes

martinholters added 2 commits January 29, 2018 09:54

Make Symmetric/Hermitian recursive

429ce96

Add docstrings for symmetric/hermitian

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.

GPG key ID: 4AEE18F83AFDEB23
Expired

Verified
Learn about vigilant mode

04ea22e

martinholters force-pushed the mh/symmetric branch from 9ab815d to 04ea22e Compare January 29, 2018 08:55

martinholters changed the title ~~RFC: make Symmetric/Hermitian recursive~~ make Symmetric/Hermitian recursive Jan 29, 2018

andreasnoack merged commit c62558c into master Jan 30, 2018

andreasnoack deleted the mh/symmetric branch January 30, 2018 14:59

fredrikekre mentioned this pull request Feb 9, 2018

fix getindex(::HermOrSym, i, j) when i == j #25976

Merged

KristofferC mentioned this pull request Sep 24, 2018

avoid conversion from Char to String to Symbol in some LinearAlgebra routines #29352

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

make Symmetric/Hermitian recursive #25688

make Symmetric/Hermitian recursive #25688

martinholters commented Jan 22, 2018 •

edited

Loading

andreasnoack Jan 22, 2018

tkoolen Jan 23, 2018 •

edited

Loading

andyferris Jan 23, 2018

tkoolen Jan 23, 2018

martinholters Jan 23, 2018

Sacha0 Jan 24, 2018

tkoolen commented Jan 23, 2018

martinholters Jan 23, 2018

martinholters commented Jan 23, 2018

martinholters commented Jan 24, 2018 •

edited

Loading

tkoolen commented Jan 24, 2018

andreasnoack commented Jan 24, 2018

martinholters commented Jan 24, 2018

tkoolen commented Jan 24, 2018

Sacha0 Jan 24, 2018

martinholters Jan 25, 2018

Sacha0 left a comment

andreasnoack commented Jan 25, 2018

martinholters commented Jan 25, 2018

tkoolen commented Jan 25, 2018

andreasnoack commented Jan 25, 2018

andyferris commented Jan 29, 2018

martinholters commented Jan 29, 2018

martinholters commented Jan 30, 2018

make Symmetric/Hermitian recursive #25688

make Symmetric/Hermitian recursive #25688

Conversation

martinholters commented Jan 22, 2018 • edited Loading

Choose a reason for hiding this comment

tkoolen Jan 23, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tkoolen commented Jan 23, 2018

Case 1: A = Symmetric([rand(T, n, n) for i = 1:m, j = 1:m])

Case 2: A = Symmetric([rand(SMatrix{n, n, T}) for i = 1 : m, j = 1 : m])

Case 3: same as case 2, but hypothetically, transpose(::SMatrix} uses Transpose

Case 4: A = Symmetric([rand(SVector{n, T}) for i = 1 : m, j = 1 : m])

Choose a reason for hiding this comment

martinholters commented Jan 23, 2018

martinholters commented Jan 24, 2018 • edited Loading

tkoolen commented Jan 24, 2018

andreasnoack commented Jan 24, 2018

martinholters commented Jan 24, 2018

tkoolen commented Jan 24, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Sacha0 left a comment

Choose a reason for hiding this comment

andreasnoack commented Jan 25, 2018

martinholters commented Jan 25, 2018

tkoolen commented Jan 25, 2018

andreasnoack commented Jan 25, 2018

andyferris commented Jan 29, 2018

martinholters commented Jan 29, 2018

martinholters commented Jan 30, 2018

martinholters commented Jan 22, 2018 •

edited

Loading

tkoolen Jan 23, 2018 •

edited

Loading

Case 1: `A = Symmetric([rand(T, n, n) for i = 1:m, j = 1:m])`

Case 2: `A = Symmetric([rand(SMatrix{n, n, T}) for i = 1 : m, j = 1 : m])`

Case 3: same as case 2, but hypothetically, `transpose(::SMatrix}` uses `Transpose`

Case 4: `A = Symmetric([rand(SVector{n, T}) for i = 1 : m, j = 1 : m])`

martinholters commented Jan 24, 2018 •

edited

Loading