Julia JuMP StackOverflow Error when defining a constraint on linear matrix inequalities

Question

I was inspired by this post to switch from MATLAB's LMI toolbox to a Julia implementation and I tried to adapt the suggested script in the post to a problem I am trying to solve. However, I have failed to do this, consistently causing Julia to crash or give me a StackOverflowError while trying to solve the following LMI problem with JuMP,

Lyapunov stability of hybrid systems

For the unknown symmetric matrices Ui, Wi, M (and Pi) where Ui, Wi are greater than or equal to zero.

I have run the model set-up line by line and found that it occurs when I define the last constraint:

@SDconstraint(model, Ai'*Pi + Pi*Ai + Ei'*Ui*Ei ⪯ 0)

I saw other posts about this with JuMP and they usually were a result of improper variable definition but that knowledge has not allowed me to figure out the issue. Here is a minimal form of the code that produces the error,

using JuMP
using MosekTools
using LinearAlgebra

E1p = [1 0; 0 1];
A1 = [-1 10; -100 -1];

Ei = E1p
Ai = A1
n=2

model = Model(Mosek.Optimizer)
@variable(model, Pi[i=1:n, j=1:n], Symmetric)
@variable(model, Ui[i=1:n, j=1:n], Symmetric)
@variable(model, Wi[i=1:n, j=1:n], Symmetric)
@variable(model, M[i=1:n, j=1:n], Symmetric)
@SDconstraint(model, Ui ⪰ 0)
@SDconstraint(model, Wi ⪰ 0)
@SDconstraint(model, Ei'*M*Ei ⪰ 0)
@SDconstraint(model, Pi - Ei'*Wi*Ei ⪰ 0)
@SDconstraint(model, Ai'*Pi + Pi*Ai + Ei'*Ui*Ei ⪯ 0)
optimize!(model)
Pi = value.(Pi)

The first line of the error which repeats to great length is:

in mutable_operate! at MutableArithmetics/bPWR4/src/linear_algebra.jl:132

Answer 1

That looks like a bug.

Here's a MWE

using JuMP
model = Model()
E = [1 1; 1 1]
@variable(model, P[1:size(E, 1), 1:size(E, 2)], Symmetric)
julia> @SDconstraint(model, E' * P + P * E + E' * P * E <= 0)
ERROR: StackOverflowError:
Stacktrace:
 [1] mutable_operate!(::typeof(MutableArithmetics.sub_mul), ::Matrix{AffExpr}, ::LinearAlgebra.Symmetric{VariableRef, Matrix{VariableRef}}, ::Matrix{Int64}) (repeats 79984 times)
   @ MutableArithmetics ~/.julia/packages/MutableArithmetics/bPWR4/src/linear_algebra.jl:132

If I drop the first term, I get

julia> @SDconstraint(model, P * E + E' * P * E <= 0)
ERROR: MethodError: no method matching *(::Matrix{AffExpr})
Closest candidates are:
  *(::StridedMatrix{T} where T, ::LinearAlgebra.AbstractQ) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/qr.jl:680
  *(::StridedMatrix{T} where T, ::LinearAlgebra.Adjoint{var"#s832", var"#s831"} where {var"#s832", var"#s831"<:LinearAlgebra.AbstractQ}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/qr.jl:720
  *(::StridedVecOrMat{T} where T, ::LinearAlgebra.Adjoint{var"#s832", var"#s831"} where {var"#s832", var"#s831"<:LinearAlgebra.LQPackedQ}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lq.jl:254
  ...
Stacktrace:
 [1] mutable_operate!(::typeof(MutableArithmetics.sub_mul), ::Matrix{AffExpr}, ::LinearAlgebra.Adjoint{Int64, Matrix{Int64}}, ::Matrix{AffExpr}) (repeats 2 times)
   @ MutableArithmetics ~/.julia/packages/MutableArithmetics/bPWR4/src/linear_algebra.jl:132
 [2] operate_fallback!(::MutableArithmetics.IsMutable, ::Function, ::Matrix{AffExpr}, ::LinearAlgebra.Adjoint{Int64, Matrix{Int64}}, ::LinearAlgebra.Symmetric{VariableRef, Matrix{VariableRef}}, ::Matrix{Int64})
   @ MutableArithmetics ~/.julia/packages/MutableArithmetics/bPWR4/src/interface.jl:330
 [3] operate!(::typeof(MutableArithmetics.sub_mul), ::Matrix{AffExpr}, ::LinearAlgebra.Adjoint{Int64, Matrix{Int64}}, ::LinearAlgebra.Symmetric{VariableRef, Matrix{VariableRef}}, ::Matrix{Int64})
   @ MutableArithmetics ~/.julia/packages/MutableArithmetics/bPWR4/src/rewrite.jl:80
 [4] macro expansion
   @ ~/.julia/packages/MutableArithmetics/bPWR4/src/rewrite.jl:276 [inlined]
 [5] macro expansion
   @ ~/.julia/packages/JuMP/y5vgk/src/macros.jl:447 [inlined]
 [6] top-level scope
   @ REPL[15]:1

So it's probably related to https://github.com/jump-dev/MutableArithmetics.jl/issues/84.

I'll open an issue. (Edit: https://github.com/jump-dev/MutableArithmetics.jl/issues/86)

As a work-around, split it into an expression first:

julia> @expression(model, ex, E' * P + P * E + E' * P * E)
2×2 Matrix{AffExpr}:
 3 P[1,1] + 4 P[1,2] + P[2,2]    4 P[1,2] + 2 P[2,2] + 2 P[1,1]
 2 P[1,1] + 4 P[1,2] + 2 P[2,2]  4 P[1,2] + 3 P[2,2] + P[1,1]

julia> @SDconstraint(model, ex <= 0)
[-3 P[1,1] - 4 P[1,2] - P[2,2]    -2 P[1,1] - 4 P[1,2] - 2 P[2,2];
 -2 P[1,1] - 4 P[1,2] - 2 P[2,2]  -P[1,1] - 4 P[1,2] - 3 P[2,2]] ∈ PSDCone()