Second-order cone constraint in JuMP

Specific Domains Optimization (Mathematical)
, ,
1

Hi guys!
How to solve this constraint in Julia in JuMP, where the Q is a 2×2 Hessian matrix of a convex function, so it is a positive semidefinite matrix.

I know it can be reformulated as a second-order cone constraint, but I don’t know how to transfer.

Thank you very much!

image
image628×202 5.9 KB

2

If you are using Gurobi, you can just write it without reformulation.

If your Q is real:

  1. if it is positive definite, then it is easy—a cholesky decomposition is sufficient.
  2. if it is positive semidefinite, then maybe you can do a spectral decomposition

something like

# This `A` is a PSD matrix
A = [36 -26 -18 -24 -4; -26 59 -32 81 -15; -18 -32 62 -60 37; -24 81 -60 117 -30; -4 -15 37 -30 93]

using LinearAlgebra
vals, Q = eigen(A); # spectral decomposition
D = Diagonal(vals);
Q = Matrix(Q');
# @assert (Q'D)Q == A

# suppose you want to add this constr
t ≥ (x'A)x 
# you can define
JuMP.@expression(model, y, (Q)x)
# then the constr becomes
t ≥ (y'D)y
# since `D` is Diagonal, this constr can be converted to SOC form

Then after removing 0 entries in D, then followed by a cholesky decomposition, you can write the last constraint as a rotated second order cone here.

Then you can in turn transform that rotated second order cone to a standard SOC.

PS to know whether a given symmetric matrix is positive definite, you can do cholesky directly, if success, then it is PD, otherwise you can do spectral decomposition as a recourse.

3

Hi @Z_M,

Here’s the Cholesky version:

using JuMP
Q = [1.0 0.5; 0.5 3.0]
model = Model()
@variable(model, t)
@variable(model, x[i in 1:2])
@objective(model, Min, t)
QQ = LinearAlgebra.cholesky(Q)
# 2 * t * 1 >= || QQ.U * x ||_2^2
@constraint(model, [t; 1; QQ.U * x] in RotatedSecondOrderCone())
4

Thank you!
I use Gurobi in JuMP, so how to write it without reformulation.
In addition, if I have the JuMP.@expression(model, y, (Q)x) and then how to write a soc constraint in JuMP.

Thank you very much!

5

Thank you for your reply.
When I run this code, for the constraint: @constraint(model, [t; 1; QQ.U * x] in RotatedSecondOrderCone())

It always has an error:
MethodError: Cannot convert an object of type QuadExpr to an object of type AffExpr

In addition, I use Gurobi in JuMP, so I wonder if Gurobi can solve the RotatedSecondOrderCone constraint.

6

Gurobi supports quadratic constraints. Just write it as:

using JuMP
Q = [1.0 0.5; 0.5 3.0]
model = Model()
@variable(model, t)
@variable(model, x[i in 1:2])
@objective(model, Min, t)
@constraint(model, t >= 0.5 * x' * Q * x)

The second order cone equivalent is:

@constraint(model, [(t + 1) / sqrt(2); (t - 1) / sqrt(2); QQ.U * x] in SecondOrderCone())

(If it isn’t obvious why, try expanding out the SOC constraint to the quadratic form.)

7

This error is not coming from the code I wrote. Can you provide a reproducible example?

8

I run the following code and get the error.

Q = [1.0 0.5; 0.5 3.0]

QQ = LinearAlgebra.cholesky(Q)

model = Model(Gurobi.Optimizer)

@variable(model, t)

@variable(model, x[i in 1:2])

@objective(model, Min, t)

@constraint(model, t >= 0.5 * x’ * Q * x)

@constraint(model, [(t + 1) / sqrt(2); (t - 1) / sqrt(2); QQ.U * x] in SecondOrderCone())

optimize!(model)

image
image1728×969 97.3 KB

9

Thank you very much! I got it.

10

I can’t reproduce:

julia> using JuMP, Gurobi, LinearAlgebra

julia> Q = [1.0 0.5; 0.5 3.0]
2×2 Matrix{Float64}:
 1.0  0.5
 0.5  3.0

julia> model = Model(Gurobi.Optimizer)
Set parameter LicenseID to value 890341
A JuMP Model
├ solver: Gurobi
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none

julia> @variable(model, t)
t

julia> @variable(model, x[i in 1:2])
2-element Vector{VariableRef}:
 x[1]
 x[2]

julia> @objective(model, Min, t)
t

julia> QQ = LinearAlgebra.cholesky(Q)
Cholesky{Float64, Matrix{Float64}}
U factor:
2×2 UpperTriangular{Float64, Matrix{Float64}}:
 1.0  0.5
  ⋅   1.65831

julia> @constraint(model, [(t + 1) / sqrt(2); (t - 1) / sqrt(2); QQ.U * x] in SecondOrderCone())
[0.7071067811865475 t + 0.7071067811865475, 0.7071067811865475 t - 0.7071067811865475, x[1] + 0.5 x[2], 1.6583123951777 x[2]] ∈ MathOptInterface.SecondOrderCone(4)

julia> optimize!(model)
Gurobi Optimizer version 12.0.2 build v12.0.2rc0 (mac64[x86] - Darwin 24.1.0 24B83)

CPU model: Intel(R) Core(TM) i5-8259U CPU @ 2.30GHz
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 4 rows, 7 columns and 9 nonzeros
Model fingerprint: 0x303d1a5d
Model has 1 quadratic constraint
Coefficient statistics:
  Matrix range     [5e-01, 2e+00]
  QMatrix range    [1e+00, 1e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [7e-01, 7e-01]
Presolve removed 3 rows and 5 columns
Presolve time: 0.00s
Presolved: 1 rows, 2 columns, 2 nonzeros
Presolved model has 1 second-order cone constraint
Ordering time: 0.00s

Barrier statistics:
 AA' NZ     : 0.000e+00
 Factor NZ  : 1.000e+00
 Factor Ops : 1.000e+00 (less than 1 second per iteration)
 Threads    : 1

                  Objective                Residual
Iter       Primal          Dual         Primal    Dual     Compl     Time
   0   1.20979837e-02  0.00000000e+00  8.55e-03 7.04e-03  5.54e-03     0s
   1   1.31297546e-05 -1.79908716e-07  9.18e-06 7.46e-06  6.05e-06     0s
   2   1.44738910e-08 -1.63817071e-10  1.01e-08 8.26e-09  6.66e-09     0s

Barrier solved model in 2 iterations and 0.00 seconds (0.00 work units)
Optimal objective 1.44738910e-08


User-callback calls 62, time in user-callback 0.00 sec