Fix non-rigorous crossed bound results by bernalde · Pull Request #4013 · Pyomo/pyomo

bernalde · 2026-08-15T17:53:12Z

Fixes #3947.

Summary/Motivation:

When GDPopt LOA is applied to a nonconvex problem, its linearizations are not valid relaxations, so the "dual bound" it accumulates is not a rigorous bound. If those bounds cross, GDPopt currently snaps the dual bound to the primal bound and reports optimal with LB == UB. That presents a non-rigorous result as a certified global optimum, which is the failure reported in #3947.

The fix reports what is actually known in that case: the incumbent is feasible, and no certified dual bound is available.

This is a replacement for #3948, which GitHub would not allow me to reopen after it was closed. The prior review discussion and history remain available there.

Changes proposed in this PR:

Add a _crossed_bounds_are_certified class flag to _GDPoptAlgorithm (default True), set to False on GDP_LOA_Solver. GLOA keeps the existing certified behavior.
When bounds cross on a solver whose bounds are not certified and the model is detected as possibly non-rigorous, terminate with feasible instead of optimal, and report the uncertified side of the bound as infinite rather than snapping it to the primal bound.
Add _problem_may_have_nonrigorous_dual_bound() in pyomo/contrib/gdpopt/loa.py, which inspects the active objective and constraints: anything non-polynomial or of degree > 2 is treated as possibly non-rigorous, quadratic equalities are treated as possibly non-rigorous, and quadratic inequalities are classified by curvature against their bound sense.
Classify quadratic curvature from the full symmetric Q matrix rather than diagonal terms only, so convex forms with cross terms such as x**2 + x*y + y**2 <= 4 are still certified.
Use NumPy eigvalsh when available, with a pure-Python Jacobi eigenvalue fallback so the classification also works in Pyomo's NumPy-free test environments.
Apply the same crossed-bound result semantics to MindtPy GOA.
Add regressions in pyomo/contrib/gdpopt/tests/test_gdpopt.py and pyomo/contrib/mindtpy/tests/test_mindtpy_no_discrete.py covering crossed-bound reporting, non-polynomial detection, quadratic convexity classification (including the no-NumPy path), and preservation of certified GLOA behavior.
Validation performed locally on this branch, rebased onto current main:
- python -m pytest -q pyomo/contrib/gdpopt/tests/test_gdpopt.py pyomo/contrib/mindtpy/tests/test_mindtpy_no_discrete.py
  - Result: 96 passed, 19 skipped, 3 deselected in 43.01s
- python -m black --check on the six changed files
  - Result: passed, 6 files would be left unchanged
- typos --config ./.github/workflows/typos.toml on the six changed files
  - Result: passed, no findings

AI-Use Disclosure

AI tools were NOT used during the preparation of this PR

or

Notes for reviewers (optional): This replacement PR carries over the implementation and tests from #3948 unchanged, then refreshes the branch against current main. The replacement PR description and branch-refresh workflow were prepared with AI assistance and reviewed before posting.

This is the largest of the four GDPopt fixes I have open and it changes reported termination conditions, so it deserves the closest look. Specific points where I would like reviewer judgement:

Behavior change in reported status. Nonconvex LOA runs that previously returned optimal with LB == UB will now return feasible with one bound infinite. That is the intent of GDPopt LOA reports non-rigorous nonconvex bounds as optimal LB=UB #3947, but it is user-visible and may affect downstream code that keys on optimal.
Asymmetry between GDPopt and MindtPy. GDPopt gates the new behavior on both not _crossed_bounds_are_certified and _problem_may_have_nonrigorous_dual_bound(model). The MindtPy GOA path in bounds_converged() gates only on not _crossed_bounds_are_certified, so any crossed bound there becomes feasible. I would like a decision on whether MindtPy should carry the same model-structure gate.
GLOA vs GOA. GDPopt GLOA keeps _crossed_bounds_are_certified = True while MindtPy GOA is set to False, even though they are the analogous global algorithms. That was deliberate — GLOA's certified behavior is covered by a regression test — but the inconsistency is worth an explicit maintainer decision.
Conservatism of the classifier. _problem_may_have_nonrigorous_dual_bound() is intentionally conservative: quadratic equality constraints and anything of degree > 2 or non-polynomial are all treated as possibly non-rigorous. Convex problems expressed in ways the classifier cannot certify will lose their optimal status.
The pure-Python Jacobi fallback. It exists so the classification works without NumPy. Sampled eigenvalues were checked against NumPy and a no-NumPy regression test is included, but the rotation cap (50 * order**2) returning None on non-convergence — which is then treated as "not certified" — is a heuristic worth confirming.

Legal Acknowledgement

By contributing to this software project, I have read the contribution guide and agree to the following terms and conditions for my contribution:

I agree my contributions are submitted under the BSD license.
I represent I am authorized to make the contributions and grant the license. If my employer has rights to intellectual property that includes these contributions, I represent that I have received permission to make contributions and grant the required license on behalf of that employer.

Pyomo#3947) Extract the convexity/curvature classification out of GDP_LOA_Solver into a shared pyomo/contrib/gdpopt/nonrigorous_bounds.py module, and drive it from a single _crossed_bounds_are_certified class flag on both the GDPopt and MindtPy algorithm base classes. An algorithm is certified when its relaxation is valid for any model it accepts. The McCormick-based global algorithms (GDPopt GLOA, MindtPy GOA) are certified. The algorithms that linearize at trial points (GDPopt LOA, MindtPy OA and ECP) are only valid on convex models, so they are not certified and the model is checked for convexity before crossed bounds are trusted. This corrects the MindtPy classification, which previously marked GOA (the McCormick-based algorithm) as uncertified while leaving OA and ECP certified, and adds MindtPy the model-structure gate that GDPopt already had.

bernalde · 2026-08-15T18:15:59Z

Pushed f4c13fd52 — this resolves review points 2 and 3 in the PR description above, which flagged an inconsistency between the GDPopt and MindtPy halves of this PR. Rather than leaving them for a reviewer to arbitrate, I traced the cause and generalized the mechanism.

The inconsistency

The two toolkits have matching algorithm pairs, and the distinguishing property is how each builds its relaxation:

Algorithm	Relaxation	Rigorous dual bound?	Previously in this PR
GDPopt LOA	linearize at trial points	only if convex	`certified = False` — correct
GDPopt GLOA	McCormick envelopes	yes	`certified = True` — correct
MindtPy OA	linearize at trial points (`add_oa_cuts`)	only if convex	`certified = True` — wrong
MindtPy ECP	linearize at trial points (`add_ecp_cuts`)	only if convex	`certified = True` — wrong
MindtPy GOA	McCormick envelopes (`add_affine_cuts`, `use_mcpp = True`)	yes	`certified = False` — wrong

The flag was on the wrong MindtPy solver. GOA is the McCormick-based algorithm — the direct analogue of GLOA, which this PR deliberately leaves certified — while OA and ECP are the analogues of LOA and were left untouched.

What changed

Extracted the curvature/convexity classification out of GDP_LOA_Solver into a new shared module, pyomo/contrib/gdpopt/nonrigorous_bounds.py.
Both algorithm base classes now implement _problem_may_have_nonrigorous_dual_bound() once: return False when the algorithm is certified, otherwise run the model check. Individual solvers only declare _crossed_bounds_are_certified, so no algorithm carries its own copy of this logic.
Moved the MindtPy flag: GOA is certified, OA and ECP are not.
Gave MindtPy the model-structure gate GDPopt already had, so a crossed bound only downgrades to feasible when the model itself could not be certified convex. Previously any crossed bound in MindtPy downgraded unconditionally.

A side benefit of routing through the flag: certified algorithms now skip the convexity scan entirely instead of relying on a base method that returned False.

Effect on reported results

For MindtPy this is a behavior change in both directions, and both directions look like the intended fix:

GOA no longer downgrades to feasible on crossed bounds. Its bounds are rigorous, so a crossing is numerical and optimal is correct.
OA and ECP now downgrade to feasible on crossed bounds, but only for models that fail the convexity check. Convex models are unaffected.

Tests

TestMindtPyGOAResults became TestMindtPyCrossedBoundResults and now covers OA and ECP downgrading, OA on a convex model staying optimal, GOA preserving certified behavior, the flag matrix itself, and certified algorithms skipping the model scan. The GDPopt no-NumPy test was repointed at the new module.

pytest -q pyomo/contrib/mindtpy/tests/test_mindtpy_no_discrete.py pyomo/contrib/gdpopt/tests/ → 115 passed, 35 skipped, 5 deselected
pytest -q pyomo/contrib/mindtpy/tests/ → 58 passed, 32 skipped
black --check on all nine changed files → 9 files would be left unchanged
typos --config ./.github/workflows/typos.toml pyomo/contrib/gdpopt/ pyomo/contrib/mindtpy/ → clean

Still worth a maintainer decision

Module location. nonrigorous_bounds.py sits under pyomo/contrib/gdpopt/ because MindtPy already imports from pyomo.contrib.gdpopt in eight modules, so this follows existing precedent. If a neutral location is preferred, it moves cleanly.
MindtPy_FP_Solver is left at the default (certified). The feasibility pump is a primal heuristic rather than a bounding algorithm, and it does not reference bounds_converged or the dual bound directly, so I did not classify it. Worth a second opinion.
The original crossing evidence. The GOA test in the previous revision used primal 1011.6577899409375 against dual 1000.0 — a crossing of roughly 11.7 against a 1e-6 tolerance, which is far too large to be numerical noise. If that came from a real GOA run, then GOA's dual bound was not rigorous in that instance, and the underlying cause is a separate bug in the McCormick cut generation rather than something this flag should paper over. I could not reproduce the run to confirm it, so I have preserved the number in the tests but classified GOA on the basis of what the algorithm is designed to do. If anyone recognizes that instance, it is worth chasing separately.

bernalde · 2026-08-15T18:24:31Z

Following up on open point 3 from my previous comment, where I flagged that the crossed-bound figures used in the tests (primal 1011.6577899409375 against dual 1000.0, a gap of roughly 11.7 against a 1e-6 tolerance) looked too large to be numerical noise, and asked whether they came from a real MindtPy GOA run that would indicate a non-rigorous GOA dual bound.

They did not. The point is resolved, and no code change is needed.

Tracing the provenance:

The value originates in the reproducer in GDPopt LOA reports non-rigorous nonconvex bounds as optimal LB=UB #3947, which is a GDPopt LOA reproducer. The issue text labels it explicitly — "This mimics a nonconvex LOA state" — and constructs the state by hand:
```
solver.LB = 1011.6577899409375   # feasible incumbent from a local NLP subproblem
solver.UB = 1000.0               # non-rigorous OA/discrete bound that has crossed it
```
The precise LB is an incumbent objective from GDPLib's HDA model; the UB is an arbitrary round number chosen to force the crossing. It is a synthetic fixture for exercising the result-finalization path, not a measurement.
The same pair is used by _setup_crossed_bound_state() in the GDPopt tests here. The MindtPy test added in the previous revision reused those same two numbers in a GOA-shaped unit test. So it was the LOA fixture copied across, not an independent GOA observation.
The number appears nowhere else in Pyomo — only in tests introduced by this PR.
Searching Pyomo issues for reports of MindtPy GOA, or MindtPy generally, producing crossed bounds returns nothing.
Consistent with this, the GDPLib HDA tracking issue (Refactor hda to fix benchmark issues SECQUOIA/gdplib#67, referenced from GDPopt LOA reports non-rigorous nonconvex bounds as optimal LB=UB #3947) lists LOA among the strategies causing benchmark problems, while GLOA is not checked off.

So there is no evidence that a McCormick-based algorithm ever produced a genuinely crossed bound, and no separate cut-generation bug implied. The reclassification in f4c13fd52 — GOA certified, OA and ECP not — rests on what each algorithm's relaxation actually guarantees, and nothing in the record contradicts it.

The figures remain in the tests as a synthetic fixture, which is all they ever were. Open points 1 (module location) and 2 (MindtPy_FP_Solver left unclassified) from my previous comment still stand for maintainer input.

emma58

Thank you: This is a nice addition to keep users using from shooting themselves in the foot when we can guarantee a problem is convex. I had quite a few comments: I find all the names confusing--I'd vote for names that actually say something about certified or not certified convexity rather than about the crossed bounds. I also would strongly prefer that the implementation not depend on generate_standard_repn and instead use the QuadraticRepnVisitor. Also I'm curious why it is imperative that this work without numpy?

emma58 · 2026-09-01T23:13:20Z

+from pyomo.common.dependencies import numpy as np, numpy_available
+from pyomo.core import Block, Constraint, Objective, minimize, value
+from pyomo.gdp import Disjunct
+from pyomo.repn import generate_standard_repn


Please replace use of generate_standard_repn with the LinearRepnVisitor or the QuadraticRepnVisitor. It has known issues and we are hoping to eventually deprecate it. In this case, it looks like you want the QuadraticRepnVisitor.

Replaced generate_standard_repn with QuadraticRepnVisitor in convexity.py.

emma58 · 2026-09-01T23:15:58Z

+    Uses NumPy when it is available. Otherwise falls back to a cyclic Jacobi
+    iteration so that curvature classification still works in environments
+    without NumPy. Returns None if the Jacobi iteration does not converge.
+    """


Is it crucial that this work in an environment without numpy? Implementing a fallback seems like a lot of technical debt, and I'm not sure for what gain. It would not be unprecedented in Pyomo to have a piece of functionality that does truly depend on numpy--that would be fine for GDPopt.

Removed the custom Jacobi fallback. Quadratic curvature checks now use numpy.linalg.eigvalsh

emma58 · 2026-09-01T23:19:36Z

+    if q_matrix is None:
+        return None
+
+    eigenvalue_tolerance = 1e-10


This should be configurable in the options for the global algorithms

Added eigenvalue_tolerance to the GDPopt and MindtPy solver options, with a default of 1e-10.

emma58 · 2026-09-01T23:22:07Z

+    return None
+
+
+def model_may_have_nonrigorous_dual_bound(model):


I'd propose renaming this to model_is_not_certified_convex or something along those lines since that's what you're actually checking.

Renamed the check to model_is_not_certified_convex.

emma58 · 2026-09-01T23:24:08Z

+    approximation dual bound computed for it must not be treated as rigorous.
+    """
+    for obj in model.component_data_objects(Objective, active=True, descend_into=True):
+        degree = obj.expr.polynomial_degree()


polynomial_degree has some known issues also. It may be safer to walk this with the QuadraticRepnVisitor and draw your own conclusion from the result.

Removed the polynomial_degree() checks. QuadraticRepnVisitor now determines whether each expression is quadratic or contains remaining nonlinear terms.

emma58 · 2026-09-01T23:35:53Z

+        self._nonrigorous_crossed_bounds = False
+        self._nonrigorous_dual_bound_possible = False


I'd advocate renaming these to something about whether or not convexity is certified, but if you prefer this, that's okay.

Renamed the state to _model_is_not_certified_convex and _bounds_crossed_without_certified_convexity.

emma58 · 2026-09-01T23:37:44Z

 class _GDPoptAlgorithm:
    CONFIG = ConfigBlock("GDPopt")
    _add_common_configs(CONFIG)
+    _crossed_bounds_are_certified = True


This name is confusing to me. Certified what? Seems like this should be _treat_problem_as_convex or something like that.

Replaced the ambiguous certification flag with _requires_model_convexity. LOA, OA, and ECP set it to True; GLOA and GOA retain the default False.

emma58 · 2026-09-02T13:18:06Z

+        try:
+            self.assertFalse(
+                GDP_LOA_Solver()._problem_may_have_nonrigorous_dual_bound(convex)
+            )
+            self.assertTrue(
+                GDP_LOA_Solver()._problem_may_have_nonrigorous_dual_bound(nonconvex)
+            )
+        finally:
+            nonrigorous_bounds_module.numpy_available = original_numpy_available
+


This is unnecessary: Just import numpy_available here and guard the whole test to only run if it is available.

Removed the Jacobi-specific test and guarded the NumPy-dependent curvature tests with numpy_available.

codecov · 2026-09-02T14:09:31Z

Codecov Report

❌ Patch coverage is 97.19626% with 3 lines in your changes missing coverage. Please review.
✅ Project coverage is 89.99%. Comparing base (ce205a5) to head (4f51f8e).
⚠️ Report is 61 commits behind head on main.

Files with missing lines	Patch %	Lines
pyomo/contrib/gdpopt/convexity.py	95.23%	3 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #4013      +/-   ##
==========================================
+ Coverage   89.93%   89.99%   +0.05%     
==========================================
  Files         917      918       +1     
  Lines      109135   109866     +731     
==========================================
+ Hits        98148    98869     +721     
- Misses      10987    10997      +10

Flag	Coverage Δ
builders	`29.07% <16.82%> (-0.02%)`	⬇️
default	`86.14% <97.19%> (?)`
expensive	`35.12% <61.68%> (?)`
linux	`87.69% <97.19%> (-1.95%)`	⬇️
linux_other	`87.69% <97.19%> (+0.02%)`	⬆️
oldsolvers	`27.54% <16.82%> (-0.01%)`	⬇️
osx	`83.13% <97.19%> (+0.03%)`	⬆️
win	`85.45% <97.19%> (+0.03%)`	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

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

bernalde · 2026-09-02T23:00:53Z

Addressed in 4f51f8e: the naming now describes convexity certification, expression inspection uses QuadraticRepnVisitor, the custom NumPy fallback was removed, the eigenvalue tolerance is configurable, and an end-to-end nonconvex LOA regression test was added.

emma58

One last comment on import statements, but otherwise this looks good, thank you!

emma58 · 2026-09-05T14:54:39Z

+
+class TestMindtPyCrossedBoundResults(unittest.TestCase):
+    def test_oa_crossed_bounds_are_not_reported_as_global_optimal(self):
+        from pyomo.contrib.mindtpy.outer_approximation import MindtPy_OA_Solver


Can you please move this import and the others below to the top of the file?

bernalde added 4 commits August 15, 2026 13:48

Fix nonrigorous crossed bound results (Pyomo#3947)

2bb78c7

Certify LOA convex quadratic cross terms (Pyomo#3947)

2a1737a

Support LOA quadratic certification without NumPy (Pyomo#3947)

b37c09f

This was referenced Aug 15, 2026

Fix non-rigorous crossed bound results #3948

Closed

Epic: 2020 GDP/MINLP roadmap revisited — status audit and open ideas (meeting_20200527) bernalde/pyomo#2

Open

blnicho requested a review from emma58 September 1, 2026 19:26

Test Jacobi fallback without dependency patching

b0cc3a4

emma58 reviewed Sep 2, 2026

View reviewed changes

emma58 added the AT: PRE-TEST INSPECTED label Sep 2, 2026

pyomo-autotest removed the AT: PRE-TEST INSPECTED label Sep 2, 2026

Address convexity review feedback

4f51f8e

emma58 reviewed Sep 5, 2026

View reviewed changes

emma58 added the AT: PRE-TEST INSPECTED label Sep 5, 2026

pyomo-autotest removed the AT: PRE-TEST INSPECTED label Sep 5, 2026

		return None


		def model_may_have_nonrigorous_dual_bound(model):

		self._nonrigorous_crossed_bounds = False
		self._nonrigorous_dual_bound_possible = False

Uh oh!

Conversation

bernalde commented Aug 15, 2026

Fixes #3947.

Summary/Motivation:

Changes proposed in this PR:

AI-Use Disclosure

Legal Acknowledgement

Uh oh!

bernalde commented Aug 15, 2026

The inconsistency

What changed

Effect on reported results

Tests

Still worth a maintainer decision

Uh oh!

bernalde commented Aug 15, 2026

Uh oh!

emma58 left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

codecov Bot commented Sep 2, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

bernalde commented Sep 2, 2026

Uh oh!

emma58 left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

codecov Bot commented Sep 2, 2026 •

edited

Loading