Approximating a $L^2$ Ball with Linear Constraints¶

using JuMP
# We need to have a solver that supports lazy
# constraints. Options include GLPK, Gurobi, and
# CPLEX - others may be come available in the future.
#using GLPKMathProgInterface
using Gurobi

# solve_ball
# Optimizes the linear function c'x
# where ||x|| ≤ Γ + ϵ, where ϵ is an
# absolute feasibility tolerance.
function solve_ball(c, Γ, ϵ=1e-6)
    # Size of vector
    n = length(c)
    
    # Create the model
    # Make sure to change this line!
    #m = Model(solver=GLPKSolverMIP())
    m = Model(solver=GurobiSolver(OutputFlag=0))
    
    # Put the variables in a box
    # This ensures the initial solution, before any
    # linear constraints are generated, is bounded.
    @variable(m, -Γ ≤ x[1:n] ≤ +Γ, Int)
    
    # Objective is simply the inner product
    @objective(m, Max, dot(c,x))
    
    # We'll now build a callback function
    # Note that the callback is an inner function,
    # or closure. This means that it inherits the
    # scope of the enclosing function, including
    # the definition of the variables x
    # A callback must take a single argument, the
    # callback handle. When we want to communicate
    # back to the solver, we do so through the
    # callback handle
    # We'll keep track of how many times the
    # callback was called, just for interests sake
    num_callbacks = 0
    function norm_callback(cb)
        num_callbacks += 1
        # First, lets get the value of x at the
        # current solution. Since we the linear
        # constraint we are adding is going to be
        # a tangent hyperplane of the L2 ball at
        # this point, we can think of the value of
        # x as the normal of that hyperplane
        N = getvalue(x)
        
        # Lets get the length of the normal
        L = norm(N)
        
        # If the length is small enough...
        if L ≤ Γ + ϵ
            # We can stop right now
            return
        end
        
        # The constraint is violated, so we need to
        # add a tangent hyperplane to the L2 ball
        # at this point. We have the normal of the
        # plane, and the point we touch the sphere
        # is at (Γ/L)*N, so the plane equation is
        # <N,x> = <N,x_0> = <N,(Γ/L)*N>
        #                 = <N,(Γ/||N||)*N>
        #                 = Γ*||N||
        # Note we don't use @constraint(m, ...)!
        @lazyconstraint(cb, dot(N,x) ≤ Γ*L)
    end
    
    # We now tell JuMP/the solver to use this callback
    addlazycallback(m, norm_callback)
    
    # Now solve it!
    solve(m)
    
    # Return the solution
    return getvalue(x), num_callbacks
end

# Lets generate a random instance
srand(1234)
# GLPK is very slow with anything bigger
# than 2 or 3, Gurobi should be fine up
# to at least 10
n = 2
c = rand(n)
Γ = 50.0

# Solve and display
sol, num_callbacks = solve_ball(c, Γ)
println(sol)
println(norm(sol))
println(num_callbacks)

# If the problem is small enough, and 2D, we can even
# visualize the solution using Compose.jl and set
# the vector c with Interact.jl's sliders
using Compose, Interact
set_default_graphic_size(8cm, 8cm)

@manipulate for c1 in -1:0.1:+1, c2 in -1:0.1:+1, logϵ in -4:2
    sol, _ = solve_ball([c1,c2], 100, 10.0^logϵ)

    compose(context(),
    # Draw the solution vector
    compose(context(),
            line([(0.5,0.5),(0.5+sol[1]/300,0.5+sol[2]/300)]),
            stroke("black")),
    # Draw the intersection point
    compose(context(),
            circle((0.5 + (100/norm(sol))*sol/300)...,0.02),
            fill("red")),
    # Draw the feasible region
    compose(context(),circle(0.5,0.5,0.333),fill("lightblue"))
    )
end