From 4fd958233cecf66d267c45afadba52dd49afe80b Mon Sep 17 00:00:00 2001
From: Oscar Dowson <odow@users.noreply.github.com>
Date: Mon, 28 Aug 2023 10:56:26 +1200
Subject: [PATCH] Add upper_triangle and docs for
 {Log,Root}DetCone{Triangle,Square} (#3456)

---
 docs/src/tutorials/conic/ellipse_approx.jl    |  3 +-
 docs/src/tutorials/conic/experiment_design.jl |  5 +-
 docs/src/tutorials/conic/tips_and_tricks.jl   | 74 +++++++++++++++++++
 src/sd.jl                                     | 24 +++++-
 test/test_constraint.jl                       | 11 +++
 5 files changed, 110 insertions(+), 7 deletions(-)

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