Exponential Cone

Project.toml

@@ -23,7 +23,9 @@ julia = "1"
 
 [extras]
 ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "ChainRulesTestUtils"]
+test = ["Test", "ChainRulesTestUtils", "JuMP", "SCS"]

src/MathOptSetDistances.jl

@@ -11,6 +11,7 @@ import FiniteDifferences
 
 export distance_to_set, projection_on_set, projection_gradient_on_set
 
+include("utils.jl")
 include("distances.jl")
 include("distance_sets.jl")
 include("projections.jl")

src/distance_sets.jl

@@ -131,6 +131,11 @@ function distance_to_set(::NormedEpigraphDistance{p}, v::AbstractVector{<:Real},
     x = v[1]
     y = v[2]
     z = v[3]
+
+    if x <= 0 && isapprox(y, 0, atol=1e-10) && z >= 0
+        return zero(eltype(v))
+    end
+
     result = y * exp(x/y) - z
     return LinearAlgebra.norm(
         (max(-y, zero(result)), max(result, zero(result))),
@@ -147,6 +152,11 @@ function distance_to_set(::NormedEpigraphDistance{p}, vs::AbstractVector{<:Real}
     u = vs[1]
     v = vs[2]
     w = vs[3]
+
+    if isapprox(u, 0, atol=1e-10) && v >= 0 && w >= 0
+        return zero(result)
+    end
+
     result = -u*exp(v/u) - ℯ * w
     return LinearAlgebra.norm(
         (max(u, zero(result)), max(result, zero(result))),

src/projections.jl

@@ -1,6 +1,7 @@
 # find expression of projections on cones and their derivatives here:
 #   https://stanford.edu/~boyd/papers/pdf/cone_prog_refine.pdf
 
+
 """
     projection_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.Zeros) where {T}
 
@@ -152,6 +153,106 @@ function vec_symm(X)
     return X[LinearAlgebra.triu(trues(size(X)))]
 end
 
+"""
+    projection_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.ExponentialCone) where {T}
+
+projection of vector `v` on closure of the exponential cone
+i.e. `cl(Kexp) = {(x,y,z) | y e^(x/y) <= z, y>0 } U {(x,y,z)| x <= 0, y = 0, z >= 0}`.
+
+References:
+* [Proximal Algorithms, 6.3.4](https://web.stanford.edu/~boyd/papers/pdf/prox_algs.pdf)
+by Neal Parikh and Stephen Boyd.
+* [Projection, presolve in MOSEK: exponential, and power cones]
+(https://docs.mosek.com/slides/2018/ismp2018/ismp-friberg.pdf) by Henrik Friberg
+"""
+function projection_on_set(::DefaultDistance, v::AbstractVector{T}, s::MOI.ExponentialCone) where {T}
+    _check_dimension(v, s)
+
+    if _in_exp_cone(v; dual=false)
+        return v
+    end
+    if _in_exp_cone(-v; dual=true)
+        # if in polar cone Ko = -K*
+        return zeros(T, 3)
+    end
+    if v[1] <= 0 && v[2] <= 0
+        return [v[1]; 0.0; max(v[3],0)]
+    end
+
+    return _exp_cone_proj_case_4(v)
+end
+
+function _in_exp_cone(v::AbstractVector{T}; dual=false, tol=1e-8) where {T}
+    if dual
+        return (
+            (isapprox(v[1], 0, atol=tol) && v[2] >= 0 && v[3] >= 0) ||
+            (v[1] < 0 && v[1]*exp(v[2]/v[1]) + ℯ*v[3] >= tol)
+        )
+    else
+        return (
+            (v[1] <= 0 && isapprox(v[2], 0, atol=tol) && v[3] >= 0) ||
+            (v[2] > 0 && v[2] * exp(v[1] / v[2]) - v[3] <= tol)
+        )
+    end
+end
+
+function _exp_cone_proj_case_4(v::AbstractVector{T}) where {T}
+    # Ref: https://docs.mosek.com/slides/2018/ismp2018/ismp-friberg.pdf, p47-48
+    # Thm: h(x) is smooth, strictly increasing, and changes sign on domain
+    r, s, t = v[1], v[2], v[3]
+    h(x) = (((x-1)*r + s) * exp(x) - (r - x*s)*exp(-x))/(x^2 - x + 1) - t
+
+    # Note: won't both be Inf by case 3 of projection
+    lb = r > 0 ? 1 - s/r : -Inf
+    ub = s > 0 ? r/s : Inf
+
+    # Deal with ±Inf bounds
+    if isinf(lb)
+        lb = min(ub-0.125, -0.125)
+        for _ in 1:10
+            h(lb) < 0 && break
+            ub = lb
+            lb *= 2
+        end
+    end
+    if isinf(ub)
+        ub = max(lb+0.125, 0.125)
+        for _ in 1:10
+            h(ub) > 0 && break
+            lb = ub
+            ub *= 2
+        end
+    end
+
+    # Check bounds
+    if !(h(lb) < 0 && h(ub) > 0)
+        error("Failure to find bracketing interval for exp cone projection.")
+    end
+
+    x, code = _bisection(h, lb, ub)
+    if code == 2
+        error("Failured to solve root finding problem in exp cone projection")
+    end
+
+    return ((x - 1)*r + s)/(x^2 - x + 1) * [x; 1.0; exp(x)]
+end
+
+"""
+    projection_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.DualExponentialCone) where {T}
+
+projection of vector `v` on the dual exponential cone
+i.e. `Kexp^* = {(u,v,w) | u < 0, -u*exp(v/u) <= ew } U {(u,v,w)| u == 0, v >= 0, w >= 0}`.
+
+References:
+* [Proximal Algorithms, 6.3.4](https://web.stanford.edu/~boyd/papers/pdf/prox_algs.pdf)
+by Neal Parikh and Stephen Boyd.
+* [Projection, presolve in MOSEK: exponential, and power cones](https://docs.mosek.com/slides/2018/ismp2018/ismp-friberg.pdf)
+by Henrik Friberg
+"""
+function projection_on_set(d::DefaultDistance, v::AbstractVector{T}, ::MOI.DualExponentialCone) where {T}
+    return v + projection_on_set(d, -v, MOI.ExponentialCone())
+end
+
 """
     projection_on_set(::DefaultDistance, v::AbstractVector{T}, sets::Array{<:MOI.AbstractSet})
 
@@ -295,6 +396,61 @@ function projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, ::M
     return D
 end
 
+"""
+    projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.ExponentialCone) where {T}
+
+derivative of projection of vector `v` on closure of the exponential cone,
+i.e. `cl(Kexp) = {(x,y,z) | y e^(x/y) <= z, y>0 } U {(x,y,z)| x <= 0, y = 0, z >= 0}`.
+
+References:
+* [Solution Refinement at Regular Points of Conic Problems](https://stanford.edu/~boyd/papers/cone_prog_refine.html)
+by Enzo Busseti, Walaa M. Moursi, and Stephen Boyd
+"""
+function projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, s::MOI.ExponentialCone) where {T}
+    _check_dimension(v, s)
+    Ip(z) = z >= 0 ? 1 : 0
+
+    if distance_to_set(DefaultDistance(), v, MOI.ExponentialCone()) < 1e-8
+        return Matrix{Float64}(I, 3, 3)
+    end
+    if distance_to_set(DefaultDistance(), -v, MOI.DualExponentialCone()) < 1e-8
+        # if in polar cone Ko = -K*
+        return zeros(3,3) #FillArrays.Zeros(3, 3)
+    end
+    if v[1] <= 0 && v[2] <= 0
+        return Diagonal([1; Ip(v[2]); Ip(v[3])])
+    end
+
+    z1, z2, z3 = _exp_cone_proj_case_4(v)
+    nu = z3 - v[3]
+    rs = z1/z2
+    exp_rs = exp(rs)
+
+    mat = inv([
+        1+nu*exp_rs/z2     -nu*exp_rs*rs/z2       0     exp_rs;
+        -nu*exp_rs*rs/z2   1+nu*exp_rs*rs^2/z2    0     (1-rs)*exp_rs;
+        0                  0                      1     -1
+        exp_rs             (1-rs)*exp_rs          -1    0
+    ])
+    return @view(mat[1:3,1:3])
+end
+
+"""
+    projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.DualExponentialCone) where {T}
+
+derivative of projection of vector `v` on the dual exponential cone,
+i.e. `Kexp^* = {(u,v,w) | u < 0, -u*exp(v/u) <= ew } U {(u,v,w)| u == 0, v >= 0, w >= 0}`.
+
+References:
+* [Solution Refinement at Regular Points of Conic Problems]
+(https://stanford.edu/~boyd/papers/cone_prog_refine.html)
+by Enzo Busseti, Walaa M. Moursi, and Stephen Boyd
+"""
+function projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.DualExponentialCone) where {T}
+    # from Moreau decomposition: x = P_K(x) + P_-K*(x)
+    return I - projection_gradient_on_set(DefaultDistance(), -v, MOI.ExponentialCone())
+end
+
 """
     projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, sets::Array{<:MOI.AbstractSet})
 

src/utils.jl

@@ -0,0 +1,37 @@
+
+function _bisection(f, left, right; max_iters=10000, tol=1e-10)
+    # STOP CODES:
+    #   0: Success (floating point limit or exactly 0)
+    #   1: Max iters but within tol
+    #   2: Failure
+
+    for _ in 1:max_iters
+        f_left, f_right = f(left), f(right)
+        sign(f_left) == sign(f_right) && error("Interval became non-bracketing.")
+
+        mid = (left + right) / 2
+        if left == mid || right == mid
+            return mid, 0
+        end
+
+        f_mid = f(mid)
+        if f_mid == 0
+            return mid, 0
+        end
+        if sign(f_mid) == sign(f_left)
+            left = mid
+            continue
+        end
+        if  sign(f_mid) == sign(f_right)
+            right = mid
+            continue
+        end
+    end
+
+    mid = (left + right) / 2
+    if abs(f(mid)) < tol
+        return mid, 1
+    end
+
+    return nothing, 2
+end

test/projection_gradients.jl

@@ -169,4 +169,59 @@ end
             @test ≈(dfor, dΠ, atol=1e-5) || ≈(dback, dΠ, atol=1e-5)
         end
     end
+
+    @testset "Exp Cone" begin
+        function det_case_exp_cone(v; dual=false)
+            v = dual ? -v : v
+            if MOD.distance_to_set(DD, v, MOI.ExponentialCone()) < 1e-8
+                return 1
+            elseif MOD.distance_to_set(DD, -v, MOI.DualExponentialCone()) < 1e-8
+                return 2
+            elseif v[1] <= 0 && v[2] <= 0 #TODO: threshold here??
+                return 3
+            else
+                return 4
+            end
+        end
+
+        Random.seed!(0)
+        s = MOI.ExponentialCone()
+        sd = MOI.DualExponentialCone()
+        case_p = zeros(4)
+        case_d = zeros(4)
+        # Adjust tolerance down because a 1-2 errors when projection ends up
+        #   very close to the z axis
+        # For intuition, see Fig 5.1 https://docs.mosek.com/modeling-cookbook/expo.html
+        #   Note that their order is reversed: (x, y, z) = (x3, x2, x1) [theirs]
+        tol = 1e-6
+        for ii in 1:100
+            v = 5*randn(3)
+            @testset "Primal Cone" begin
+                case_p[det_case_exp_cone(v; dual=false)] += 1
+                dΠ = MOD.projection_gradient_on_set(MOD.DefaultDistance(), v, s)
+                grad_fdm1 = FiniteDifferences.jacobian(ffdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
+                grad_fdm2 = FiniteDifferences.jacobian(bfdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
+                @test size(grad_fdm1) == size(grad_fdm2) == size(dΠ)
+                @test ≈(dΠ, grad_fdm1,atol=tol) || ≈(dΠ, grad_fdm2, atol=tol)
+                if !(≈(dΠ, grad_fdm1,atol=tol) || ≈(dΠ, grad_fdm2, atol=tol))
+                    println("v = $v")
+                    println("n1: $(norm(dΠ - grad_fdm1))\nn2: $(norm(dΠ - grad_fdm2))")
+                end
+            end
+
+            @testset "Dual Cone" begin
+                case_d[det_case_exp_cone(v; dual=true)] += 1
+                dΠ = MOD.projection_gradient_on_set(MOD.DefaultDistance(), v, sd)
+                grad_fdm1 = FiniteDifferences.jacobian(ffdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, sd), v)[1]'
+                grad_fdm2 = FiniteDifferences.jacobian(bfdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, sd), v)[1]'
+                @test size(grad_fdm1) == size(grad_fdm2) == size(dΠ)
+                @test ≈(dΠ, grad_fdm1,atol=tol) || ≈(dΠ, grad_fdm2, atol=tol)
+                if !(≈(dΠ, grad_fdm1,atol=tol) || ≈(dΠ, grad_fdm2, atol=tol))
+                    println("v = $v")
+                    println("n1: $(norm(dΠ - grad_fdm1))\nn2: $(norm(dΠ - grad_fdm2))")
+                end
+            end
+        end
+        @test all(case_p .> 0) && all(case_d .> 0)
+    end
 end

test/projections.jl

@@ -1,3 +1,4 @@
+using JuMP, SCS
 const DD = MOD.DefaultDistance()
 
 @testset "Test projections distance on vector sets" begin
@@ -76,3 +77,56 @@ end
     output_joint = MOD.projection_gradient_on_set(DD, [v1, v2], [c1, c2])
     @test output_joint ≈ BlockDiagonal([output_1, output_2])
 end
+
+
+@testset "Exponential Cone Projections" begin
+    function det_case_exp_cone(v; dual=false)
+        v = dual ? -v : v
+        if MOD.distance_to_set(DD, v, MOI.ExponentialCone()) < 1e-8
+            return 1
+        elseif MOD.distance_to_set(DD, -v, MOI.DualExponentialCone()) < 1e-8
+            return 2
+        elseif v[1] <= 0 && v[2] <= 0 #TODO: threshold here??
+            return 3
+        else
+            return 4
+        end
+    end
+
+    function _test_proj_exp_cone_help(x, tol; dual=false)
+        cone = dual ? MOI.DualExponentialCone() : MOI.ExponentialCone()
+        model = Model()
+        set_optimizer(model, optimizer_with_attributes(
+            SCS.Optimizer, "eps" => 1e-10, "max_iters" => 10000, "verbose" => 0))
+        @variable(model, z[1:3])
+        @variable(model, t)
+        @objective(model, Min, t)
+        @constraint(model, sum((x-z).^2) <= t)
+        @constraint(model, z in cone)
+        optimize!(model)
+        z_star = value.(z)
+        px = MOD.projection_on_set(DD, x, cone)
+        if !isapprox(px, z_star, atol=tol)
+            # error("Exp cone projection failed:\n x = $x\nMOD: $px\nJuMP: $z_star
+            #        \nnorm: $(norm(px - z_star))")
+            return false
+       end
+       return true
+    end
+
+    Random.seed!(0)
+    n = 3
+    atol = 1e-7
+    case_p = zeros(4)
+    case_d = zeros(4)
+    for _ in 1:100
+        x = randn(3)
+
+        case_p[det_case_exp_cone(x; dual=false)] += 1
+        @test _test_proj_exp_cone_help(x, atol; dual=false)
+
+        case_d[det_case_exp_cone(x; dual=true)] += 1
+        @test _test_proj_exp_cone_help(x, atol; dual=true)
+    end
+    @test all(case_p .> 0) && all(case_d .> 0)
+end

test/runtests.jl

@@ -6,6 +6,7 @@ const MOD = MathOptSetDistances
 const MOI = MathOptSetDistances.MOI
 
 using LinearAlgebra
+using Random
 import BlockDiagonals: BlockDiagonal
 
 include("distances.jl")