How to construct and plot a polyhedron with absolute values?

[WuSiren]

First I would like to say thanks to this useful package and its developers!

Now I'm faced with a question on the construction and visualization of a polyhedron with absolute value structure, e.g.,

$x\in\mathbb{R}^2: |w_1^\top x-b_1| + |w_2^\top x-b_2| \leq c$ or generally, $x\in\mathbb{R}^n: \sum_i k_i |w_i^\top x-b_i| \leq c$

In the following example, I've tried three method to construct and visualize the polyhedron but none of them satisfied me. As for Method-1, it can only draw the contour of the polyhedron and it seems not to be quite exact. Method-2 is more straightforward with the help of the packages JuMP.jl and Polyhedra.jl, but it will report an error when converting the JuMP model into a Polyhedron:

ERROR: MathOptInterface.UnsupportedConstraint{MathOptInterface.ScalarNonlinearFunction, MathOptInterface.LessThan{Float64}}: `MathOptInterface.ScalarNonlinearFunction`-in-`MathOptInterface.LessThan{Float64}` constraint is not supported by the model.

Then I tried Method-3 by introducing two new auxiliary variables. It does work, but I have to eliminate the irrelevant auxiliary variables before visualizing the polyhedron, which will be much more time consuming especially when there are lots of absolute value structures. (And by the way, can it display the contour of the polyhedron with another color?)

So what will be the most recommended method to deal with this kind of polyhedra? Or is there any other better method?

Looking forward to your reply at your earliest convenience!

using Polyhedra
using CairoMakie
using JuMP

w₁ = [1, 2]
w₂ = [3, 1]
b₁ = -1
b₂ = 2
f(x) = abs(w₁' * x - b₁) + abs(w₂' * x - b₂)
c = 3

fig = Figure()

# Method-1:
ax = Axis(fig[1, 1]; aspect=4/6)
xlims!(ax, [-1, 3])
ylims!(ax, [-4, 2])
xs = LinRange(-1, 3, 100)
ys = LinRange(-4, 2, 100)
zs = [f([x, y]) for x in xs, y in ys]
contour!(xs, ys, zs, levels = [c], color=:red)
fig

# Method-2:
model = Model()
@variable(model, x[1:2])
@constraint(model, f(x) <= c)
p = polyhedron(model)
plot(p, ratio=:equal)
m = Polyhedra.Mesh(p)
mesh(m, color=:blue)

# Method-3:
model = Model()
@variable(model, x[1:2])
@variable(model, v₁)
@variable(model, v₂)
@constraint(model, -v₁ <= w₁' * x - b₁)
@constraint(model, w₁' * x - b₁ <= v₁)
@constraint(model, -v₂ <= w₂' * x - b₂)
@constraint(model, w₂' * x - b₂ <= v₂)
@constraint(model, v₁ + v₂ <= c)
all_variables(model)
P = polyhedron(model)
p = eliminate(P, [3,4])

ax = Axis(fig[1, 2]; aspect=4/6)
xlims!(ax, [-1, 3])
ylims!(ax, [-4, 2])
m = Polyhedra.Mesh(p)
mesh!(ax, m, color=:blue)
fig

[blegat]

As for the support of nonlinear functions, you could use Convex.jl to rewrite the JuMP nonlinear model once v0.16 is released, see https://github.com/jump-dev/Convex.jl/?tab=readme-ov-file#using-with-jump, I should write an example how to do it. But this will also create auxiliary variables. You could also use the MOI.NormOneCone:

@constraint(model, [c, w₁' * x - b₁, w₂' * x - b₂] in MOI.NormOneCone(3))

but that will also create auxiliary variables.

You can directly create the projection and not create any auxiliary variables by doing:

@constraint(model, w₁' * x - b₁ + w₂' * x - b₂ <= c)
@constraint(model, -(w₁' * x - b₁) + w₂' * x - b₂ <= c)
@constraint(model, w₁' * x - b₁ - (w₂' * x - b₂) <= c)
@constraint(model, -(w₁' * x - b₁) - (w₂' * x - b₂) <= c)

This creates 2^n half-spaces so that doesn't scale so well but that's how the cross-polytope is :) The hypercube has 2n half-spaces and 2^n points while its polar, the cross-polytope has 2n points but 2^n half-spaces. This is precisely why in optimization, we use the fact that by adding a few variables, this lifted cross-polytope has a small H-representation instead of using the H-representation of the cross-polytope which has 2^n. As for the color, doesn't the Makie argument for color work ?

[WuSiren]

Thank you so much, @blegat ! I'm very excited to see your reply. It's very professional and comprehensive, and I feel inspired after reading it. Please accept my respect!

[blegat]

Thanks for your enthusiasm! Let me know if there is still something we can do to help or if the issue can be closed.

[WuSiren]

No thanks! This issue can be closed.