-
Notifications
You must be signed in to change notification settings - Fork 4
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Exponential Cone #29
Exponential Cone #29
Changes from 19 commits
2fb3701
906e784
45bac58
f21cded
052dbe1
ba899b8
4bf77c4
3ce2c5e
7ead9ad
48d3fcf
52fd09f
1988359
cc830fc
2070ab1
b1e922e
773bf8c
82fc6f9
b1fa666
4e12347
a61e4ee
6c275a1
3912710
276ca15
7f37235
c97645a
25d8518
1b24cf6
f7688e1
9a1f470
dc2606a
b547f28
aee7112
edad4d4
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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,92 @@ 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) | ||
|
||
tjdiamandis marked this conversation as resolved.
Show resolved
Hide resolved
|
||
if distance_to_set(DefaultDistance(), v, MOI.ExponentialCone()) < 1e-8 | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I think distances should use projections more than the other way around. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. What's the advantage of using projections for distances? If we fall into case 4 for the exponential cone projection, i.e., (x,y,z) not in K, -K*, and (x > 0 or y > 0), then computing the projection is going to take significantly longer than computing the distance. A few microseconds vs 10s of nanoseconds. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
distance can be seen as a composition of the projection with a point-to-point distance. What we need here to be precise is the set membership for the exp cone and its polar. I advised using distance for this and checking distance approx 0. @joaquimg do you think we should use another membership function here? The motivation for re-using distance was avoiding duplication There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This might be the first version of having two distances. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. @tjdiamandis to simplify your work for now, you can just implement the quick check _in_exponential(v) that you had before. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. No worries! Sounds good. Will do today. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. We can just open an issue and let it be like this for now. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. see here: #32 |
||
return v | ||
end | ||
if distance_to_set(DefaultDistance(), -v, MOI.DualExponentialCone()) < 1e-8 | ||
# if in polar cone Ko = -K* | ||
return zeros(3) | ||
tjdiamandis marked this conversation as resolved.
Show resolved
Hide resolved
|
||
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 _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(::DefaultDistance, v::AbstractVector{T}, ::MOI.DualExponentialCone) where {T} | ||
return v + projection_on_set(DefaultDistance(), -v, MOI.ExponentialCone()) | ||
end | ||
tjdiamandis marked this conversation as resolved.
Show resolved
Hide resolved
|
||
|
||
""" | ||
projection_on_set(::DefaultDistance, v::AbstractVector{T}, sets::Array{<:MOI.AbstractSet}) | ||
|
||
|
@@ -295,6 +382,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) | ||
tjdiamandis marked this conversation as resolved.
Show resolved
Hide resolved
|
||
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) | ||
tjdiamandis marked this conversation as resolved.
Show resolved
Hide resolved
|
||
end | ||
if v[1] <= 0 && v[2] <= 0 | ||
return diagm([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}) | ||
|
||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,37 @@ | ||
|
||
function _bisection(f, left, right; max_iters=10000, tol=1e-10) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. The code in this function will yield many iterations or not converge because of strict inequality check, one should use: There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Originally designed to hit floating point limit (since intervals are relatively small), but I agree that this level of precision probably yields too many iterations. Will change to There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
yes but here you ask for strict equality right? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Yes, so the algorithm stops when the boundary is indistinguishable from the midpoint to floating point precision. E.g., |
||
# 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 |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Random note:
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Oh interesting. Admittedly, I didn't think to check the latter. |
||
@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 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Dont we need Random here?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I saw that Random was already in the package dependencies, but I don't think it's used outside of tests. I can move it to
extras
, unless there's something I'm missing.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'll do that in a following PR, we can leave it in the deps for now