Symmetric(StaticMatrix) causes several problems

[vincent-picaud]

Hello,

Thank you for this package,

Starting to use it, I encountered several problems with Symmetric(StaticMatrix). This prevents me to have some generic code (working without modification for StaticMatrices and Julia regular matrix types).

Example 1: M\v strange behavior

M\v does not work with this invertible M matrix

julia> using LinearAlgebra, StaticArrays

julia> M = Symmetric(@SMatrix Float64[1 1;1 2])
2×2 Symmetric{Float64, SMatrix{2, 2, Float64, 4}}:
 1.0  1.0
 1.0  2.0

julia> v = @SVector Float64[1,2]
2-element SVector{2, Float64} with indices SOneTo(2):
 1.0
 2.0

julia> M\v
ERROR: setindex!(::SMatrix{2, 2, Float64, 4}, value, ::Int) is not defined.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:33
  [2] setindex!(a::SMatrix{2, 2, Float64, 4}, value::Float64, i::Int64)
    @ StaticArrays ~/.julia/packages/StaticArrays/OWJK7/src/indexing.jl:3
  [3] macro expansion
    @ ~/.julia/packages/StaticArrays/OWJK7/src/indexing.jl:66 [inlined]
  [4] _setindex!_scalar
    @ ~/.julia/packages/StaticArrays/OWJK7/src/indexing.jl:46 [inlined]
  [5] setindex!
    @ ~/.julia/packages/StaticArrays/OWJK7/src/indexing.jl:42 [inlined]
  [6] copytri!
    @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:514 [inlined]
  [7] copytri!
    @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:510 [inlined]
  [8] lu!(A::Symmetric{Float64, SMatrix{2, 2, Float64, 4}}, pivot::Val{true}; check::Bool)
    @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:89
  [9] lu(A::Symmetric{Float64, SMatrix{2, 2, Float64, 4}}, pivot::Val{true}; check::Bool)
    @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:273
 [10] lu(A::Symmetric{Float64, SMatrix{2, 2, Float64, 4}}, pivot::Val{true}) (repeats 2 times)
    @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:272
 [11] \(A::Symmetric{Float64, SMatrix{2, 2, Float64, 4}}, B::SVector{2, Float64})
    @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/generic.jl:1136
 [12] top-level scope
    @ REPL[119]:1

however, the same procedure works with the same matrix type but different component values:

julia> M = Symmetric(@SMatrix Float64[1 0;0 2])
2×2 Symmetric{Float64, SMatrix{2, 2, Float64, 4}}:
 1.0  0.0
 0.0  2.0

julia> M\v
2-element SVector{2, Float64} with indices SOneTo(2):
 1.0
 1.0

nb: same observation if one uses MMatrix instead of SMatrix.

Example 2: M + I

julia> using LinearAlgebra, StaticArrays

julia> M = Symmetric(@SMatrix Float64[1 0;0 2])
2×2 Symmetric{Float64, SMatrix{2, 2, Float64, 4}}:
 1.0  0.0
 0.0  2.0

julia> M + I
ERROR: setindex!(::SMatrix{2, 2, Float64, 4}, value, ::Int) is not defined.
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:33
 [2] setindex!(a::SMatrix{2, 2, Float64, 4}, value::Float64, i::Int64)
   @ StaticArrays ~/.julia/packages/StaticArrays/OWJK7/src/indexing.jl:3
 [3] macro expansion
   @ ~/.julia/packages/StaticArrays/OWJK7/src/indexing.jl:66 [inlined]
 [4] _setindex!_scalar
   @ ~/.julia/packages/StaticArrays/OWJK7/src/indexing.jl:46 [inlined]
 [5] setindex!
   @ ~/.julia/packages/StaticArrays/OWJK7/src/indexing.jl:42 [inlined]
 [6] setindex!(A::Symmetric{Float64, SMatrix{2, 2, Float64, 4}}, v::Float64, i::Int64, j::Int64)
   @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/symmetric.jl:226
 [7] +(A::Symmetric{Float64, SMatrix{2, 2, Float64, 4}}, J::UniformScaling{Bool})
   @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/uniformscaling.jl:219
 [8] top-level scope
   @ REPL[132]:1

Looking at the code I understand the origin of the problem:

function (+)(A::AbstractMatrix, J::UniformScaling)
    checksquare(A)
    B = copy_oftype(A, Base._return_type(+, Tuple{eltype(A), typeof(J)}))
    @inbounds for i in axes(A, 1)
        B[i,i] += J
    end
    return B
end

The function copy_oftype() create a copy of of M.data of type SMatrix then try to modify its diagonal.

As expected, this does not cause an issue with MMatrix:

Julia> M = Symmetric(@MMatrix Float64[1 1;1 2])
2×2 Symmetric{Float64, MMatrix{2, 2, Float64, 4}}:
 1.0  1.0
 1.0  2.0

julia> M + I
2×2 Symmetric{Float64, MMatrix{2, 2, Float64, 4}}:
 2.0  1.0
 1.0  3.0

Sorry to bring this to the table, Thanks. Vincent

[mateuszbaran]

I think this shouldn't be hard to fix. Implementations in LinearAlgebra usually assume that array types are mutable and thus MArray often just works while SArray needs custom methods.

[vincent-picaud]

Thank you for this quick answer.

I am currently using this as patch:

using LinearAlgebra: Symmetric, UniformScaling 

import LinearAlgebra

function LinearAlgebra. +(A::Symmetric{T,SMatrix{N,N,T,L}},J::UniformScaling) where {T,N,L}
    Symmetric(A.data + J,A.uplo)
end

function LinearAlgebra. \(A::Union{Symmetric{T,SMatrix{N,N,T,L}},Symmetric{T,MMatrix{N,N,T,L}}}, B::AbstractVecOrMat) where {T,N,L}
    M = convert(SMatrix{N,N,T,L},A)

    M\B
end

However, I still do not understand the origin of the problem in "Example 1". I am surprised that the behavior depends on values (with the same matrix type). Maybe a different code path due to pivot selection? Also please note that in this "Example 1", we get the same exact behavior with MMatrix instead of SMatrix.

[mateuszbaran]

Yes, in example 1 a different code path in LinearAlgebra is taken, one that doesn't do any mutation. Note that the generic M\v that is used checks first whether M is triangular; your first matrix isn't but the second one is.

[vincent-picaud]

Ok, thanks. I get it now... this "triangular matrix" is detected at runtime... I'm still too infused by compiled languages like C++ for this to come to mind :)

---

I understand better the problem now. No problem for me if you want to close the issue

[mateuszbaran]

I think it's a good issue, and worth fixing in StaticArrays.jl.

[vincent-picaud]

Ok thank you for your quick response and interest. I am not sure that my quick patch is the right way to go otherwise I would have proposed a more completed PR.

[mateuszbaran]

The "right" way to fix this is a bit more complicated. I'll patch a part of this myself. Do you actually need solving M\v for such small matrices? If so, you could write dedicated solves for small sizes like here: https://github.com/JuliaArrays/StaticArrays.jl/blob/master/src/solve.jl .

[vincent-picaud]

I am developing a nonlinear squares solver, but please do not use it now... this is a WIP project :)

Usual problems will use Julia standard matrices types. However, my idea was to support StaticArrays for simple test functions like Rosenbrock, where I return approximate Hessian as a Symmetric(StaticMatrix) of size 2x2. Ideally, I want my "generic" solvers to support these two types of matrices without modification. And It is not important for me to have optimized M\v for tiny matrices, I only want to use the same code for both types :)

(my WIP code, with the previous path is there: https://github.com/vincent-picaud/NLS_Solver.jl/tree/no_more_inplace)

[mcabbott]

BTW, M + I is partially fixed by writing Base.axes(x::Symmetric) = axes(parent(x)), as this causes it to make an MMatrix to write into. Probably that should be done in any case, for all the wrappers.

[3f6a]

Another example, inv(::Symmetric) is not Symmetric for StaticArray:

X = @SMatrix randn(3,3) # StaticArray
X = Symmetric(X + X') / 2 # this Symmetric{..., SMatrix,....}
inv(X) # returns plain Matrix, but should be Symmetric!

Here inv(X) should be Symmetric.