jump-dev / MutableArithmetics.jl

Interface for arithmetics on mutable types in Julia
Mozilla Public License 2.0
50 stars 8 forks source link

Errors in complex algebra with Symmetric and Hermitian #264

Closed araujoms closed 8 months ago

araujoms commented 8 months ago

The algebra with complex numbers and Symmetric or Hermitian matrices doesn't work as it should. First of all, multiplying a complex number by a Hermitian matrix gives an incorrect result wrapped in Hermitian:

using JuMP
julia> model = Model();
julia> @variable(model, m[1:2,1:2] in HermitianPSDCone());
julia> im*m
2×2 LinearAlgebra.Hermitian{GenericAffExpr{ComplexF64, VariableRef}, Matrix{GenericAffExpr{ComplexF64, VariableRef}}}:
 -real(m[1,1]) im                 real(m[1,2]) im - imag(m[1,2])
 -real(m[1,2]) im - imag(m[1,2])  -real(m[2,2]) im

Also, if I create a real variable and multiply it by a Symmetric or Hermitian matrix, I get the correct result, but it's no longer wrapped in Symmetric or Hermitian:

using JuMP
julia> model = Model();
julia> @variable(model, x);
julia> m = randn(ComplexF64,2,2);
julia> x*Symmetric(m)
2×2 Matrix{GenericAffExpr{ComplexF64, VariableRef}}:
 (-1.2483493211093275 + 0.5094508206087239im) x  (-0.6295524575806658 + 1.0006355385104637im) x
 (-0.6295524575806658 + 1.0006355385104637im) x  (-0.8701716643747957 - 0.27288163539446847im) x
julia> x*Hermitian(m)
2×2 Matrix{GenericAffExpr{ComplexF64, VariableRef}}:
 -1.2483493211093275 x                           (-0.6295524575806658 + 1.0006355385104637im) x
 (-0.6295524575806658 - 1.0006355385104637im) x  -0.8701716643747957 x

(Of course, Symmetric matrices multiplied by complex variables should still be Symmetric, but I don't think this will ever matter.)

blegat commented 8 months ago

The issue is that we shouldn't wrap in Hermitian something that's not Hermitian ^^

julia> B = MA._mult_upper(im, m)
2×2 Matrix{GenericAffExpr{ComplexF64, VariableRef}}:
 real(m[1,1]) im                 real(m[1,2]) im - imag(m[1,2])
 real(m[1,2]) im + imag(m[1,2])  real(m[2,2]) im

julia> h = Hermitian(B, :U)
2×2 Hermitian{GenericAffExpr{ComplexF64, VariableRef}, Matrix{GenericAffExpr{ComplexF64, VariableRef}}}:
 -real(m[1,1]) im                 real(m[1,2]) im - imag(m[1,2])
 -real(m[1,2]) im - imag(m[1,2])  -real(m[2,2]) im

julia> h.data
2×2 Matrix{GenericAffExpr{ComplexF64, VariableRef}}:
 real(m[1,1]) im                 real(m[1,2]) im - imag(m[1,2])
 real(m[1,2]) im + imag(m[1,2])  real(m[2,2]) im
araujoms commented 8 months ago

That's weird. Shouldn't Hermitian(B) return instead

2×2 Hermitian{GenericAffExpr{ComplexF64, VariableRef}, Matrix{GenericAffExpr{ComplexF64, VariableRef}}}:
 0.0                              real(m[1,2]) im - imag(m[1,2])
 -real(m[1,2]) im - imag(m[1,2])  0.0
blegat commented 8 months ago

Yes, looks like another bug, it calls LinearAlgebra.hermitian(h.data[1, 1], :U) which doesn't seem to do the right thing see https://github.com/jump-dev/JuMP.jl/issues/3692

odow commented 8 months ago

With #265, we now get:

julia> using JuMP

julia> model = Model();

julia> @variable(model, m[1:2,1:2] in HermitianPSDCone());

julia> im*m
2×2 Matrix{GenericAffExpr{ComplexF64, VariableRef}}:
 real(m[1,1]) im                 real(m[1,2]) im - imag(m[1,2])
 real(m[1,2]) im + imag(m[1,2])  real(m[2,2]) im
araujoms commented 8 months ago

Wait, please, don't close the issue! What about getting a new dispatch so that a Hermitian matrix multiplied by a real variable stays Hermitian? It is the only reason I was even testing this stuff.

odow commented 8 months ago

What about getting a new dispatch so that a Hermitian matrix multiplied by a real variable stays Hermitian?

I'll open a separate issue.

araujoms commented 8 months ago

Great, thanks.