Exponential Cone

Project.toml

@@ -10,6 +10,7 @@ FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
@@ -23,7 +24,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

@@ -8,6 +8,7 @@ using LinearAlgebra
 import ChainRulesCore
 const CRC = ChainRulesCore
 import FiniteDifferences
+import NonlinearSolve
 
 export distance_to_set, projection_on_set, projection_gradient_on_set
 

src/projections.jl

@@ -1,5 +1,11 @@
 # find expression of projections on cones and their derivatives here:
 #   https://stanford.edu/~boyd/papers/pdf/cone_prog_refine.pdf
+# See also
+#   https://github.com/tjdiamandis/ConeProgramDiff.jl/blob/main/cone_ref.pdf
+#   and references therein
+const EXP_CONE_THRESH = 1e-8
+const POW_CONE_THRESH = 1e-8
+
 
 """
     projection_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.Zeros) where {T}
@@ -152,6 +158,116 @@ 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
+    elseif _in_exp_cone(-v; dual=true)
+        # if in polar cone Ko = -K*
+        return zeros(3)
+    elseif v[1] <= 0 && v[2] <= 0
+        return [v[1]; 0.0; max(v[3],0)]
+    else
+        return _exp_cone_proj_case_4(v)
+    end
+end
+
+function _in_exp_cone(v::AbstractVector{T}; dual=false) where {T}
+    # See pg. 184 https://web.stanford.edu/~boyd/papers/pdf/prox_algs.pdf
+    # TODO: Tol for == 0 to avoid denom blowing up? see case 4 in deriv
+    if dual
+        return (
+            (v[1] == 0 && v[2] >= 0 && v[3] >= 0) ||
+            (v[1] < 0 && v[1]*exp(v[2]/v[1]) + ℯ*v[3] >= EXP_CONE_THRESH)
+        )
+    else
+        return (
+            (v[1] <= 0 && v[2] == 0 && v[3] >= 0) ||
+            (v[2] > 0 && v[2] * exp(v[1] / v[2]) - v[3] <= EXP_CONE_THRESH)
+        )
+    end
+end
+
+function _exp_cone_proj_case_4(v::AbstractVector{T}) where {T}
+    # https://docs.mosek.com/slides/2018/ismp2018/ismp-friberg.pdf
+    # Thm: h(x) is smooth, strictly increasing, and changes sign on domain
+    r, s, t = v[1], v[2], v[3]
+    h(x,p) = (((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, nothing) < 0 && break
+            ub = lb
+            lb *= 2
+        end
+    elseif isinf(ub)
+        ub = max(lb+0.125, 0.125)
+        for _ in 1:10
+            h(ub, nothing) > 0 && break
+            lb = ub
+            ub *= 2
+        end
+    end
+
+    if !(h(lb, nothing) < 0 && h(ub, nothing) > 0)
+        error("Failure to find bracketing interval for exp cone")
+    end
+
+    # Solve with Bisection
+    prob = NonlinearSolve.NonlinearProblem(h, (lb, ub))
+    sol = NonlinearSolve.solve(prob, NonlinearSolve.Bisection())
+
+    if sol.retcode == NonlinearSolve.MAXITERS_EXCEED
+        error("Numerical error in exp cone projection")
+    elseif sol.retcode == NonlinearSolve.FLOATING_POINT_LIMIT
+        # left == mid or right == mid, and (left, right) still bracketing
+        x = (sol.left + sol.right)/2
+    elseif sol.retcode == NonlinearSolve.EXACT_SOLUTION_LEFT
+        x = sol.left
+    elseif sol.retcode == NonlinearSolve.EXACT_SOLUTION_RIGHT
+        x = sol.right
+    else
+        error("NonlinearSolve.jl retcode not recognized in exp cone")
+    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(::DefaultDistance, v::AbstractVector{T}, ::MOI.DualExponentialCone) where {T}
+    return v + projection_on_set(DefaultDistance(), -v, MOI.ExponentialCone())
+end
+
 """
     projection_on_set(::DefaultDistance, v::AbstractVector{T}, sets::Array{<:MOI.AbstractSet})
 
@@ -295,6 +411,59 @@ 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 _in_exp_cone(v; dual=false)
+        return Matrix{Float64}(I, 3, 3)
+    elseif _in_exp_cone(-v; dual=true)
+        # if in polar cone Ko = -K*
+        return zeros(3,3) #FillArrays.Zeros(3, 3)
+    elseif v[1] <= 0 && v[2] <= 0
+        return diagm([1; Ip(v[2]); Ip(v[3])])
+    else
+        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
+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})
 

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._in_exp_cone(v; dual=false)
+                return 1
+            elseif MOD._in_exp_cone(-v; dual=true)
+                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._in_exp_cone(v; dual=false)
+            return 1
+        elseif MOD._in_exp_cone(-v; dual=true)
+            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")