Stack overflow error with @SDconstraint

Specific Domains Optimization (Mathematical)
1

I have the following simple code snippet which is causing me issues:

function find_common_lyapunov(system::PiecewiseLinearSystem, gamma::Float64)
    @assert(gamma > 0.0)
    dim = system_dim(system)
    model = Model(Mosek.Optimizer)
    @variable(model, P[1:dim, 1:dim], PSD)
    @SDconstraint(model, P >= I)

    @show(length(system.matrices))
    for Ai in system.matrices
        @show(Ai)
        @SDconstraint(model, Ai' * P + P * Ai + gamma * P <= 0)
    end 
                                                                                                                                       
    @objective(model, Min, sum(P .^ 2)) 
                                                                                                                                       
    optimize!(model)
    status = termination_status(model)
    if status != MOI.OPTIMAL
        println("TERMINATION WAS NOT OPTIMAL: $status")
        return nothing
    end 
                                                                                                                                       
    value.(P)
end

When I run this, I get the following:

julia> P_star = find_common_lyapunov(system, 1e-3)
length(system.matrices) = 4
Ai = [-0.3216088732749837 0.2456225325321437 0.0 0.0 0.0 0.0 0.0 0.0; -4.550258232817038 -1.0264957875870517 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 -0.01 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 -0.01 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 -0.01 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 -0.01 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 -0.01 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.01]
ERROR: StackOverflowError:
Stacktrace:
 [1] mutable_operate!(::typeof(MutableArithmetics.sub_mul), ::Array{GenericAffExpr{Float64,VariableRef},2}, ::Symmetric{VariableRef,Array{VariableRef,2}}, ::Array{Float64,2}) at /Users/stephentu/.julia/packages/MutableArithmetics/ZGFsK/src/linear_algebra.jl:69 (repeats 79984 times)

This seems like too small of an SDP to run into memory issues. Any ideas?

2

This looks like a bug. Here is a simpler reproducer:

using JuMP
model = Model()
@variable(model, P[1:2, 1:2], PSD)
A = [1 0; 0 1]
@SDconstraint(model, A' * P .+ P * A <= 0)  # Works
@SDconstraint(model, A' * P + P * A <= 0)  # Bus 10 error

I’ll open an issue. Already reported: https://github.com/jump-dev/MutableArithmetics.jl/issues/48

1 Like
3

Thanks @odow. Does the .+ line you showed have the same semantics as the + (can I use it as a proper workaround for now?)

4

It seemed to work for me, but you’ll have to try it out to make sure.

In this case, .+ is equivalent to +. It is just native Julia broadcasting: Mathematical Operations and Elementary Functions · The Julia Language

(So note that P .* A is _not) the same as P * A.)