Bug Report: Error when Subtracting AbstractTriangular wrapping Symmetric/Hermitian Matrices in LinearAlgebra

[singularitti]

I encountered the following bug in the LinearAlgebra standard library:

julia> using LinearAlgebra

julia> A = hermitianpart(rand(4, 4))
4×4 Hermitian{Float64, Matrix{Float64}}:
 0.853248  0.49251   0.351079  0.671892
 0.49251   0.462168  0.578388  0.291257
 0.351079  0.578388  0.267023  0.40795
 0.671892  0.291257  0.40795   0.382695

julia> B = UpperTriangular(A)
4×4 UpperTriangular{Float64, Hermitian{Float64, Matrix{Float64}}}:
 0.853248  0.49251   0.351079  0.671892
  ⋅        0.462168  0.578388  0.291257
  ⋅         ⋅        0.267023  0.40795
  ⋅         ⋅         ⋅        0.382695

julia> B - B'
ERROR: ArgumentError: Cannot set a non-diagonal index in a Hermitian matrix
Stacktrace:
 [1] setindex!
   @ ~/.asdf/installs/julia/1.11.0/share/julia/stdlib/v1.11/LinearAlgebra/src/symmetric.jl:264 [inlined]
 [2] _setindex!
   @ ./abstractarray.jl:1443 [inlined]
 [3] setindex!
   @ ./abstractarray.jl:1413 [inlined]
 [4] copyto_unaliased!(deststyle::IndexCartesian, dest::Hermitian{…}, srcstyle::IndexCartesian, src::UpperTriangular{…})
   @ Base ./abstractarray.jl:1095
 [5] copyto!
   @ ./abstractarray.jl:1061 [inlined]
 [6] -(A::UpperTriangular{Float64, Hermitian{…}}, B::LowerTriangular{Float64, Hermitian{…}})
   @ LinearAlgebra ~/.asdf/installs/julia/1.11.0/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:849
 [7] top-level scope
   @ REPL[5]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> @which B - B'
-(A::LinearAlgebra.AbstractTriangular, B::LinearAlgebra.AbstractTriangular)
     @ LinearAlgebra ~/.asdf/installs/julia/1.11.0/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:849

The error seems to originate from the following line:

-(A::AbstractTriangular, B::AbstractTriangular) = copyto!(similar(parent(A)), A) - copyto!(similar(parent(B)), B)

This fails because you cannot copyto! a Hermitian, Symmetric, SymTridiagonal, etc.:

julia> A = Symmetric(rand(4, 4))
4×4 Symmetric{Float64, Matrix{Float64}}:
 0.412561   0.535526    0.989538   0.0304181
 0.535526   0.219967    0.0725816  0.00634629
 0.989538   0.0725816   0.44987    0.456327
 0.0304181  0.00634629  0.456327   0.781835

julia> B = UnitLowerTriangular(A)
4×4 UnitLowerTriangular{Float64, Symmetric{Float64, Matrix{Float64}}}:
 1.0         ⋅           ⋅         ⋅
 0.535526   1.0          ⋅         ⋅
 0.989538   0.0725816   1.0        ⋅
 0.0304181  0.00634629  0.456327  1.0

julia> B - B'
ERROR: ArgumentError: Cannot set a non-diagonal index in a symmetric matrix
Stacktrace:
 [1] setindex!
   @ ~/.asdf/installs/julia/1.11.0/share/julia/stdlib/v1.11/LinearAlgebra/src/symmetric.jl:258 [inlined]
 [2] _setindex!
   @ ./abstractarray.jl:1443 [inlined]
 [3] setindex!
   @ ./abstractarray.jl:1413 [inlined]
 [4] copyto_unaliased!(deststyle::IndexCartesian, dest::Symmetric{…}, srcstyle::IndexCartesian, src::UnitLowerTriangular{…})
   @ Base ./abstractarray.jl:1095
 [5] copyto!
   @ ./abstractarray.jl:1061 [inlined]
 [6] -(A::UnitLowerTriangular{Float64, Symmetric{…}}, B::UnitUpperTriangular{Float64, Symmetric{…}})
   @ LinearAlgebra ~/.asdf/installs/julia/1.11.0/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:849
 [7] top-level scope
   @ REPL[10]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> A = SymTridiagonal(dv, ev)
4×4 SymTridiagonal{Int64, Vector{Int64}}:
 1  7  ⋅  ⋅
 7  2  8  ⋅
 ⋅  8  3  9
 ⋅  ⋅  9  4

julia> B = UnitLowerTriangular(A);
 1  ⋅  ⋅  ⋅
 7  1  ⋅  ⋅
 0  8  1  ⋅
 0  0  9  1

julia> B - B'
ERROR: ArgumentError: cannot set off-diagonal entry (2, 1)
Stacktrace:
 [1] setindex!
   @ ~/.asdf/installs/julia/1.11.0/share/julia/stdlib/v1.11/LinearAlgebra/src/tridiag.jl:466 [inlined]
 [2] _setindex!
   @ ./abstractarray.jl:1443 [inlined]
 [3] setindex!
   @ ./abstractarray.jl:1413 [inlined]
 [4] copyto_unaliased!(deststyle::IndexCartesian, dest::SymTridiagonal{Int64, Vector{…}}, srcstyle::IndexCartesian, src::UnitLowerTriangular{Int64, SymTridiagonal{…}})
   @ Base ./abstractarray.jl:1095
 [5] copyto!
   @ ./abstractarray.jl:1061 [inlined]
 [6] -(A::UnitLowerTriangular{Int64, SymTridiagonal{Int64, Vector{Int64}}}, B::UnitUpperTriangular{Int64, SymTridiagonal{Int64, Vector{Int64}}})
   @ LinearAlgebra ~/.asdf/installs/julia/1.11.0/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:849
 [7] top-level scope
   @ REPL[15]:1
Some type information was truncated. Use `show(err)` to see complete types.

The issue arises because I'm trying to check if a matrix is symmetric or Hermitian using an atol:

julia> isapprox(B, B')
ERROR: ArgumentError: cannot set off-diagonal entry (2, 1)
Stacktrace:
 [1] setindex!
   @ ~/.asdf/installs/julia/1.11.0/share/julia/stdlib/v1.11/LinearAlgebra/src/tridiag.jl:466 [inlined]
...
Some type information was truncated. Use `show(err)` to see complete types.

This test is implemented in the IsApprox.jl package.

My initial suggestion is to check whether A and B have underlying symmetric parents and handle those cases explicitly, but I am uncertain if this is the best solution.

System info:

julia> versioninfo()
Julia Version 1.11.0
Commit 501a4f25c2b (2024-10-07 11:40 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 16 × Apple M3 Max
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, apple-m3)