In-place broadcast on Symmetric and Hermition fails against constants

[ChrisRackauckas]

A = Symmetric(Diagonal(rand(2)))
A .= A ./ 2

ERROR: ArgumentError: Cannot set a non-diagonal index in a symmetric matrix
Stacktrace:
 [1] setindex! at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\symmetric.jl:208 [inlined]
 [2] _setindex! at .\abstractarray.jl:1050 [inlined]
 [3] setindex! at .\abstractarray.jl:1020 [inlined]
 [4] macro expansion at .\broadcast.jl:843 [inlined]
 [5] macro expansion at .\simdloop.jl:73 [inlined]
 [6] copyto! at .\broadcast.jl:842 [inlined]
 [7] copyto! at .\broadcast.jl:797 [inlined]
 [8] materialize!(::Symmetric{Float64,Diagonal{Float64,Array{Float64,1}}}, ::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2},Nothing,typeof(/),Tuple{Symmetric{Float64,Diagonal{Float64,Array{Float64,1}}},Int64}}) at .\broadcast.jl:756
 [9] top-level scope at none:0

though that should be harmless since if it was symmetric before, that would have to keep the same values in the end. Same happens with Hermitian.

[jarlebring]

Two debatable results in your example are at play:

1. A ./ 2 materializes a Matrix if A isa Symmetric, although it could be a Symmetric. Note: A / 2 does not materialize a Matrix.
2. If B is a diagonal Matrix and A::Symmetric{Diagonal} the assignment A .= B will fail

Point 1 in more minimal example:

julia> A=Symmetric(randn(2,2))
2×2 Symmetric{Float64,Array{Float64,2}}:
 0.174594  0.441573
 0.441573  1.67689 

julia> A ./2
2×2 Array{Float64,2}:
 0.087297  0.220786
 0.220786  0.838447

julia> A / 2
2×2 Symmetric{Float64,Array{Float64,2}}:
 0.087297  0.220786
 0.220786  0.838447

Point 2 in more minimal example:

julia> A = Symmetric(Diagonal(rand(2)))
julia> A .= diagm(0=>randn(2))
ERROR: ArgumentError: Cannot set a non-diagonal index in a symmetric matrix
Stacktrace:
 [1] setindex! at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/symmetric.jl:208 [inlined]
 [2] _setindex! at ./abstractarray.jl:1050 [inlined]
 [3] setindex! at ./abstractarray.jl:1020 [inlined]
 [4] copyto!(::Symmetric{Float64,Diagonal{Float64,Array{Float64,1}}}, ::Array{Float64,2}) at ./multidimensional.jl:877
 [5] copyto! at ./broadcast.jl:838 [inlined]
 [6] copyto! at ./broadcast.jl:797 [inlined]
 [7] materialize!(::Symmetric{Float64,Diagonal{Float64,Array{Float64,1}}}, ::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2},Nothing,typeof(identity),Tuple{Array{Float64,2}}}) at ./broadcast.jl:756
 [8] top-level scope at none:0

If one of the two is fixed, your problem should be solved. Which of the two is the most unnatural? I would say point 2 is more difficult to fix without considerable additional computation, i.e. checking if a Matrix is diagonal although it is somehow handled when A isa Diagonal and A .= diagm(0=>...).

[jishnub]

It is also a bit surprising that while broadcast assignment works for identical Symmetric matrices on the LHS and the RHS:

julia> A = Symmetric(Diagonal(rand(2)))
2×2 Symmetric{Float64, Diagonal{Float64, Vector{Float64}}}:
 0.846732  0.0
 0.0       0.0850426

julia> A .= A
2×2 Symmetric{Float64, Diagonal{Float64, Vector{Float64}}}:
 0.846732  0.0
 0.0       0.0850426

an in-place sum doesn't:

julia> A .+= A
ERROR: ArgumentError: Cannot set a non-diagonal index in a symmetric matrix
Stacktrace:
  [1] setindex!(A::Symmetric{Float64, Diagonal{Float64, Vector{Float64}}}, v::Float64, i::Int64, j::Int64)
    @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/symmetric.jl:225
  [2] _setindex!
    @ ./abstractarray.jl:1297 [inlined]
  [3] setindex!
    @ ./abstractarray.jl:1267 [inlined]
  [4] macro expansion
    @ ./broadcast.jl:983 [inlined]
  [5] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [6] copyto!
    @ ./broadcast.jl:982 [inlined]
  [7] copyto!
    @ ./broadcast.jl:935 [inlined]
  [8] materialize!
    @ ./broadcast.jl:893 [inlined]
  [9] materialize!(dest::Symmetric{Float64, Diagonal{Float64, Vector{Float64}}}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2}, Nothing, typeof(+), Tuple{Symmetric{Float64, Diagonal{Float64, Vector{Float64}}}, Symmetric{Float64, Diagonal{Float64, Vector{Float64}}}}})
    @ Base.Broadcast ./broadcast.jl:890
 [10] top-level scope
    @ REPL[14]:1

Surely the same indices are being set?

[LorenzoFioroni]

It seems to me that Symmetric matrices would nicely fit within the StructuredMatrix framework

LinearAlgebra.jl/src/structuredbroadcast.jl

Line 11 in 2c3fe9b

const StructuredMatrix{T} = Union{Diagonal{T},Bidiagonal{T},SymTridiagonal{T},Tridiagonal{T},LowerTriangular{T},UnitLowerTriangular{T},UpperTriangular{T},UnitUpperTriangular{T}}

In this way broadcasting could be extended in the same way it is done, e.g., for LowerTriangular matrices.
Is there a drawback to this approach? Why isn't Symmetric included in the definition of StructuredMatrix?

[araujoms]

I think StructuredMatrix was originally conceived to handle matrices that have structural zeros. But no, there's nothing stopping the extension, in fact @mbauman said that this is how it should be done: #773 (comment)
Someone just needs to do it.