Incorrect summation with Hermitian matrices

[odow]

This took me a while to track down!!!

julia> A = LinearAlgebra.Hermitian([1.0 2.0; 2.0 3.0] .+ [0 -4.0; 4.0 0.0] .* im)
2×2 LinearAlgebra.Hermitian{ComplexF64, Matrix{ComplexF64}}:
 1.0+0.0im  2.0-4.0im
 2.0+4.0im  3.0+0.0im

julia> model = Model();

julia> @variable(model, X[1:2, 1:2] in HermitianPSDCone())
2×2 Matrix{GenericAffExpr{ComplexF64, VariableRef}}:
 real(X[1,1])                                real(X[1,2]) + (0.0 + 1.0im) imag(X[1,2])
 real(X[1,2]) + (-0.0 - 1.0im) imag(X[1,2])  real(X[2,2])

julia> Y = LinearAlgebra.I - X
2×2 Matrix{GenericAffExpr{ComplexF64, VariableRef}}:
 (-1.0 - 0.0im) real(X[1,1]) + (1.0 - 0.0im)               (-1.0 - 0.0im) real(X[1,2]) + (-0.0 - 1.0im) imag(X[1,2])
 (-1.0 - 0.0im) real(X[1,2]) + (0.0 + 1.0im) imag(X[1,2])  (-1.0 - 0.0im) real(X[2,2]) + (1.0 - 0.0im)

julia> Z = LinearAlgebra.Hermitian(Y)
2×2 LinearAlgebra.Hermitian{GenericAffExpr{ComplexF64, VariableRef}, Matrix{GenericAffExpr{ComplexF64, VariableRef}}}:
 (-1.0 + 0.0im) real(X[1,1]) + (1.0 + 0.0im)                (-1.0 - 0.0im) real(X[1,2]) + (-0.0 - 1.0im) imag(X[1,2])
 (-1.0 + 0.0im) real(X[1,2]) + (-0.0 + 1.0im) imag(X[1,2])  (-1.0 + 0.0im) real(X[2,2]) + (1.0 + 0.0im)

julia> Y == Z
true

julia> real(LinearAlgebra.dot(A, X) + LinearAlgebra.dot(A, Y))
4

julia> real(LinearAlgebra.dot(A, X) + LinearAlgebra.dot(A, Z))
-16 imag(X[1,2]) + 4

[odow]

julia> A = [0 1; 1 0] .+ [0 -1; 1 0] .* im
2×2 Matrix{Complex{Int64}}:
 0+0im  1-1im
 1+1im  0+0im

julia> H = LinearAlgebra.Hermitian(A)
2×2 LinearAlgebra.Hermitian{Complex{Int64}, Matrix{Complex{Int64}}}:
 0+0im  1-1im
 1+1im  0+0im

julia> model = Model();

julia> @variable(model, X[1:2, 1:2] in HermitianPSDCone())
2×2 Matrix{GenericAffExpr{ComplexF64, VariableRef}}:
 real(X[1,1])                                real(X[1,2]) + (0.0 + 1.0im) imag(X[1,2])
 real(X[1,2]) + (-0.0 - 1.0im) imag(X[1,2])  real(X[2,2])

julia> Y = LinearAlgebra.Hermitian(X)
2×2 LinearAlgebra.Hermitian{GenericAffExpr{ComplexF64, VariableRef}, Matrix{GenericAffExpr{ComplexF64, VariableRef}}}:
 real(X[1,1])                               real(X[1,2]) + (0.0 + 1.0im) imag(X[1,2])
 real(X[1,2]) + (0.0 - 1.0im) imag(X[1,2])  real(X[2,2])

julia> LinearAlgebra.dot(A, X)
(2.0 + 0.0im) real(X[1,2]) + (-2.0 + 0.0im) imag(X[1,2])

julia> LinearAlgebra.dot(A, Y)
(2.0 + 0.0im) real(X[1,2]) + (-2.0 + 0.0im) imag(X[1,2])

julia> LinearAlgebra.dot(H, X)
(2.0 + 0.0im) real(X[1,2]) + (-2.0 + 0.0im) imag(X[1,2])

julia> LinearAlgebra.dot(H, Y)
(2.0 + 0.0im) real(X[1,2]) + (2.0 + 0.0im) imag(X[1,2])

Note the sign change

[odow]

Looks ike this might be a problem in MutableArithmetics

julia> @which LinearAlgebra.dot(A, X)
dot(lhs::AbstractArray, rhs::AbstractArray{var"#s36", N} where {var"#s36"<:MutableArithmetics.AbstractMutable, N}) in MutableArithmetics at /Users/oscar/.julia/packages/MutableArithmetics/geMUn/src/dispatch.jl:61

julia> @which LinearAlgebra.dot(A, Y)
dot(lhs::AbstractArray, rhs::AbstractArray{var"#s36", N} where {var"#s36"<:MutableArithmetics.AbstractMutable, N}) in MutableArithmetics at /Users/oscar/.julia/packages/MutableArithmetics/geMUn/src/dispatch.jl:61

julia> @which LinearAlgebra.dot(H, X)
dot(lhs::AbstractArray, rhs::AbstractArray{var"#s36", N} where {var"#s36"<:MutableArithmetics.AbstractMutable, N}) in MutableArithmetics at /Users/oscar/.julia/packages/MutableArithmetics/geMUn/src/dispatch.jl:61

julia> @which LinearAlgebra.dot(H, Y)
dot(A::LinearAlgebra.Hermitian, B::LinearAlgebra.Hermitian) in LinearAlgebra at /Users/oscar/.julia/juliaup/julia-1.6.7+0.x64.apple.darwin14/share/julia/stdlib/v1.6/LinearAlgebra/src/symmetric.jl:422

[odow]

This is a bug in MA:

julia> using JuMP

julia> const MA = JuMP._MA
MutableArithmetics

julia> A = [0 1; 1 0] .+ [0 -1; 1 0] .* im
2×2 Matrix{Complex{Int64}}:
 0+0im  1-1im
 1+1im  0+0im

julia> model = Model();

julia> @variable(model, X[1:2, 1:2] in HermitianPSDCone())
2×2 Matrix{GenericAffExpr{ComplexF64, VariableRef}}:
 real(X[1,1])                                real(X[1,2]) + (0.0 + 1.0im) imag(X[1,2])
 real(X[1,2]) + (-0.0 - 1.0im) imag(X[1,2])  real(X[2,2])

julia> y = zero(GenericAffExpr{ComplexF64, VariableRef})
(0.0 + 0.0im)

julia> MA.operate!(MA.add_dot, y, A[2], adjoint(X[2]))
(1.0 - 1.0im) real(X[1,2]) + (1.0 + 1.0im) imag(X[1,2])

julia> y = zero(GenericAffExpr{ComplexF64, VariableRef})
(0.0 + 0.0im)

julia> MA.add_dot(y, A[2], X[2])
(1.0 + 1.0im) real(X[1,2]) + (1.0 - 1.0im) imag(X[1,2])

[blegat]

If I'm not mistaken,

julia> LinearAlgebra.dot(H, Y)
(2.0 + 0.0im) real(X[1,2]) + (2.0 + 0.0im) imag(X[1,2])

is the incorrect output while this is the one using LinearAlgebra

julia> @which LinearAlgebra.dot(H, Y)
dot(A::LinearAlgebra.Hermitian, B::LinearAlgebra.Hermitian) in LinearAlgebra at /Users/oscar/.julia/juliaup/julia-1.6.7+0.x64.apple.darwin14/share/julia/stdlib/v1.6/LinearAlgebra/src/symmetric.jl:422

This is because LinearAlgebra.dot is incorrect with JuMP types (fixed in #3244) and MA does not use LinearAlgebra.dot.

[odow]

D'oh. I just assumed the Base was correct. No wonder I was getting confused!