Can this NormOneCone be made easier? as SecondOrderCone is?
It will find it use in linear programming.
In my code, I now write something like
nm1e = JuMP.@variable(outer)
nm1c = JuMP.@constraint(outer, [nm1e; vec(outer_β - β)] ∈ JuMP.MOI.NormOneCone(1+length(β)))
I wish I could write
nm1c = JuMP.@constraint(outer, [nm1e; vec(outer_β - β)] ∈ JuMP.NormOneCone())
in the future. Will this be okay?
There are arguments for and against this.
For:
Against:
PSDCone() etc, so that we can convert between the vector in MOI and the symmetric matrix in JuMP. That doesn’t apply here.SecondOrderCone because we did. Mainly because it is quite commonly used. I don’t know that I would add it if I could choose again today.JuMP.NormOneCone, what about all the other MOI cones? I don’t really want to have so much duplication, especially just to save a few characters of typing.We have NormOneCone documented here: Tips and tricks · JuMP. I think I’d prefer to make it really clear that you can use MOI cones directly in JuMP. They’re not second class citizens.
Note that if you do import MathOptInterface as MOI you don’t need JuMP.MOI.
Perhaps one could redefine the MOI cones to not need to give the size of the input as a argument then?
That is a reasonable suggestion but it would be very breaking so its not going to happen.
How would it be breaking? I’m suggesting to make giving the input size optional, not removing it entirely. One can still use it for error checking.
It is breaking because various codes assume you can call MOI.dimension(set).
You can add a new optional data attribute. You can’t take away a previously required attribute.
One would still create the set with the appropriate dimension, and it would still be readable with MOI.dimension. It would just be computed from the given input vector instead of given explicitly as an argument.
Thats exactly the behavior of the JuMP Abstract Set. Which we can do. But then there are two copies of every set. One nondimensional in JuMP and a dimensional set in MOI. For the sake of a few characters, id prefer we had less code in JuMP.
I still don’t see the difficulty with creating a dimensional set without giving it the dimension explicitly. @variable and @constraint are macros, they can handle the necessary rewriting.
Otherwise MOI sets are in fact second-class citizens.
There is no difficulty. It’s what JuMP.AbstractSet already supports: Extensions · JuMP
The code required would be:
struct NormOneCone <: JuMP.AbstractVectorSet end
moi_set(set::NormOneCone, dim::Int) = MOI.NormOneCone(dim)
but then I need to implement, test, document and maintain two copies of every set. Just so people don’t have to specify the dimension in their JuMP model. That’s certainly a worthy trade-off up for discussion. But the costs (to me) don’t seem to outweigh the minimal costs to users.
Otherwise MOI sets are in fact second-class citizens.
My point is that they aren’t.
You just need to write MOI.NormOneCone(3) instead of JuMP.NormOneCone().
We should teach people that they can use any MOI set, not just a subset of ones that have a corresponding JuMP.AbstractSet.
Having two copies of each set would indeed be silly. But I’m certain that there is a sensible way to implement this feature, either by redefining the MOI sets themselves, or redefining the macros.
In any case, there is no point in arguing, since you’re clearly not interested in doing it and I’m not going to do it myself.
I find a new (perhaps better) method to write this NormOneCone constraint, when the model is intended to be revise-and-solved iteratively. @odow is this good? Give me some advice
import LinearAlgebra.norm as norm
import JuMP, Gurobi
I = 2; # `I` is the length of the decision vector `x`
x_old = rand(2); # The new solution `value(x)` after `optimize!` should stay close to `x_old`
# In this code file, I propose to use formulation `reg2` ✅, instead of the existing `reg0`
# `reg0` is the original version of regularization subprogram
reg0 = JuMP.Model(Gurobi.Optimizer); JuMP.set_silent(reg0);
JuMP.@variable(reg0, x[1:I])
JuMP.@variable(reg0, adx[1:I]) # adx[i] attempts to = abs(x[i] - x_old[i]), ∀i
JuMP.@constraint(reg0, c1[i = 1:I], x[i] - x_old[i] ≤ adx[i])
JuMP.@constraint(reg0, c2[i = 1:I], x_old[i] - x[i] ≤ adx[i])
JuMP.@objective(reg0, Min, 2 * sum(adx)) # it's immaterial to mul a positive coefficiant `2`
JuMP.optimize!(reg0); JuMP.assert_is_solved_and_feasible(reg0; allow_local = false)
distance = norm(JuMP.value.(x) - x_old, 1)
println("distance = $distance")
# We introduce 2 new decisions (m, p) based on the above formulation
reg1 = JuMP.Model()
JuMP.@variable(reg1, x[1:I])
JuMP.@variable(reg1, adx[1:I])
JuMP.@variable(reg1, m[1:I])
JuMP.@variable(reg1, p[1:I])
JuMP.@constraint(reg1, [i = 1:I], adx[i] == m[i] + x[i]) # m[i] := adx[i] - x[i], ∀i
JuMP.@constraint(reg1, [i = 1:I], adx[i] == p[i] - x[i]) # p[i] := adx[i] + x[i], ∀i
# update these 2 constrs
JuMP.@constraint(reg1, c1[i = 1:I], -x_old[i] ≤ m[i])
JuMP.@constraint(reg1, c2[i = 1:I], x_old[i] ≤ p[i])
JuMP.@objective(reg1, Min, sum(adx + adx))
# Finally, we eliminate the free decision `adx` and derive the following ✅
reg2 = JuMP.Model(Gurobi.Optimizer); JuMP.set_silent(reg2);
JuMP.@variable(reg2, x[1:I])
JuMP.@variable(reg2, m[1:I])
JuMP.@variable(reg2, p[1:I])
JuMP.@expression(reg2, new_obj_expr, sum(m + p))
JuMP.@constraint(reg2, [i = 1:I], m[i] + x[i] == p[i] - x[i]) # p[i] := adx[i] + x[i], ∀i
JuMP.set_lower_bound.(m, -x_old) # in place of `c1`
JuMP.set_lower_bound.(p, x_old) # in place of `c2`
JuMP.@objective(reg2, Min, new_obj_expr)
# 🖥️ Let's do some tests on `reg2`
for t in 1:100
x_old = rand(2)
let # update `reg2`
JuMP.set_lower_bound.(m, -x_old)
JuMP.set_lower_bound.(p, x_old)
JuMP.@objective(reg2, Min, new_obj_expr)
end
JuMP.optimize!(reg2); JuMP.assert_is_solved_and_feasible(reg2; allow_local = false)
distance = norm(JuMP.value.(x) - x_old, 1)
iszero(distance) || error("distance is not zero!")
end
The two most common ways to formulate an L1-norm cone are:
using JuMP
model = Model()
@variable(model, x[1:3])
@variable(model, t)
# Option: 1
@variable(model, abs_x[1:3] >= 0)
@constraint(model, abs_x .>= x)
@constraint(model, abs_x .>= -x)
@constraint(model, t == sum(abs_x))
# Option: 2
@variable(model, x_p[1:3] >= 0)
@variable(model, x_n[1:3] >= 0)
@constraint(model, x .== x_p .- x_n)
@constraint(model, t == sum(x_p) + sum(x_n))
Your “n” means the negative part, “p” means positive part.
My “m” means minus, “p” means plus.
My formulation is indeed a generalization of your Option 2, where x_old and @objective is revisable. If the objective is changed to something rather than new_obj_expr, e.g. @objective(reg2, Max, 0), then m, p and new_obj_expr enforce no restrictions to x—meaning that they can be safely left in the reg2 model.