Adding (some) conic constraints.

[dnguyen227]

Adding some conic constraints that can be handled with BigM and Hull. Some things I found odd that maybe you can comment on or correct me.

1. Hull reformulation is standard throughout all cones, but cones in MOI inherit from AbstractVectorSet, so I can't just ask for any cone. Right now I'm stuck naming each type of cone rather than a generic MOI.Cone datatype.

2. BigM reformulation looks like it needs to be adapted per cone type, and right now I'm kind of having to reason out the correct direction based on the table provided in the paper. Snippet from the paper below.

Overall a good introduction to cones but I just can't get past these two things while I was trying to generalize the code.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 99.71%. Comparing base (471b92a) to head (55fa83f).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #136      +/-   ##
==========================================
+ Coverage   99.69%   99.71%   +0.01%     
==========================================
  Files          21       21              
  Lines        2298     2456     +158     
==========================================
+ Hits         2291     2449     +158     
  Misses          7        7              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[pulsipher]

Ya, we'll have to simply enumerate the conic sets. You could could make your own internal Union type though:

_ConicSets =Union{ExponentialCone, GeometricMeanCone, ...}

[pulsipher]

Thoughts @bernalde?

[bernalde]

Folks, this is exciting. I’ll do a more in depth review once you consider it useful, but one observation from my paper was the fact that the same cone appears in the BigM and HR reformulations, just involving the binary variables. If you already have access to the type of cone and how it was built, then you have everything to complete the recipe. Happy to talk more about this and provide code that we used (the repo in the paper should have all). Cheers!

[dnguyen227]

I had a look at using MOI for the hessian evaluation and in my opinion because it's quadratic we actually have simpler version since I can read off the coefficents.

What we do now. I can just read off whats given.

Q = zeros(n, n)
for (pair, c) in quad.terms
    i, j = index[pair.a], index[pair.b]
    i == j ? (Q[i, i] += c) : (Q[i, j] += c / 2; Q[j, i] += c / 2)
end
eigmin(Symmetric(Q)) >= -tol

For MOI, I end up having to build something, then evaluate (this version works locally, so I can swap to this if we'd prefer).

nl = MOI.Nonlinear.Model()
MOI.Nonlinear.add_constraint(nl, expr, MOI.LessThan(0.0))  # expr rebuilt in MOI's index space
ev = MOI.Nonlinear.Evaluator(nl, MOI.Nonlinear.SparseReverseMode(),  MOI.VariableIndex.(1:n))
MOI.initialize(ev, [:Hess])
rows = MOI.hessian_lagrangian_structure(ev)
H = zeros(length(rows))
MOI.eval_hessian_lagrangian(ev, H, zeros(n), 0.0, [1.0])
Q = zeros(n, n)
for (k, (i, j)) in enumerate(rows)
    Q[i, j] += H[k] / 2
    i == j || (Q[j, i] += H[k] / 2)
end
eigmin(Symmetric(Q)) >= -tol   # same test either way

[pulsipher]

I had a look at using MOI for the hessian evaluation and in my opinion because it's quadratic we actually have simpler version since I can read off the coefficents.

What we do now. I can just read off whats given.

Q = zeros(n, n)
for (pair, c) in quad.terms
    i, j = index[pair.a], index[pair.b]
    i == j ? (Q[i, i] += c) : (Q[i, j] += c / 2; Q[j, i] += c / 2)
end
eigmin(Symmetric(Q)) >= -tol

For MOI, I end up having to build something, then evaluate (this version works locally, so I can swap to this if we'd prefer).

nl = MOI.Nonlinear.Model()
MOI.Nonlinear.add_constraint(nl, expr, MOI.LessThan(0.0))  # expr rebuilt in MOI's index space
ev = MOI.Nonlinear.Evaluator(nl, MOI.Nonlinear.SparseReverseMode(),  MOI.VariableIndex.(1:n))
MOI.initialize(ev, [:Hess])
rows = MOI.hessian_lagrangian_structure(ev)
H = zeros(length(rows))
MOI.eval_hessian_lagrangian(ev, H, zeros(n), 0.0, [1.0])
Q = zeros(n, n)
for (k, (i, j)) in enumerate(rows)
    Q[i, j] += H[k] / 2
    i == j || (Q[j, i] += H[k] / 2)
end
eigmin(Symmetric(Q)) >= -tol   # same test either way

I would say that it is preferred to directly compute it if you can easily do so based on the type. If it made since, we could make the AD approach a fallback, but I believe you have to enumerate each conic type anyway in MOI.

[dnguyen227]

@pulsipher ready for review

[pulsipher]

Looks good to me, I just would like to see the docs reflect that these capabilities have been added.

@bernalde: Do you approve?

[bernalde]

Review of 6e87fa2 (7 commits on top of a4ebe8a; local diff file set matches GitHub's 11 files).

The mathematics checks out. The conic hull A·ν - b·y ∈ K is exact by homogeneity of the cone; GEHR (ν'Qν + (a'ν)y + d·y² ≤ 0) and CEHR (ν'Qν ≤ t·y, t + a'ν + d·y ≤ 0, t ≥ 0) reduce to the perspective for y ∈ (0, 1] and to ν = t = 0 at y = 0; the :cehr_conic rows (t/2, y, √λ_k v_k'ν) reproduce ν'Qν ≤ t·y with zero eigenvalues dropped; and the four Big-M interior directions produce no violations at M = 100 on a [-10, 10]³ grid. The docs Josh asked for are in at this head (README, docs/src/index.md, Hull docstring).

The one blocking item is on the InfiniteGDPModel side: the CEHR epigraph variable is finite even when the constraint is infinite, which over-restricts the model (counterexample inline). Details are in the inline comments.

Findings

Inline (6):

• Blocking — src/hull.jl: finite CEHR epigraph variable on infinite models cuts off feasible points.
• Nonblocking — src/bigm.jl: conic Big-M skips the finite-M guard.
• Nonblocking — src/bigm.jl: Float64 cone directions promote coefficients on Float32 models.
• Nonblocking — README.md: document that M is used as given for cones and that the curved cones need M above the variable range.
• Nonblocking — src/constraints.jl: Indicator() on a conic disjunct fails with a raw MethodError.
• Nonblocking — test/solve.jl: no solve-level test for RSOC / exp / power cones or the README example.

Review body (1):

Nonblocking — the description and commit history cover half the change.
The title and body describe conic sets, but the PR also adds the exact quadratic hull reformulations (GEHR / CEHR / :cehr_conic) behind a new public Hull(quadratic = …) keyword with five options, a LinearAlgebra dependency, and Hypatia / Pajarito test dependencies; the commit messages are Initial, cehr_conic tag, Fix tests. Please update the body to describe the quadratic option (with the Gusev & Bernal Neira reference the docs already cite) and squash-merge with a message naming both features — that text is what the release notes and git blame will point at.

Coordination with open PRs

git merge-tree reports textual conflicts with #139 and #140 (Project.toml [compat] / [targets]) and with #141 (README.md, src/constraints.jl). #141 also adds _set_support(::Union{Nonnegatives, Nonpositives, Zeros}) = true to gate which vector sets basic steps accept; after both land that whitelist rejects the conic sets added here, so whichever merges second should decide explicitly whether basic steps accept _ConicSets.

Tests run

• PR head, local: Pkg.test() at 6e87fa2 on Julia 1.10.11 / Linux — passed (32 testsets, all Pass).
• PR head, GitHub: Julia 1.10 on ubuntu / macOS / windows, docs build, codecov/patch, codecov/project — all success for 6e87fa2.
• PR head, scratch environment: cone-direction grid check; GDPModel{Float32} reformulation; Indicator / MBM / PSplit / CuttingPlanes on a conic disjunct; the infinite-model counterexample above under BigM, Hull(), :gehr, :exact.
• Merge result: merge-tree against master is clean; conflicts against #139 / #140 / #141 as listed (no tests run on merge results).

Summary

Blocking: 1. Nonblocking: 6. Questions: 0. Total: 7 (6 inline, 1 in this body).

I would not merge this until the blocking issues above are addressed.

[dnguyen227]

@pulsipher This PR requires the fix from infiniteopt/InfiniteOpt.jl#442. I'm happy with this initial version otherwise.

[pulsipher]

I tagged a new release of InfiniteOpt.