jump-dev / Convex.jl

A Julia package for disciplined convex programming
https://jump.dev/Convex.jl/stable/
Other
568 stars 123 forks source link

Modify problem in its canonical convex form #702

Closed tcovert closed 3 months ago

tcovert commented 3 months ago

Is there a way to define a problem using the Convex.jl operators, and then somehow modify it via its "canonical" convex form? By this I mean, suppose I've got a problem like this:

  beta = Variable(size(X, 2))
  gamma = Variable(size(Z, 2))
  prb = maximize(
    y' * ((X * beta) .* (Z * gamma)) -
    sum(logsumexp((X[g, :] * beta) .* (Z[g, :] * gamma)) for g in gs),
    gamma >= 0
  )

where X is a matrix of float64s, Z is a matrix of float64s that are all 0 or 1, X and Z have the same number of rows, and gs is vector of int64 vectors (a set of indices into the rows of X/Z). This program is not convex, but holding gamma fixed, its convex in beta, and holding beta fixed its convex in gamma. I would like to use fix! and free! to do block-coordinate-descent on this problem. The documentation for this works exactly as advertised, so thank you for that.

However, I am also interested in trying out the "slack variable" version of this alternating coordinate descent heuristic algorithm, as suggested in this paper. It seems that I need to do this by way of the canonical form: write out the problem with a linear objective + a set of convex cone constraints, and then add slack variables to the cone constraints. In my head, I imagine there is some function canonical(prb) which delivers the problem in its canonical form, and then perhaps add_convex_constraint! and add_convex_variable!? I don't see these in documentation so they probably don't exist, but maybe I've read the docs wrong. However, it feels like something in this spirit must exist, because Clarabel/SCS and other conic solvers expect problems in this form, and I would assume Convex is creating that form under the hood.

I'm guessing there is a way to work backward on pen and paper and figure out what the equivalent "original" problem would look like, including these slack variables, but I'm hoping there is an easier way.

Thanks in advance for any suggestions you've got.

ericphanson commented 3 months ago

I don't think I have any specific suggestions, but I can at least help with this bit:

In my head, I imagine there is some function canonical(prb) which delivers the problem in its canonical form, and then perhaps add_convex_constraint! and add_convex_variable!? I don't see these in documentation so they probably don't exist, but maybe I've read the docs wrong. However, it feels like something in this spirit must exist, because Clarabel/SCS and other conic solvers expect problems in this form, and I would assume Convex is creating that form under the hood.

Convex.jl uses MathOptInterface (MOI) as a translation layer between the user input and the solver. We apply extended formulations to the problem and create a MOI model. Then MOI reformulates it into whatever the solver can handle (the solvers likewise interface with MOI to express what kind of formulations they can accept).

MOI is very expressive and supports a lot of possibilities, so there isn’t really 1 canonical model, but rather we express the model it it’s language and it reformulates it as efficiently as it can to the specific solvers needs at solve time.

All that to say: there isn’t quite such a canonical function, and we express membership-in-set constraints directly to MOI, which can choose to formulate them. You can see this problem after solving by doing print(prb.model). For example:


julia> print(prb.model)
Maximize ScalarAffineFunction{Float64}:
 0.0 + 6.237063806526346 v[1] + 4.886300601090418 v[2] + 6.2978507788998135 v[3] + 4.5992429510859445 v[4] + 4.469419614740134 v[5] + 6.594732380173294 v[6] + 5.239946951111149 v[7] + 5.574661751076638 v[8] + 5.426078130556097 v[9] + 4.00900537675019 v[10] - 1.0 v[11]

Subject to:

VectorAffineFunction{Float64}-in-Nonnegatives
 ┌                                                               ┐
 │1.0 - 1.0 v[12] - 1.0 v[13] - 1.0 v[14] - 1.0 v[15] - 1.0 v[16]│
 └                                                               ┘ ∈ Nonnegatives(1)

VectorAffineFunction{Float64}-in-ExponentialCone
 ┌                                                                                                                                                                                                                                                                                     ┐
 │0.0 + 0.6020956881780196 v[1] + 0.989640132821234 v[2] + 1.4115648659655993 v[3] + 0.9397309305735646 v[4] + 1.9986722705699311 v[5] + 1.2178945067117224 v[6] + 0.7174192769987721 v[7] + 0.6522603257230224 v[8] + 1.3825336872267182 v[9] + 0.050200856468531334 v[10] - 1.0 v[11]│
 │1.0                                                                                                                                                                                                                                                                                  │
 │0.0 + 1.0 v[12]                                                                                                                                                                                                                                                                      │
 └                                                                                                                                                                                                                                                                                     ┘ ∈ ExponentialCone()
 ┌                                                                                                                                                                                                                                                                                     ┐
 │0.0 + 1.7587700523004042 v[1] + 1.3796867946400135 v[2] + 1.3286545789995903 v[3] + 0.4534184241823647 v[4] + 1.1941900433888082 v[5] + 1.8353598267734466 v[6] + 0.6783330861439868 v[7] + 1.6628698625204938 v[8] + 1.2714881622050405 v[9] + 0.06355482986590741 v[10] - 1.0 v[11]│
 │1.0                                                                                                                                                                                                                                                                                  │
 │0.0 + 1.0 v[13]                                                                                                                                                                                                                                                                      │
 └                                                                                                                                                                                                                                                                                     ┘ ∈ ExponentialCone()
 ┌                                                                                                                                                                                                                                                                                     ┐
 │0.0 + 1.3326247884339864 v[1] + 1.6770984642757902 v[2] + 0.7749914863607495 v[3] + 0.509492574143402 v[4] + 0.42723564557585625 v[5] + 0.28650556342022176 v[6] + 0.8351489443717468 v[7] + 0.5636154250788334 v[8] + 0.2115226902614478 v[9] + 0.2948131986437591 v[10] - 1.0 v[11]│
 │1.0                                                                                                                                                                                                                                                                                  │
 │0.0 + 1.0 v[14]                                                                                                                                                                                                                                                                      │
 └                                                                                                                                                                                                                                                                                     ┘ ∈ ExponentialCone()
 ┌                                                                                                                                                                                                                                                                                    ┐
 │0.0 + 0.9284195939049624 v[1] + 0.4424947206777909 v[2] + 0.8767693532773614 v[3] + 2.36215516348266 v[4] + 0.040640121534430064 v[5] + 1.4430410099107005 v[6] + 1.3051005396047795 v[7] + 1.0880221224906905 v[8] + 0.7912394433405268 v[9] + 0.9390738684104556 v[10] - 1.0 v[11]│
 │1.0                                                                                                                                                                                                                                                                                 │
 │0.0 + 1.0 v[15]                                                                                                                                                                                                                                                                     │
 └                                                                                                                                                                                                                                                                                    ┘ ∈ ExponentialCone()
 ┌                                                                                                                                                                                                                                                                                       ┐
 │0.0 + 0.33068717611703524 v[1] + 0.1291860725841855 v[2] + 0.7460500824528628 v[3] + 0.33813835499767597 v[4] + 0.4666657777219453 v[5] + 0.17550801391769325 v[6] + 0.7877036737913496 v[7] + 0.606140934670409 v[8] + 0.15747849986805068 v[9] + 0.5642840440743772 v[10] - 1.0 v[11]│
 │1.0                                                                                                                                                                                                                                                                                    │
 │0.0 + 1.0 v[16]                                                                                                                                                                                                                                                                        │
 └                                                                                                                                                                                                                                                                                       ┘ ∈ ExponentialCone()

(BTW: next time include some random data so it's easier to copy-paste to run your code).

If you don't want to solve first, you can do

import MathOptInterface as MOI
const MOIU = MOI.Utilities
context = Convex.Context(prb, MOI.Utilities.Model{Float64});
print(context.model)

This does Convex's reformulations but does not solve the model.

odow commented 3 months ago

write out the problem with a linear objective + a set of convex cone constraints, and then add slack variables to the cone constraints

This is one possible standard form, but it is not the only one :smile:

The upside of using MathOptInterface is that it handles the reformulation to standard form and is very flexible, letting us match what is required by the solver.

The downside is that there isn't one-true-standard-form Ax - s = b, s in K, which can make "simple" things like "give me the matrices" more complicated.

The good thing though is that if you write your algorithm as a MOI optimizer with declared support for a few cones, then all models in Convex.jl that can be reformulated to your solver will work.

Since the documentation/tutorials can be a bit sparse and there are a number of possible ways to implement a MOI optimizer, perhaps the best suggestion is to follow: https://github.com/blegat/Zaphod.jl.

You may find these bits helpful:

For a similar problem with LPs, see:

odow commented 3 months ago

Since this is not a bug in Convex.jl, and it is not a feature that we are going to add, I think I am going to close this issue as "won't-fix".

@tcovert: if you have questions after reading the links, please make a post on our community forum https://discourse.julialang.org/c/domain/opt/13 and I am happy keep discussing over there :smile: