Skip to content
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

Merged
merged 33 commits into from
Dec 27, 2020
Merged

Exponential Cone #29

merged 33 commits into from
Dec 27, 2020

Conversation

tjdiamandis
Copy link
Contributor

Projection onto exponential cone and its derivative.

The function _in_exp_cone(v::AbstractVector{T}; dual=false) may be necessary and incorporated into the distance function. However, in the projection I allow some tolerance for the points close to the boundary.

@codecov
Copy link

codecov bot commented Dec 23, 2020

Codecov Report

Merging #29 (edad4d4) into master (3af4edb) will increase coverage by 1.10%.
The diff coverage is 93.33%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master      #29      +/-   ##
==========================================
+ Coverage   88.16%   89.26%   +1.10%     
==========================================
  Files           5        6       +1     
  Lines         355      438      +83     
==========================================
+ Hits          313      391      +78     
- Misses         42       47       +5     
Impacted Files Coverage Δ
src/MathOptSetDistances.jl 100.00% <ø> (ø)
src/distance_sets.jl 72.25% <50.00%> (-0.59%) ⬇️
src/utils.jl 94.11% <94.11%> (ø)
src/projections.jl 98.78% <96.29%> (-1.22%) ⬇️
src/chainrules.jl 100.00% <0.00%> (ø)

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 3af4edb...edad4d4. Read the comment docs.

src/projections.jl Outdated Show resolved Hide resolved
src/projections.jl Outdated Show resolved Hide resolved
tjdiamandis and others added 2 commits December 23, 2020 15:29
Co-authored-by: Mathieu Besançon <mathieu.besancon@gmail.com>
Co-authored-by: Mathieu Besançon <mathieu.besancon@gmail.com>
Project.toml Outdated Show resolved Hide resolved
src/projections.jl Outdated Show resolved Hide resolved
@matbesancon
Copy link
Owner

A major drawback is the addition of a dependency here on NonlinearSolve. Ideally this package will be used by others requiring set distances and projections so we should avoid bringing in a lot of direct and transitive dependencies)

src/projections.jl Outdated Show resolved Hide resolved
@tjdiamandis tjdiamandis marked this pull request as draft December 23, 2020 21:26
src/projections.jl Outdated Show resolved Hide resolved
Copy link
Collaborator

@joaquimg joaquimg left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great job!

function projection_on_set(::DefaultDistance, v::AbstractVector{T}, s::MOI.ExponentialCone) where {T}
_check_dimension(v, s)

if distance_to_set(DefaultDistance(), v, MOI.ExponentialCone()) < 1e-8
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think distances should use projections more than the other way around.

Copy link
Contributor Author

Choose a reason for hiding this comment

The 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.

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What's the advantage of using projections for distances?

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

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This might be the first version of having two distances.
Because a very common distance is the euclidean.
But euclidean requires the projection first.
Our current distance is some epigraphical distance which is good for boolean checking if the point is the the cone.
I don't we should remove the current distance, I would just rename them and have both.
One could be used inside the the other, but I would use euclidean by default. Actually, whenever I possible I think we should use euclidean to have is as uniform as possible.

Copy link
Owner

Choose a reason for hiding this comment

The 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.
Sorry for making you backtrack on this. I'll refactor it with the additional distance in another PR

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No worries! Sounds good. Will do today.

Copy link
Collaborator

Choose a reason for hiding this comment

The 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.

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

see here: #32

@@ -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"
Copy link
Collaborator

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?

Copy link
Contributor Author

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.

Copy link
Owner

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

src/projections.jl Outdated Show resolved Hide resolved
src/projections.jl Outdated Show resolved Hide resolved
src/utils.jl Outdated
@@ -0,0 +1,37 @@

function _bisection(f, left, right; max_iters=10000, tol=1e-10)
Copy link
Owner

Choose a reason for hiding this comment

The 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:
a ≈ b atol=tol instead of equality checks for floating-point numbers

Copy link
Contributor Author

@tjdiamandis tjdiamandis Dec 25, 2020

Choose a reason for hiding this comment

The 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 abs(f(mid)) < tol for early stop (similar to lines 32-33 for max_iter stop check).

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Originally designed to hit floating point limit

yes but here you ask for strict equality right?

Copy link
Contributor Author

Choose a reason for hiding this comment

The 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., 1 == 1 + eps()/2 evaluates to true.

tjdiamandis and others added 4 commits December 25, 2020 15:44
Co-authored-by: Mathieu Besançon <mathieu.besancon@gmail.com>
Co-authored-by: Mathieu Besançon <mathieu.besancon@gmail.com>
Co-authored-by: Mathieu Besançon <mathieu.besancon@gmail.com>
Co-authored-by: Mathieu Besançon <mathieu.besancon@gmail.com>
@variable(model, z[1:3])
@variable(model, t)
@objective(model, Min, t)
@constraint(model, sum((x-z).^2) <= t)
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Random note:
interesting, I didn't think this way of writing the squared norm would be faster than dot:

julia> @btime @constraint($model, sum(($x-$z).^2) <= $t)
  4.522 μs (112 allocations: 9.39 KiB)
julia> @btime @constraint($model, ($x-$z) ⋅ ($x-$z)  <= $t)
  6.419 μs (168 allocations: 13.28 KiB)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh interesting. Admittedly, I didn't think to check the latter.

src/projections.jl Outdated Show resolved Hide resolved
src/projections.jl Outdated Show resolved Hide resolved
src/projections.jl Outdated Show resolved Hide resolved
tjdiamandis and others added 3 commits December 25, 2020 17:10
Co-authored-by: Mathieu Besançon <mathieu.besancon@gmail.com>
Co-authored-by: Mathieu Besançon <mathieu.besancon@gmail.com>
Co-authored-by: Mathieu Besançon <mathieu.besancon@gmail.com>
@tjdiamandis
Copy link
Contributor Author

Sorry for the type return errors. Thanks so much for the thorough review! I understand the general design principle now.

@matbesancon
Copy link
Owner

no problem, thanks a lot for the work. This is also me failing to catch everything at once, hence the back and forth which can feel tedious.

For diagm, the issue is that Julia 1.0 did not have the default diagm(v) = diagm(0 => v) yet, with 0 the convention to indicate the main diagonal of the matrix

@matbesancon
Copy link
Owner

Tests look green, we'll let the PR there for a day if other folks want to take a look or add something, but looks good on my side!

src/utils.jl Outdated Show resolved Hide resolved
tjdiamandis and others added 3 commits December 25, 2020 17:43
Co-authored-by: Mathieu Besançon <mathieu.besancon@gmail.com>
Merge branch 'master' of github.com:tjdiamandis/MathOptSetDistances.jl into master
@tjdiamandis tjdiamandis marked this pull request as ready for review December 25, 2020 23:38
@matbesancon matbesancon merged commit 133212c into matbesancon:master Dec 27, 2020
@matbesancon
Copy link
Owner

Thanks again for the PR!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants