WIP: First version of Polyhedral with QPDAS solver

Project.toml

@@ -8,6 +8,7 @@ OSQP = "ab2f91bb-94b4-55e3-9ba0-7f65df51de79"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e"
+QPDAS = "297e508c-b2a7-11e8-1524-5dccae6124f6"
 
 [compat]
 IterativeSolvers = "0.8.0"

REQUIRE

@@ -2,3 +2,4 @@ julia 1.0
 IterativeSolvers 0.8.0
 TSVD 0.3.0
 OSQP 0.3.0
+QPDAS 0.1.0

src/functions/indPolyhedral.jl

@@ -20,6 +20,8 @@ it is assumed to be (plus or minus) infinity.
 function IndPolyhedral(args...; solver=:osqp)
     if solver == :osqp
         IndPolyhedralOSQP(args...)
+    elseif solver == :qpdas
+        IndPolyhedralQPDAS(args...)
     else
         error("unknown solver")
     end
@@ -28,3 +30,4 @@ end
 # including concrete types
 
 include("indPolyhedralOSQP.jl")
+include("indPolyhedralQPDAS.jl")

src/functions/indPolyhedralQPDAS.jl

@@ -0,0 +1,146 @@
+# IndPolyhedral: QPDAS implementation
+
+import QPDAS
+
+"""
+**Indicator of a polyhedral set**
+
+    IndPolyhedralQPDAS(A, C, b, d)
+
+```math
+S = \\{x : \\langle A, x \\rangle = b\\ ∧ \\langle C, x \\rangle \\le b\\}.
+```
+"""
+struct IndPolyhedralQPDAS{R<:Real, MT<:AbstractMatrix{R}, VT<:AbstractVector{R}, QP<:QPDAS.QuadraticProgram} <: IndPolyhedral
+    A::MT
+    b::VT
+    C::MT
+    d::VT
+    z::VT
+    qp::QP
+    first_prox::Ref{Bool}
+    function IndPolyhedralQPDAS{R}(A::MT, b::VT, C::MT, d::VT) where {R<:Real, MT<:AbstractMatrix{R}, VT<:AbstractVector{R}, QP<:QPDAS.QuadraticProgram}
+        @assert size(A,1) == size(b,1)
+        qp = QPDAS.QuadraticProgram(A, b, C, d, smartstart=false)
+        new{R, MT, VT, typeof(qp)}(A, b, C, d, zeros(R,size(A,2)), qp, Ref(true))
+    end
+end
+
+# properties
+
+is_prox_accurate(::IndPolyhedralQPDAS) = true
+
+# constructors
+
+function IndPolyhedralQPDAS(
+    l::AbstractVector{R}, A::AbstractMatrix{R}, u::AbstractVector{R}
+) where R
+    if !all(l .<= u)
+        error("function is improper (are some bounds inverted?)")
+    end
+    eqinds = (l .== u) .& .!isnothing.(l)
+    Aeq = A[eqinds,:]
+    beq = l[eqinds]
+
+    _islower(l::T) where T =
+        l != typemin(T) && !isnan(l) && !isnothing(l)
+    _isupper(u::T) where T =
+        u != typemax(T) && !isnan(u) && !isnothing(u)
+
+    lower = _islower.(l) .& (.!eqinds)
+    upper = _isupper.(u) .& (.!eqinds)
+
+    lower_only = lower .& (.! upper)
+    upper_only = upper .& (.! lower)
+    upper_and_lower = upper .& lower
+
+
+    Cieq = [-A[lower_only, :];
+             A[upper_only, :];
+            -A[upper_and_lower, :];
+             A[upper_and_lower, :] ]
+    dieq = [-l[lower_only];
+             u[upper_only];
+            -l[upper_and_lower];
+             u[upper_and_lower] ]
+
+    IndPolyhedralQPDAS{R}(Aeq, beq, Cieq, dieq)
+end
+
+IndPolyhedralQPDAS(
+    l::AbstractVector{R}, A::AbstractMatrix{R}, u::AbstractVector{R},
+    xmin::AbstractVector{R}, xmax::AbstractVector{R}
+) where R =
+    IndPolyhedralQPDAS([l; xmin], [A; I], [u; xmax])
+
+IndPolyhedralQPDAS(
+    l::AbstractVector{R}, A::AbstractMatrix{R}, args...
+) where R =
+    IndPolyhedralQPDAS(
+        l, A, R(Inf).*ones(R, size(A, 1)), args...
+    )
+
+IndPolyhedralQPDAS(
+    A::AbstractMatrix{R}, u::AbstractVector{R}, args...
+) where R =
+    IndPolyhedralQPDAS(
+        R(-Inf).*ones(R, size(A, 1)), A, u, args...
+    )
+
+# function evaluation
+
+function (f::IndPolyhedralQPDAS{R})(x::AbstractVector{R}) where R
+    Ax = f.A * x
+    Cx = f.C * x
+    return all(Ax .<= f.b .& Cx .<= f.d) ? R(0) : Inf
+end
+
+# prox
+
+function prox!(y::AbstractVector{R}, f::IndPolyhedralQPDAS{R}, x::AbstractVector{R}, gamma::R=R(1)) where R
+    # Linear term in qp is -x
+    f.z .= .- x
+    # Update the problem
+    QPDAS.update!(f.qp, z=f.z)
+
+    if f.first_prox[]
+        # This sets the initial active set based on z, should only be run once
+        QPDAS.run_smartstart(f.qp.boxQP)
+        f.first_prox[] = false
+    end
+    sol, val = QPDAS.solve!(f.qp)
+    y .= sol
+    return R(0)
+end
+
+# naive prox
+
+# we want to compute the projection p of a point x
+#
+# primal problem is: minimize_p (1/2)||p-x||^2 + g(Ap)
+# where g is the indicator of the box [l, u]
+#
+# dual problem is: minimize_y (1/2)||-A'y||^2 - x'A'y + g*(y)
+# can solve with (fast) dual proximal gradient method
+
+function prox_naive(f::IndPolyhedralQPDAS{R}, x::AbstractVector{R}, gamma::R=R(1)) where R
+    # Rewrite to l ≤ Ax ≤ u
+    A = [f.A; f.C]
+    l = [f.b; fill(R(-Inf), length(f.d))]
+    u = [f.b; f.d]
+    y = zeros(R, size(A, 1)) # dual vector
+    y1 = y
+    g = IndBox(l, u)
+    gstar = Conjugate(g)
+    gstar_y = R(0)
+    stepsize = R(1)/opnorm(Matrix(A*A'))
+    for it = 1:1e6
+        w = y + (it-1)/(it+2)*(y - y1)
+        y1 = y
+        z = w - stepsize * (A * (A'*w - x))
+        y, = prox(gstar, z, stepsize)
+        if norm(y-w)/(1+norm(w)) <= 1e-12 break end
+    end
+    p = -A'*y + x
+    return p, R(0)
+end

test/benchmarkpolyhedral.jl

@@ -0,0 +1,83 @@
+using ProximalOperators, Random
+# Number of variables
+n = 1000
+# Number of halfspaces
+mi = 50 # Inequalities with C
+me = 50 # Equalities with A
+
+Random.seed!(1)
+# One point in polytope
+x0 = randn(n)
+
+# Create polytope containing x0
+# Inequality
+C = Matrix{Float64}(undef, mi, n)
+d = randn(mi)
+
+# Make sure x0 is in polytope by setting sign of inequality, random C part
+for i = 1:mi
+    v = randn(n)
+    b = randn()
+    if v'x0  <= b
+        C[i,:] .= v
+    else
+        C[i,:] .= -v
+    end
+    d[i] = b
+end
+
+# Create equality
+A = randn(me, n)
+b = A*x0
+
+l = [b;fill(-Inf, mi)]
+u = [b;d]
+AC = [A;C]
+
+# Precompile
+polyOSQP  = IndPolyhedral(l, AC, u; solver=:osqp)
+polyQPDAS = IndPolyhedral(l, AC, u; solver=:qpdas)
+x = randn(n)
+y = similar(x0)
+prox!(y, polyOSQP, x)
+prox!(y, polyQPDAS, x)
+# Run tests
+println("Setup OSQP")
+@time polyOSQP  = IndPolyhedral(l, AC, u; solver=:osqp)
+println("Setup QPDAS")
+@time polyQPDAS = IndPolyhedral(l, AC, u; solver=:qpdas)
+
+Random.seed!(2)
+x = randn(n)
+y = similar(x0)
+println("First prox OSQP:")
+@time prox!(y, polyOSQP, x)
+
+Random.seed!(2)
+x = randn(n)
+println("First prox QPDAS:")
+@time prox!(y, polyQPDAS, x)
+
+Random.seed!(3)
+N = 100
+xs = randn(n, N)
+x = xs[:,1]
+
+println("100 prox! OSQP:")
+@time for i = 1:100
+    # Project from here
+    x .= xs[:,i]
+    prox!(y, polyOSQP, x)
+end
+
+Random.seed!(3)
+N = 100
+xs = randn(n, N)
+x = xs[:,1]
+
+println("100 prox! QPDAS:")
+@time for i = 1:100
+    # Project from here
+    x .= xs[:,i]
+    prox!(y, polyQPDAS, x)
+end

test/test_indPolyhedral.jl

@@ -23,32 +23,36 @@ p = similar(x)
 # define test cases
 
 constructors_positive = [
-    () -> IndPolyhedral(l, A),
-    () -> IndPolyhedral(l, A, xmin, xmax),
-    () -> IndPolyhedral(A, u),
-    () -> IndPolyhedral(A, u, xmin, xmax),
-    () -> IndPolyhedral(l, A, u),
-    () -> IndPolyhedral(l, A, u, xmin, xmax),
+    (solver) -> IndPolyhedral(l, A, solver=solver),
+    (solver) -> IndPolyhedral(l, A, xmin, xmax, solver=solver),
+    (solver) -> IndPolyhedral(A, u, solver=solver),
+    (solver) -> IndPolyhedral(A, u, xmin, xmax, solver=solver),
+    (solver) -> IndPolyhedral(l, A, u, solver=solver),
+    (solver) -> IndPolyhedral(l, A, u, xmin, xmax, solver=solver),
 ]
 
 constructors_negative = [
-    () -> IndPolyhedral(l, A, xmax, xmin),
-    () -> IndPolyhedral(A, u, xmax, xmin),
-    () -> IndPolyhedral(l, A, u, xmax, xmin),
+    (solver) -> IndPolyhedral(l, A, xmax, xmin, solver=solver),
+    (solver) -> IndPolyhedral(A, u, xmax, xmin, solver=solver),
+    (solver) -> IndPolyhedral(l, A, u, xmax, xmin, solver=solver),
 ]
 
 # run positive tests
 
 for constr in constructors_positive
-    f = constr()
-    @test ProximalOperators.is_convex(f) == true
-    @test ProximalOperators.is_set(f) == true
-    fx = call_test(f, x)
-    p, fp = prox_test(f, x)
+    for solver in [:osqp, :qpdas]
+        f = constr(solver)
+        @test ProximalOperators.is_convex(f) == true
+        @test ProximalOperators.is_set(f) == true
+        fx = call_test(f, x)
+        p, fp = prox_test(f, x)
+    end
 end
 
 # run negative tests
 
 for constr in constructors_negative
-    @test_throws ErrorException constr()
+    for solver in [:osqp, :qpdas]
+        @test_throws ErrorException constr(solver)
+    end
 end