How to add two symmetric matrices inplace?

General Usage
,
1

Symmetric matrices do not permit editing off-diagonal elements. I have two matrices S1 and S2 and I’m trying to evaluate S1 + S2 and save the result into S2. Typically S2 .+= S1 would achieve this but this doesn’t work for Symmetric matrices.

julia> S1 = Symmetric(ones(3,3), :U);

julia> S2 = similar(S1);

julia> S2 .+= S1
ERROR: ArgumentError: Cannot set a non-diagonal index in a symmetric matrix
Stacktrace:
  [1] setindex!
    @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/symmetric.jl:225 [inlined]
  [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, Matrix{Float64}}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2}, Nothing, typeof(+), Tuple{Symmetric{Float64, Matrix{Float64}}, Symmetric{Float64, Matrix{Float64}}}})
    @ Base.Broadcast ./broadcast.jl:890
 [10] top-level scope
    @ REPL[24]:1

Strangely assignment works:

julia> S2 .= S1
3×3 Symmetric{Float64, Matrix{Float64}}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

julia> all(isone, S2)
true

so perhaps one may set off-diagonal elements after all?

I realize that there’s a way around this by using the parent arrays:

julia> parent(S2) .= 0

julia> parent(S2) .+= parent(S1);

julia> S2
3×3 Symmetric{Float64, Matrix{Float64}}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

but this seems to be an uncessary foray into the internals.

Also

julia> S2 += S1
3×3 Symmetric{Float64, Matrix{Float64}}:
 3.0  2.0  2.0
 2.0  2.0  2.0
 2.0  2.0  2.0

works as expected, but this creates a new binding. This is not what I want.

2

Related issues: #38055, #32722