diff --git a/docs/src/tutorials/conic/ellipse_approx.jl b/docs/src/tutorials/conic/ellipse_approx.jl index 81fbbf5082d..6b6f1a2f2ed 100644 --- a/docs/src/tutorials/conic/ellipse_approx.jl +++ b/docs/src/tutorials/conic/ellipse_approx.jl @@ -209,8 +209,7 @@ set_silent(model) f = [1 - S[i, :]' * Z * S[i, :] + 2 * S[i, :]' * z - s for i in 1:m] @constraint(model, f in MOI.Nonnegatives(m)) ## The former constraint [t; vec(Z)] in MOI.RootDetConeSquare(n) -Z_upper = [Z[i, j] for j in 1:n for i in 1:j] -@constraint(model, 1 * vcat(t, Z_upper) .+ 0 in MOI.RootDetConeTriangle(n)) +@constraint(model, 1 * [t; triangle_vec(Z)] .+ 0 in MOI.RootDetConeTriangle(n)) ## The former @objective(model, Max, t) @objective(model, Max, 1 * t + 0) optimize!(model) diff --git a/docs/src/tutorials/conic/experiment_design.jl b/docs/src/tutorials/conic/experiment_design.jl index 88edd583dc4..f8763b6b76f 100644 --- a/docs/src/tutorials/conic/experiment_design.jl +++ b/docs/src/tutorials/conic/experiment_design.jl @@ -211,10 +211,7 @@ set_silent(dOpt) @objective(dOpt, Max, t) @constraint(dOpt, sum(np) <= n) E = V * LinearAlgebra.diagm(0 => np ./ n) * V' -@constraint( - dOpt, - [t, 1, (E[i, j] for i in 1:q for j in 1:i)...] in MOI.LogDetConeTriangle(q) -) +@constraint(dOpt, [t; 1; triangle_vec(E)] in MOI.LogDetConeTriangle(q)) optimize!(dOpt) objective_value(dOpt) #- diff --git a/docs/src/tutorials/conic/tips_and_tricks.jl b/docs/src/tutorials/conic/tips_and_tricks.jl index 2471a19029b..6d558a2116c 100644 --- a/docs/src/tutorials/conic/tips_and_tricks.jl +++ b/docs/src/tutorials/conic/tips_and_tricks.jl @@ -338,6 +338,80 @@ set_silent(model) optimize!(model) value(t), value.(x) +# ## RootDetCone + +# The [`MOI.RootDetConeSquare`](@ref) is a cone of the form: +# ```math +# K = \{ (t, X) \in \mathbb{R}^{1+d^2} : t \le \det(X)^{\frac{1}{d}} \} +# ``` + +model = Model(SCS.Optimizer) +set_silent(model) +@variable(model, t) +@variable(model, X[1:2, 1:2]) +@objective(model, Max, t) +@constraint(model, [t; vec(X)] in MOI.RootDetConeSquare(2)) +@constraint(model, X .== [2 1; 1 3]) +optimize!(model) +value(t), sqrt(LinearAlgebra.det(value.(X))) + +# If `X` is symmetric, then you can use [`MOI.RootDetConeTriangle`](@ref) +# instead. This can be more efficient because the solver does not need to add +# additional constraints to ensure `X` is symmetric. + +# When forming the function, use [`triangle_vec`](@ref) to obtain the +# column-wise upper triangle of the matrix as a vector in the order that JuMP +# requires. + +model = Model(SCS.Optimizer) +set_silent(model) +@variable(model, t) +@variable(model, X[1:2, 1:2], Symmetric) +@objective(model, Max, t) +@constraint(model, [t; triangle_vec(X)] in MOI.RootDetConeTriangle(2)) +@constraint(model, X .== [2 1; 1 3]) +optimize!(model) +value(t), sqrt(LinearAlgebra.det(value.(X))) + +# ## LogDetCone + +# The [`MOI.LogDetConeSquare`](@ref) is a cone of the form: +# ```math +# K = \{ (t, u, X) \in \mathbb{R}^{2+d^2} : t \le u \log(\det(X / u)) \} +# ``` + +model = Model(SCS.Optimizer) +set_silent(model) +@variable(model, t) +@variable(model, u) +@variable(model, X[1:2, 1:2]) +@objective(model, Max, t) +@constraint(model, [t; u; vec(X)] in MOI.LogDetConeSquare(2)) +@constraint(model, X .== [2 1; 1 3]) +@constraint(model, u == 0.5) +optimize!(model) +value(t), 0.5 * log(LinearAlgebra.det(value.(X) ./ 0.5)) + +# If `X` is symmetric, then you can use [`MOI.LogDetConeTriangle`](@ref) +# instead. This can be more efficient because the solver does not need to add +# additional constraints to ensure `X` is symmetric. + +# When forming the function, use [`triangle_vec`](@ref) to obtain the +# column-wise upper triangle of the matrix as a vector in the order that JuMP +# requires. + +model = Model(SCS.Optimizer) +set_silent(model) +@variable(model, t) +@variable(model, u) +@variable(model, X[1:2, 1:2], Symmetric) +@objective(model, Max, t) +@constraint(model, [t; u; triangle_vec(X)] in MOI.LogDetConeTriangle(2)) +@constraint(model, X .== [2 1; 1 3]) +@constraint(model, u == 0.5) +optimize!(model) +value(t), 0.5 * log(LinearAlgebra.det(value.(X) ./ 0.5)) + # ## Other Cones and Functions # For other cones supported by JuMP, check out the diff --git a/src/sd.jl b/src/sd.jl index d409771d324..77388c7392c 100644 --- a/src/sd.jl +++ b/src/sd.jl @@ -171,11 +171,33 @@ function reshape_set( return PSDCone() end -function vectorize(matrix, ::SymmetricMatrixShape) +""" + triangle_vec(matrix::Matrix) + +Return the upper triangle of a matrix concatenated into a vector in the order +required by JuMP and MathOptInterface for `Triangle` sets. + +## Example + +```jldoctest +julia> model = Model(); + +julia> @variable(model, X[1:3, 1:3], Symmetric); + +julia> @variable(model, t) +t + +julia> @constraint(model, [t; triangle_vec(X)] in MOI.RootDetConeTriangle(3)) +[t, X[1,1], X[1,2], X[2,2], X[1,3], X[2,3], X[3,3]] ∈ MathOptInterface.RootDetConeTriangle(3) +``` +""" +function triangle_vec(matrix::AbstractMatrix) n = LinearAlgebra.checksquare(matrix) return [matrix[i, j] for j in 1:n for i in 1:j] end +vectorize(matrix, ::SymmetricMatrixShape) = triangle_vec(matrix) + """ SkewSymmetricMatrixShape diff --git a/test/test_constraint.jl b/test/test_constraint.jl index 16f256baade..eaed17cbb43 100644 --- a/test/test_constraint.jl +++ b/test/test_constraint.jl @@ -1710,4 +1710,15 @@ function test_eltype_for_constraint_primal_complex_float64() return end +function test_triangle_vec() + model = Model() + x = [1 2 3; 4 5 6; 7 8 9] + @test triangle_vec(x) == [1, 2, 5, 3, 6, 9] + @variable(model, y[1:2, 1:2]) + @test triangle_vec(y) == [y[1, 1], y[1, 2], y[2, 2]] + @variable(model, z[1:2, 1:2], Symmetric) + @test triangle_vec(z) == [z[1, 1], z[1, 2], z[2, 2]] + return +end + end