From bfaae9211aba8034c6f1e0fde2e7c876fa39d724 Mon Sep 17 00:00:00 2001 From: odow Date: Mon, 5 Jun 2023 16:46:17 +1200 Subject: [PATCH 01/13] [docs] add Dualization tutorial --- docs/Project.toml | 2 + docs/make.jl | 1 + docs/src/tutorials/conic/dualization.jl | 210 ++++++++++++++++++++++++ 3 files changed, 213 insertions(+) create mode 100644 docs/src/tutorials/conic/dualization.jl diff --git a/docs/Project.toml b/docs/Project.toml index d089987599b..58455461b66 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -4,6 +4,7 @@ Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +Dualization = "191a621a-6537-11e9-281d-650236a99e60" GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6" HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" @@ -30,6 +31,7 @@ CSV = "0.10" Clarabel = "0.5" DataFrames = "1" Documenter = "0.27.9, 0.28" +Dualization = "0.5" GLPK = "=1.1.2" HTTP = "1.5.4" HiGHS = "=1.5.1" diff --git a/docs/make.jl b/docs/make.jl index 4c75d008cf9..d7916352d82 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -254,6 +254,7 @@ const _PAGES = [ "tutorials/conic/start_values.md", "tutorials/conic/tips_and_tricks.md", "tutorials/conic/simple_examples.md", + "tutorials/conic/dualization.md", "tutorials/conic/logistic_regression.md", "tutorials/conic/experiment_design.md", "tutorials/conic/min_ellipse.md", diff --git a/docs/src/tutorials/conic/dualization.jl b/docs/src/tutorials/conic/dualization.jl new file mode 100644 index 00000000000..4323e16e420 --- /dev/null +++ b/docs/src/tutorials/conic/dualization.jl @@ -0,0 +1,210 @@ +# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors #src +# This Source Code Form is subject to the terms of the Mozilla Public License #src +# v.2.0. If a copy of the MPL was not distributed with this file, You can #src +# obtain one at https://mozilla.org/MPL/2.0/. #src + +# # Dualization + +# The purpose of this tutorial is to explain how to use [Dualization.j](@ref) to +# improve the performance of some conic optimization models. There are two +# important takeaways: +# +# 1. JuMP reformulates problems to meet the input requirements of the +# solver, potentially increasing the problem size by adding slack variables +# and constraints. +# 2. Solving the dual of a conic model can be more efficient than solving the +# primal. + +# This tutorial uses the following packges + +using JuMP +import Dualization +import LinearAlgebra +import SCS + +# ## Background + +# Conic optimization solvers typically accept one of two input formulations: + +# The first is the _standard_ conic form: +# ```math +# \begin{align} +# \min_{x \in \mathbb{R}^n} & c^\top x \\ +# \;\;\text{s.t.} & A x = b \\ +# & x \in \mathcal{K} +# \end{align} +# ``` +# in which we have a set of linear equality constraints $Ax = b$ and the +# variables belong to a cone $\mathcal{K}$. +# +# The second is the _geometric_ conic form: +# ```math +# \begin{align} +# \min_{x \in \mathbb{R}^n} & c^\top x \\ +# \;\;\text{s.t.} & A x - b \in \mathcal{K} +# \end{align} +# ``` +# in which an affine function $Ax - b$ belongs to a cone $\mathcal{K}$ and the +# variables are free. + +# It is trivial to convert between these two representations, for example, to go +# from the geometric conic form to the standard conic form we introduce slack +# variables $y$: +# ```math +# \begin{align} +# \min_{x \in \mathbb{R}^n} & c^\top x \\ +# \;\;\text{s.t.} & [A -I] [x; y] = b \\ +# & [x; y] \in \mathbb{R}^n \cross \mathcal{K} +# \end{align} +# ``` +# and to go from the standard conic form to the geometric conic form, we can +# rewrite the equality constraint as a function belonging to the `{0}` cone: +# ```math +# \begin{align} +# \min_{x \in \mathbb{R}^n} & c^\top x \\ +# \;\;\text{s.t.} & [A; I] x - [b; 0] \in {0}\cross \mathcal{K} +# \end{align} +# ``` + +# From a theoretical perspective, the two formulations are equivalent, and if +# you implement a model in the standard conic form and pass it to a geometric +# conic form solver (or vice versa), then JuMP will automatically reformulate +# the problem into the correct formulation. + +# From a practical perspective though, the geometric to conic reformulation is +# particularly problematic, because the additional slack variables and +# constraints can make the problem much larger and therefore harder to solve. + +# You should also note many problems contain a mix of conic constraints and +# variables, and so they do not neatly fall into one of the two formulations. In +# these cases, JuMP reformulates only the variables and constraints as necessary +# to convert the problem into the desired form. + +# ## Primal and dual formulations + +# Duality plays a large role in conic optimization. For a detailed description +# of conic duality, see [Duality](@ref). + +# A useful observation is that if the primal problem is in standard conic form, +# then the dual problem is in geometric conic form, and vice versa. Morever, the +# primal and dual may have a different number of variables and constraints, +# although which one is smaller depends on the problem. Therefore, instead of +# reformulating the problem problem from one form to the other, it can be more +# efficient to solve the dual instead of the primal. + +# To demonstrate, we use a variation of the [Maximum cut via SDP](@ref) example. + +# The primal formulation (in standard conic form) is: + +model_primal = Model() +@variable(model_primal, X[1:2, 1:2], PSD) +@objective(model_primal, Max, sum([1 -1; -1 1] .* X)) +@constraint(model_primal, primal_c[i = 1:2], 1 - X[i, i] == 0) +print(model_primal) + +# This problem has three scalar decision variables (the matrix `X` is symmetric), +# two scalar equality constraints, and a constraint that `X` is positive +# semidefinite. + +# The dual of `model_primal` is: + +model_dual = Model() +@variable(model_dual, y[1:2]) +@objective(model_dual, Min, sum(y)) +@constraint(model_dual, dual_c, [y[1]-1 1; 1 y[2]-1] in PSDCone()) +print(model_dual) + +# This problem has two scalar decision variables, and a 2x2 positive +# semidefinite matrix constraint. + +# !!! tip +# If you haven't seen conic duality before, try deriving the dual problem +# based on the description in [Duality](@ref). You'll need to know that the +# dual cone of [`PSDCone`](@ref) is the [`PSDCone`](@ref). + +# When we solve `model_primal` with `SCS.Optimizer`, SCS reports three variables +# (`variables n: 3`), five rows in the constraint matrix (`constraints m: 5`), +# and five non-zeros in the matrix (`nnz(A): 5`): + +set_optimizer(model_primal, SCS.Optimizer) +optimize!(model_primal) + +# (There are five rows in the constraint matrix because SCS expects problems in +# geometric conic form, and so JuMP has reformulated the `X, PSD` variable +# constraint into the affine constraint `X .+ 0 in PSDCone()`.) + +# The solution we obtain is: + +value.(X) + +#- + +dual.(primal_c) + +#- + +objective_value(model_primal) + +# When we solve `model_dual` with `SCS.Optimizer`, SCS reports two variables +# (`variables n: 2`), three rows in the constraint matrix (`constraints m: 3`), +# and two non-zeros in the matrix (`nnz(A): 2`): + +set_optimizer(model_dual, SCS.Optimizer) +optimize!(model_dual) + +# and the solution we obtain is: + +dual.(dual_c) + +#- + +value.(y) + +#- + +objective_value(model_dual) + +# The problems are small enough that it isn't meaningful to compare the solve +# times, but in general, we should expect the dual problem to solve faster +# because it contains fewer variables and constraints. The difference is +# particularly noticeable on large-scale optimization problems. + +# ## `dual_optimizer` + +# Manually deriving the conic dual is difficult and error-prone. The package +# [Dualization.jl](@ref) provides the `Dualization.dual_optimizer` meta-solver, +# which wraps any MathOptInterface-compatible solver in an interface that +# automatically formulates the dual of an input problem, solves the dual +# problem, and then reports the primal solution back to the user. + +# To demonstrate, we use `Dualization.dual_optimizer` to solve `model_primal`; + +set_optimizer(model_primal, Dualization.dual_optimizer(SCS.Optimizer)) +optimize!(model_primal) + +# The performance is the same as if we solved `model_dual`, and the correct +# solution is returned to `X`: + +value.(X) + +#- + +dual.(primal_c) + +# Moreover, if we use `dual_optimizer` on `model_dual`, then we get the same +# performance as if we had solved `model_primal`: + +set_optimizer(model_dual, Dualization.dual_optimizer(SCS.Optimizer)) +optimize!(model_dual) +dual.(dual_c) +value.(y) + +# ## When to use `dual_optimizer` + +# Because it can make the problem larger or smaller, depending on the problem +# and the choice of solver, there is is no definitive rule on when you should +# use `dual_optimizer`. However, you should try `dual_optimizer` if your conic +# optimization problem takes a long time to solve, or if you need to repeatedly +# solve similarly structured problems with different data. In some cases solving +# the dual instead of the primal can make a large difference. + From 79e511851eb2180f902277ec167e535bce570e1e Mon Sep 17 00:00:00 2001 From: odow Date: Mon, 5 Jun 2023 16:48:14 +1200 Subject: [PATCH 02/13] Formatting --- docs/src/tutorials/conic/dualization.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/docs/src/tutorials/conic/dualization.jl b/docs/src/tutorials/conic/dualization.jl index 4323e16e420..9747e87ee4b 100644 --- a/docs/src/tutorials/conic/dualization.jl +++ b/docs/src/tutorials/conic/dualization.jl @@ -19,7 +19,6 @@ using JuMP import Dualization -import LinearAlgebra import SCS # ## Background @@ -207,4 +206,3 @@ value.(y) # optimization problem takes a long time to solve, or if you need to repeatedly # solve similarly structured problems with different data. In some cases solving # the dual instead of the primal can make a large difference. - From 48ec447f7239a060393657570b085f83faaee291 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Mon, 5 Jun 2023 17:20:03 +1200 Subject: [PATCH 03/13] Apply suggestions from code review --- docs/src/tutorials/conic/dualization.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/tutorials/conic/dualization.jl b/docs/src/tutorials/conic/dualization.jl index 9747e87ee4b..aabc7e24618 100644 --- a/docs/src/tutorials/conic/dualization.jl +++ b/docs/src/tutorials/conic/dualization.jl @@ -5,7 +5,7 @@ # # Dualization -# The purpose of this tutorial is to explain how to use [Dualization.j](@ref) to +# The purpose of this tutorial is to explain how to use [Dualization.jl](@ref) to # improve the performance of some conic optimization models. There are two # important takeaways: # @@ -23,7 +23,7 @@ import SCS # ## Background -# Conic optimization solvers typically accept one of two input formulations: +# Conic optimization solvers typically accept one of two input formulations. # The first is the _standard_ conic form: # ```math @@ -88,7 +88,7 @@ import SCS # then the dual problem is in geometric conic form, and vice versa. Morever, the # primal and dual may have a different number of variables and constraints, # although which one is smaller depends on the problem. Therefore, instead of -# reformulating the problem problem from one form to the other, it can be more +# reformulating the problem from one form to the other, it can be more # efficient to solve the dual instead of the primal. # To demonstrate, we use a variation of the [Maximum cut via SDP](@ref) example. @@ -176,7 +176,7 @@ objective_value(model_dual) # automatically formulates the dual of an input problem, solves the dual # problem, and then reports the primal solution back to the user. -# To demonstrate, we use `Dualization.dual_optimizer` to solve `model_primal`; +# To demonstrate, we use `Dualization.dual_optimizer` to solve `model_primal`: set_optimizer(model_primal, Dualization.dual_optimizer(SCS.Optimizer)) optimize!(model_primal) From a53d5a9b91e05c9506031b8ade4b5894d16255ea Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Mon, 5 Jun 2023 20:39:51 +1200 Subject: [PATCH 04/13] Apply suggestions from code review Co-authored-by: James Foster <38274066+jd-foster@users.noreply.github.com> --- docs/src/tutorials/conic/dualization.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/tutorials/conic/dualization.jl b/docs/src/tutorials/conic/dualization.jl index aabc7e24618..a6b76df1fae 100644 --- a/docs/src/tutorials/conic/dualization.jl +++ b/docs/src/tutorials/conic/dualization.jl @@ -53,7 +53,7 @@ import SCS # \begin{align} # \min_{x \in \mathbb{R}^n} & c^\top x \\ # \;\;\text{s.t.} & [A -I] [x; y] = b \\ -# & [x; y] \in \mathbb{R}^n \cross \mathcal{K} +# & [x; y] \in \mathbb{R}^n \times \mathcal{K} # \end{align} # ``` # and to go from the standard conic form to the geometric conic form, we can @@ -61,7 +61,7 @@ import SCS # ```math # \begin{align} # \min_{x \in \mathbb{R}^n} & c^\top x \\ -# \;\;\text{s.t.} & [A; I] x - [b; 0] \in {0}\cross \mathcal{K} +# \;\;\text{s.t.} & [A; I] x - [b; 0] \in \{0\} \times \mathcal{K} # \end{align} # ``` From 58a1a56b3206c1280b162e92b392b90234387469 Mon Sep 17 00:00:00 2001 From: James Foster <38274066+jd-foster@users.noreply.github.com> Date: Mon, 5 Jun 2023 18:57:45 +1000 Subject: [PATCH 05/13] Alignment --- docs/src/tutorials/conic/dualization.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/src/tutorials/conic/dualization.jl b/docs/src/tutorials/conic/dualization.jl index a6b76df1fae..178ba9c85b4 100644 --- a/docs/src/tutorials/conic/dualization.jl +++ b/docs/src/tutorials/conic/dualization.jl @@ -28,8 +28,8 @@ import SCS # The first is the _standard_ conic form: # ```math # \begin{align} -# \min_{x \in \mathbb{R}^n} & c^\top x \\ -# \;\;\text{s.t.} & A x = b \\ +# \min_{x \in \mathbb{R}^n} \; & c^\top x \\ +# \;\;\text{s.t.} \; & A x = b \\ # & x \in \mathcal{K} # \end{align} # ``` @@ -39,8 +39,8 @@ import SCS # The second is the _geometric_ conic form: # ```math # \begin{align} -# \min_{x \in \mathbb{R}^n} & c^\top x \\ -# \;\;\text{s.t.} & A x - b \in \mathcal{K} +# \min_{x \in \mathbb{R}^n} \; & c^\top x \\ +# \;\;\text{s.t.} \; & A x - b \in \mathcal{K} # \end{align} # ``` # in which an affine function $Ax - b$ belongs to a cone $\mathcal{K}$ and the @@ -51,8 +51,8 @@ import SCS # variables $y$: # ```math # \begin{align} -# \min_{x \in \mathbb{R}^n} & c^\top x \\ -# \;\;\text{s.t.} & [A -I] [x; y] = b \\ +# \min_{x \in \mathbb{R}^n} \; & c^\top x \\ +# \;\;\text{s.t.} \; & [A -I] [x; y] = b \\ # & [x; y] \in \mathbb{R}^n \times \mathcal{K} # \end{align} # ``` From 3efeec3114d0b30faa8c7e2c4fd2d990c569591c Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Mon, 5 Jun 2023 21:15:00 +1200 Subject: [PATCH 06/13] Apply suggestions from code review Co-authored-by: James Foster <38274066+jd-foster@users.noreply.github.com> --- docs/src/tutorials/conic/dualization.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/tutorials/conic/dualization.jl b/docs/src/tutorials/conic/dualization.jl index 178ba9c85b4..b990195e0c4 100644 --- a/docs/src/tutorials/conic/dualization.jl +++ b/docs/src/tutorials/conic/dualization.jl @@ -15,7 +15,7 @@ # 2. Solving the dual of a conic model can be more efficient than solving the # primal. -# This tutorial uses the following packges +# This tutorial uses the following packages using JuMP import Dualization @@ -60,7 +60,7 @@ import SCS # rewrite the equality constraint as a function belonging to the `{0}` cone: # ```math # \begin{align} -# \min_{x \in \mathbb{R}^n} & c^\top x \\ +# \min_{x \in \mathbb{R}^n} \; & c^\top x \\ # \;\;\text{s.t.} & [A; I] x - [b; 0] \in \{0\} \times \mathcal{K} # \end{align} # ``` @@ -85,7 +85,7 @@ import SCS # of conic duality, see [Duality](@ref). # A useful observation is that if the primal problem is in standard conic form, -# then the dual problem is in geometric conic form, and vice versa. Morever, the +# then the dual problem is in geometric conic form, and vice versa. Moreover, the # primal and dual may have a different number of variables and constraints, # although which one is smaller depends on the problem. Therefore, instead of # reformulating the problem from one form to the other, it can be more From d6cd52322a7b1e9b749a23903a9ba436b7a70efd Mon Sep 17 00:00:00 2001 From: James Foster <38274066+jd-foster@users.noreply.github.com> Date: Mon, 5 Jun 2023 20:10:25 +1000 Subject: [PATCH 07/13] Update docs/src/tutorials/conic/dualization.jl --- docs/src/tutorials/conic/dualization.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/tutorials/conic/dualization.jl b/docs/src/tutorials/conic/dualization.jl index b990195e0c4..3335025db1b 100644 --- a/docs/src/tutorials/conic/dualization.jl +++ b/docs/src/tutorials/conic/dualization.jl @@ -201,7 +201,7 @@ value.(y) # ## When to use `dual_optimizer` # Because it can make the problem larger or smaller, depending on the problem -# and the choice of solver, there is is no definitive rule on when you should +# and the choice of solver, there is no definitive rule on when you should # use `dual_optimizer`. However, you should try `dual_optimizer` if your conic # optimization problem takes a long time to solve, or if you need to repeatedly # solve similarly structured problems with different data. In some cases solving From cb5ef4241f81591239388945c801af082d9328ba Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Tue, 6 Jun 2023 09:37:26 +1200 Subject: [PATCH 08/13] Update docs/src/tutorials/conic/dualization.jl --- docs/src/tutorials/conic/dualization.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/src/tutorials/conic/dualization.jl b/docs/src/tutorials/conic/dualization.jl index 3335025db1b..8e5fa4d41ce 100644 --- a/docs/src/tutorials/conic/dualization.jl +++ b/docs/src/tutorials/conic/dualization.jl @@ -127,6 +127,7 @@ print(model_dual) set_optimizer(model_primal, SCS.Optimizer) optimize!(model_primal) +nothing # (There are five rows in the constraint matrix because SCS expects problems in # geometric conic form, and so JuMP has reformulated the `X, PSD` variable From b60fc3c512d36164cb52fe2282b702fff94d78db Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Tue, 6 Jun 2023 10:50:14 +1200 Subject: [PATCH 09/13] Update dualization.jl --- docs/src/tutorials/conic/dualization.jl | 29 ++++++++++++++++++++----- 1 file changed, 23 insertions(+), 6 deletions(-) diff --git a/docs/src/tutorials/conic/dualization.jl b/docs/src/tutorials/conic/dualization.jl index 8e5fa4d41ce..4524fa754e3 100644 --- a/docs/src/tutorials/conic/dualization.jl +++ b/docs/src/tutorials/conic/dualization.jl @@ -19,7 +19,7 @@ using JuMP import Dualization -import SCS +import Clarabel # ## Background @@ -125,9 +125,11 @@ print(model_dual) # (`variables n: 3`), five rows in the constraint matrix (`constraints m: 5`), # and five non-zeros in the matrix (`nnz(A): 5`): -set_optimizer(model_primal, SCS.Optimizer) +set_optimizer(model_primal, Clarabel.Optimizer) + +#- + optimize!(model_primal) -nothing # (There are five rows in the constraint matrix because SCS expects problems in # geometric conic form, and so JuMP has reformulated the `X, PSD` variable @@ -149,7 +151,10 @@ objective_value(model_primal) # (`variables n: 2`), three rows in the constraint matrix (`constraints m: 3`), # and two non-zeros in the matrix (`nnz(A): 2`): -set_optimizer(model_dual, SCS.Optimizer) +set_optimizer(model_dual, Clarabel.Optimizer) + +#- + optimize!(model_dual) # and the solution we obtain is: @@ -179,7 +184,10 @@ objective_value(model_dual) # To demonstrate, we use `Dualization.dual_optimizer` to solve `model_primal`: -set_optimizer(model_primal, Dualization.dual_optimizer(SCS.Optimizer)) +set_optimizer(model_primal, Dualization.dual_optimizer(Clarabel.Optimizer)) + +#- + optimize!(model_primal) # The performance is the same as if we solved `model_dual`, and the correct @@ -194,9 +202,18 @@ dual.(primal_c) # Moreover, if we use `dual_optimizer` on `model_dual`, then we get the same # performance as if we had solved `model_primal`: -set_optimizer(model_dual, Dualization.dual_optimizer(SCS.Optimizer)) +set_optimizer(model_dual, Dualization.dual_optimizer(Clarabel.Optimizer)) + +#- + optimize!(model_dual) + +#- + dual.(dual_c) + +#- + value.(y) # ## When to use `dual_optimizer` From 9b4c76f6a038fe8043b9e7adb6c15c0292ea4aa6 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Tue, 6 Jun 2023 11:16:04 +1200 Subject: [PATCH 10/13] Update dualization.jl --- docs/src/tutorials/conic/dualization.jl | 26 +++++++++---------------- 1 file changed, 9 insertions(+), 17 deletions(-) diff --git a/docs/src/tutorials/conic/dualization.jl b/docs/src/tutorials/conic/dualization.jl index 4524fa754e3..937be32ff2b 100644 --- a/docs/src/tutorials/conic/dualization.jl +++ b/docs/src/tutorials/conic/dualization.jl @@ -19,7 +19,7 @@ using JuMP import Dualization -import Clarabel +import SCS # ## Background @@ -125,11 +125,9 @@ print(model_dual) # (`variables n: 3`), five rows in the constraint matrix (`constraints m: 5`), # and five non-zeros in the matrix (`nnz(A): 5`): -set_optimizer(model_primal, Clarabel.Optimizer) - -#- - +set_optimizer(model_primal, SCS.Optimizer) optimize!(model_primal) +Base.Libc.flush_cstdio() #hide # (There are five rows in the constraint matrix because SCS expects problems in # geometric conic form, and so JuMP has reformulated the `X, PSD` variable @@ -151,11 +149,9 @@ objective_value(model_primal) # (`variables n: 2`), three rows in the constraint matrix (`constraints m: 3`), # and two non-zeros in the matrix (`nnz(A): 2`): -set_optimizer(model_dual, Clarabel.Optimizer) - -#- - +set_optimizer(model_dual, SCS.Optimizer) optimize!(model_dual) +Base.Libc.flush_cstdio() #hide # and the solution we obtain is: @@ -184,11 +180,9 @@ objective_value(model_dual) # To demonstrate, we use `Dualization.dual_optimizer` to solve `model_primal`: -set_optimizer(model_primal, Dualization.dual_optimizer(Clarabel.Optimizer)) - -#- - +set_optimizer(model_primal, Dualization.dual_optimizer(SCS.Optimizer)) optimize!(model_primal) +Base.Libc.flush_cstdio() #hide # The performance is the same as if we solved `model_dual`, and the correct # solution is returned to `X`: @@ -202,11 +196,9 @@ dual.(primal_c) # Moreover, if we use `dual_optimizer` on `model_dual`, then we get the same # performance as if we had solved `model_primal`: -set_optimizer(model_dual, Dualization.dual_optimizer(Clarabel.Optimizer)) - -#- - +set_optimizer(model_dual, Dualization.dual_optimizer(SCS.Optimizer)) optimize!(model_dual) +Base.Libc.flush_cstdio() #hide #- From 22719fac931db2592f228f2481660876e8d7d7fc Mon Sep 17 00:00:00 2001 From: odow Date: Tue, 6 Jun 2023 13:55:38 +1200 Subject: [PATCH 11/13] Add tips --- docs/src/tutorials/conic/introduction.md | 6 ++++++ docs/src/tutorials/conic/quantum_discrimination.jl | 6 ++++++ 2 files changed, 12 insertions(+) diff --git a/docs/src/tutorials/conic/introduction.md b/docs/src/tutorials/conic/introduction.md index a0689769bb1..b7ca2235374 100644 --- a/docs/src/tutorials/conic/introduction.md +++ b/docs/src/tutorials/conic/introduction.md @@ -24,6 +24,12 @@ MathOptInterface, many of these solvers support a much wider range of exotic cones than they natively support. Solvers supporting discrete variables start with "(MI)" in the list of [Supported solvers](@ref). +!!! tip + Duality plays a large role in solving conic optimization models. Depending + on the solver, it can be more efficient to solve the dual instead of the + primal. If performance is an issue, see the [Dualization](@ref) tutorial for + more details. + ## How these tutorials are structured Having a high-level overview of how this part of the documentation is structured diff --git a/docs/src/tutorials/conic/quantum_discrimination.jl b/docs/src/tutorials/conic/quantum_discrimination.jl index 25550d5f3d8..d743cc0b005 100644 --- a/docs/src/tutorials/conic/quantum_discrimination.jl +++ b/docs/src/tutorials/conic/quantum_discrimination.jl @@ -114,6 +114,12 @@ objective_value(model) solution = [value.(e) for e in E] +# !!! tip +# Duality plays a large role in solving conic optimization models. Depending +# on the solver, it can be more efficient to solve the dual of this problem +# instead of the primal. If performance is an issue, see the [Dualization](@ref) +# tutorial for more details. + # ## Alternative formulation # The formulation above includes `N` Hermitian matrices and a set of linear From 4290fbebe573c9ecd64be6ce9ab01edd29c4f4a7 Mon Sep 17 00:00:00 2001 From: James Foster Date: Tue, 6 Jun 2023 17:12:28 +1000 Subject: [PATCH 12/13] Add dualization comparison --- docs/src/tutorials/conic/simple_examples.jl | 72 +++++++++++++++++---- 1 file changed, 58 insertions(+), 14 deletions(-) diff --git a/docs/src/tutorials/conic/simple_examples.jl b/docs/src/tutorials/conic/simple_examples.jl index d3015150d61..1a414f3e094 100644 --- a/docs/src/tutorials/conic/simple_examples.jl +++ b/docs/src/tutorials/conic/simple_examples.jl @@ -11,12 +11,17 @@ # This tutorial makes use of the following packages: using JuMP +import Dualization import LinearAlgebra import Plots import Random import SCS import Test +# We will define a dictionary to collect our simple models for later inspection: + +sdp_model_dict = Dict{Symbol,JuMP.Model}() + # ## Maximum cut via SDP # The [maximum cut problem](https://en.wikipedia.org/wiki/Maximum_cut) is a @@ -78,7 +83,7 @@ function solve_max_cut_sdp(weights) T = findall(.!cut) println("Solution:") println(" (S, T) = ({", join(S, ", "), "}, {", join(T, ", "), "})") - return S, T + return S, T, model end # Given the graph @@ -87,7 +92,8 @@ end # ``` # The solution is `(S, T) = ({1}, {2})` -S, T = solve_max_cut_sdp([0 5; 5 0]) +S, T, _ = solve_max_cut_sdp([0 5; 5 0]); +S, T # Given the graph # ```raw @@ -101,7 +107,8 @@ S, T = solve_max_cut_sdp([0 5; 5 0]) # ``` # The solution is `(S, T) = ({1}, {2, 3, 4})` -S, T = solve_max_cut_sdp([0 5 7 6; 5 0 0 1; 7 0 0 1; 6 1 1 0]) +S, T, _ = solve_max_cut_sdp([0 5 7 6; 5 0 0 1; 7 0 0 1; 6 1 1 0]); +S, T # Given the graph # ```raw @@ -115,7 +122,9 @@ S, T = solve_max_cut_sdp([0 5 7 6; 5 0 0 1; 7 0 0 1; 6 1 1 0]) # ``` # The solution is `(S, T) = ({1, 4}, {2, 3})` -S, T = solve_max_cut_sdp([0 1 5 0; 1 0 0 9; 5 0 0 2; 0 9 2 0]) +S, T, sdp_model_dict[:max_cut] = + solve_max_cut_sdp([0 1 5 0; 1 0 0 9; 5 0 0 2; 0 9 2 0]); +S, T # ## K-means clustering via SDP @@ -158,10 +167,10 @@ function example_k_means_clustering() end end end - return + return model end -example_k_means_clustering() +sdp_model_dict[:k_means_clustering] = example_k_means_clustering() # ## The correlation problem @@ -200,10 +209,10 @@ function example_correlation_problem() optimize!(model) println("A lower bound for ρ_AC is $(value(ρ["A", "C"]))") Test.@test value(ρ["A", "C"]) ≈ -0.978 atol = 1e-3 #src - return + return model end -example_correlation_problem() +sdp_model_dict[:correlation_problem] = example_correlation_problem(); # ## The minimum distortion problem @@ -286,7 +295,8 @@ function example_minimum_distortion() Test.@test objective_value(model) ≈ 4 / 3 atol = 1e-4 ## Recover the minimal distorted embedding: X = [zeros(3) sqrt(value.(Q)[2:end, 2:end])] - return Plots.plot( + return model, + Plots.plot( X[1, :], X[2, :], X[3, :]; @@ -304,7 +314,8 @@ function example_minimum_distortion() ) end -example_minimum_distortion() +sdp_model_dict[:minimum_distortion], plot = example_minimum_distortion(); +plot # ## Lovász numbers @@ -378,10 +389,10 @@ function example_theta_problem() Test.@test primal_status(model) == FEASIBLE_POINT Test.@test objective_value(model) ≈ sqrt(5) rtol = 1e-4 println("The Lovász number is: $(objective_value(model))") - return + return model end -example_theta_problem() +sdp_model_dict[:theta_problem] = example_theta_problem(); # ## Robust uncertainty sets @@ -416,7 +427,40 @@ function example_robust_uncertainty_sets() sqrt((1 - ɛ) / ɛ) * sqrt(c' * (Σhat + Γ2(𝛿 / 2, N) * LinearAlgebra.I) * c) Test.@test objective_value(model) ≈ exact atol = 1e-2 - return + return model +end + +sdp_model_dict[:robust_uncertainty_sets] = example_robust_uncertainty_sets(); + +# ## Comparison to dual formulation + +# It may be more efficient to solve the dual of a conic program rather than +# the original formulation. +# Comparing the numbers of variables and constraints in the primal formulations +# given here and in the dual formulation provided by the [Dualization.jl](@ref) package: + +println("Model Name \t #PrimalVars #Cons #DualVars #DualCons") +println("------------------------------------------------------------") +for (sym, model) in sdp_model_dict + # Skip models with constraints types not yet supported: + if sym in [:correlation_problem, :robust_uncertainty_sets] + continue + end + println( + rpad(string(sym), 20), + "\t", + num_variables(model), + "\t", + num_constraints(model; count_variable_in_set_constraints = true), + "\t", + num_variables(Dualization.dualize(model)), + "\t", + num_constraints( + Dualization.dualize(model); + count_variable_in_set_constraints = true, + ), + ) end -example_robust_uncertainty_sets() +# we observe that the dual formulation can have fewer variables and constraints. +# For more details see the [Dualization](@ref) tutorial . From fda087beeb0ccb2f6c7f37dd79d380f05fa1bf79 Mon Sep 17 00:00:00 2001 From: James Foster Date: Tue, 6 Jun 2023 21:33:37 +1000 Subject: [PATCH 13/13] Fixes, formatting, comment response --- docs/src/tutorials/conic/simple_examples.jl | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/docs/src/tutorials/conic/simple_examples.jl b/docs/src/tutorials/conic/simple_examples.jl index 1a414f3e094..e921a76e580 100644 --- a/docs/src/tutorials/conic/simple_examples.jl +++ b/docs/src/tutorials/conic/simple_examples.jl @@ -170,7 +170,7 @@ function example_k_means_clustering() return model end -sdp_model_dict[:k_means_clustering] = example_k_means_clustering() +sdp_model_dict[:k_means_clustering] = example_k_means_clustering(); # ## The correlation problem @@ -436,16 +436,17 @@ sdp_model_dict[:robust_uncertainty_sets] = example_robust_uncertainty_sets(); # It may be more efficient to solve the dual of a conic program rather than # the original formulation. -# Comparing the numbers of variables and constraints in the primal formulations +# We can compare the numbers of variables and constraints in the primal formulations # given here and in the dual formulation provided by the [Dualization.jl](@ref) package: println("Model Name \t #PrimalVars #Cons #DualVars #DualCons") println("------------------------------------------------------------") for (sym, model) in sdp_model_dict - # Skip models with constraints types not yet supported: + ## Skip models with constraints types not yet supported: if sym in [:correlation_problem, :robust_uncertainty_sets] continue end + dual_model = Dualization.dualize(model) println( rpad(string(sym), 20), "\t", @@ -453,14 +454,11 @@ for (sym, model) in sdp_model_dict "\t", num_constraints(model; count_variable_in_set_constraints = true), "\t", - num_variables(Dualization.dualize(model)), + num_variables(dual_model), "\t", - num_constraints( - Dualization.dualize(model); - count_variable_in_set_constraints = true, - ), + num_constraints(dual_model; count_variable_in_set_constraints = true), ) end -# we observe that the dual formulation can have fewer variables and constraints. +# We observe that the dual formulation can have fewer variables and constraints. # For more details see the [Dualization](@ref) tutorial .