From 9ef481fbf3fc231bdbab35f670f9a7ee8af81cf6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mattias=20F=C3=A4lt?= Date: Wed, 27 Mar 2019 20:49:28 +0100 Subject: [PATCH] First version of polyhedral with QPDAS --- Project.toml | 1 + REQUIRE | 1 + src/functions/indPolyhedral.jl | 3 + src/functions/indPolyhedralQPDAS.jl | 146 ++++++++++++++++++++++++++++ test/benchmarkpolyhedral.jl | 83 ++++++++++++++++ test/test_indPolyhedral.jl | 34 ++++--- 6 files changed, 253 insertions(+), 15 deletions(-) create mode 100644 src/functions/indPolyhedralQPDAS.jl create mode 100644 test/benchmarkpolyhedral.jl diff --git a/Project.toml b/Project.toml index 1f01fcab..a318f9a6 100644 --- a/Project.toml +++ b/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" diff --git a/REQUIRE b/REQUIRE index 79953b01..c8a0b5fc 100644 --- a/REQUIRE +++ b/REQUIRE @@ -2,3 +2,4 @@ julia 1.0 IterativeSolvers 0.8.0 TSVD 0.3.0 OSQP 0.3.0 +QPDAS 0.1.0 diff --git a/src/functions/indPolyhedral.jl b/src/functions/indPolyhedral.jl index 70fe2479..72b8d5d7 100644 --- a/src/functions/indPolyhedral.jl +++ b/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") diff --git a/src/functions/indPolyhedralQPDAS.jl b/src/functions/indPolyhedralQPDAS.jl new file mode 100644 index 00000000..5d4f0244 --- /dev/null +++ b/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 diff --git a/test/benchmarkpolyhedral.jl b/test/benchmarkpolyhedral.jl new file mode 100644 index 00000000..fabf96cf --- /dev/null +++ b/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 diff --git a/test/test_indPolyhedral.jl b/test/test_indPolyhedral.jl index 8656d981..a0bf3f10 100644 --- a/test/test_indPolyhedral.jl +++ b/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