Cannot convert QuadExpr - error in Julia v1.7 but not v.1.8 with PSDCone()

Specific Domains Optimization (Mathematical)
, ,
1

Hi all,

I teach a course on convex optimization and use Julia/JuMP (along with Matlab/CVX) to introduce students to implementation aspects – however, I am not a super-user of Julia/JuMP and may be missing something obvious, but I am a greatly puzzled and hope someone can help out.

Problem: I have a student implementing an SDP, who keeps getting the same error (on a Windows 10 laptop with Julia/JuMP v1.7/1.3) for code that runs fine on my computer (MacOS Julia/JuMP v1.8/1.3).

The error on the student’s end is attached in the figure: MethodError → QuadExpr / AffExpr about object type issues. The simple piece of code/script that works on my end, but fails on theirs is attached here.

If someone can help us understand why this issue happens on identical code for different computers that would be great, so I can avoid future issues with students. I have asked the student to update to Julia 1.8 (from 1.7) but am not sure that is the core issue and if so, I’d love to better understand what is happening.

Much appreciated and thank you for any responses!

using Graphs, LinearAlgebra, JuMP
using SCS


n=3;

A1=rand([0,1],n,n);  #random nxn adjaceny matrix. Values of 0 or 1.
A2=A1+A1';             #makes Adjacency matrix symmetric, implying undirected graph.
A=A2-Diagonal(diag(A2));      #makes diagonal of adjacency matrix zeros, removing self loops

model=Model(SCS.Optimizer);
# model=Model(SCS.Optimizer);

gamma=@variable(model)  
w=@variable(model, [1:n])

#constraints for fast mixing markov chain

@constraint(model, c1,  A*Diagonal(w)*A' + (1/n).*ones(n,n) <= gamma.*I(n) + I(n), PSDCone())
@constraint(model, c2,  I(n) - A*Diagonal(w)*A' - (1/n)*ones(n,n) <= gamma.*I(n), PSDCone())
@constraint(model, c3, A*Diagonal(w)*A' <= I(n), PSDCone())
@constraint(model, c4, w .>= 0)
@constraint(model, c5, gamma >= 0)

@objective(model,Min, gamma)

optimize!(model)

wVal = value.(w)
gammaVal = value(gamma)

MicrosoftTeams-image (17)
MicrosoftTeams-image (17)1420×540 42.5 KB

2

The issue is: LinearAlgebra matrix types and expressions · Issue #66 · jump-dev/MutableArithmetics.jl · GitHub

It’s hard to write methods that work when multiplying matrices with LinearAlgebra matrices because LinearAlgebra makes a number of assumptions that don’t hold in general. In particular, that the type of a result is stable under multiplication (that is, if x and y are both type T, then x * y is also type T).

Things improved a little in Julia v1.8, but issues persist in older Julia versions.

It’s a bit messy, but the solution is to wrap Matrix(...) around the LinearAlgebra type to convert it into a regular Julia Matrix type.

Here’s how I would write your model:

using JuMP
import LinearAlgebra
import SCS
n = 3
A1 = rand([0, 1], n, n)
A2 = A1 + A1'
A = A2 - LinearAlgebra.Diagonal(LinearAlgebra.diag(A2))
model = Model(SCS.Optimizer)
@variables(model, begin
    gamma >= 0
    w[1:n] >= 0
end)
@expressions(model, begin
    AwA, A * Matrix(LinearAlgebra.Diagonal(w)) * A'
end)
I = Matrix(LinearAlgebra.I(n))
@constraints(model, begin
    AwA .+ 1/n <= gamma * I + I, PSDCone()
    I - AwA .- 1/n <= gamma * I, PSDCone()
    AwA <= I, PSDCone()
end)
@objective(model, Min, gamma)
optimize!(model)
wVal, gammaVal = value.(w), value(gamma)
1 Like