Dual simplex: phase-1 basis history blocks phase-2 progress (self-contained synthetic reproducer)

[liujianing055-afk]

A cold-started, feasible, bounded synthetic LP returns Unknown after a cycling rejection immediately after the dual phase-1 to phase-2 transition. The two revisited bases occurred in phase 1 under auxiliary bounds; phase 2 has made only one accepted pivot when the second candidate is prohibited.

This seems related to the anti-cycling concerns in #1534 and #1607, but the reproducer below requires no warm start or application data.

Versions and results

Reproduced with the official Windows highspy==1.15.1 wheel (HiGHS git 04024d7, CPython 3.13.3, NumPy 2.5.1). Also reproduced through the C API on Linux using source commit 9a1dbd1e2d331f8868a7b13c6ce98700ed31afef from latest, GCC 13.3, CPython 3.12.3, NumPy 2.3.5. The Linux test reproduces when regenerating the synthetic input there as well. An older HiGHS 1.8.0 build exhibits the same failure.

All runs use serial dual simplex, presolve off, and primal/dual feasibility tolerances of 1e-10.

Source build at the commit above	Scaling strategy	Status	Iterations	Minimum x	Independently recomputed max abs(Ax-b)
Unmodified	1	Unknown	51	-2.3378374643	9.3e-15
Unmodified	2	Optimal	83	0	1.8e-14
Candidate phase-history reset below	1	Optimal	135	0	6.1e-15

The successful runs give an objective of about -0.53164328324541. The rejected vector has objective -0.54793298503713 and violates a nonnegativity bound substantially. HiGHS correctly labels that vector Unknown; this report concerns premature termination, not a false Optimal status.

Self-contained Python reproducer

The instance is constructed entirely from the literals below. It reads no files and requires only NumPy and highspy. The five equalities and x >= 0 include sum(x)=1, so the feasible region is bounded. The vector w constructs a strictly positive numerical feasibility witness.

"""Independent synthetic LP; no input files or application data are required.

Install numpy and highspy, then run: python reproduce.py
"""
import os
for name in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[name] = "1"
import json
import numpy as np
import highspy

n = 1280
t = 2 * np.pi * np.arange(n) / n
A = np.vstack((np.ones(n), np.cos(t), np.sin(t), np.cos(2*t), np.sin(2*t)))
a = np.array([-7, -9, -16, 8, -1, 8, 9, 3, -2, 4, -5, -6], dtype=float) / 16
c = sum(a[2*k-2] * np.cos(k*t) + a[2*k-1] * np.sin(k*t) for k in range(1, 7))
c = -(c - c.min()) / (c.max() - c.min())
w = np.exp(3 * np.cos(t - 2*np.pi*0.40154194513298813))
w /= w.sum()
b = A @ w

lp = highspy.HighsLp()
lp.num_col_, lp.num_row_ = n, 5
lp.col_cost_, lp.col_lower_, lp.col_upper_ = c, np.zeros(n), np.full(n, np.inf)
lp.row_lower_, lp.row_upper_ = b, b
lp.a_matrix_.num_col_, lp.a_matrix_.num_row_ = n, 5
lp.a_matrix_.format_ = highspy.MatrixFormat.kColwise
lp.a_matrix_.start_ = np.arange(n+1, dtype=np.int32) * 5
lp.a_matrix_.index_ = np.tile(np.arange(5, dtype=np.int32), n)
lp.a_matrix_.value_ = A.T.ravel()

print("NumPy", np.__version__, "HiGHS", highspy.Highs().version(), highspy.Highs().githash())
for strategy in (1, 2):
    h = highspy.Highs()
    for name, value in {
        "solver": "simplex", "simplex_strategy": 1,
        "simplex_scale_strategy": strategy, "presolve": "off",
        "threads": 1, "parallel": "off",
        "primal_feasibility_tolerance": 1e-10,
        "dual_feasibility_tolerance": 1e-10,
        "simplex_iteration_limit": 1000, "output_flag": False,
    }.items():
        assert h.setOptionValue(name, value) == highspy.HighsStatus.kOk
    assert h.passModel(lp) != highspy.HighsStatus.kError
    h.run()
    x = np.array(h.getSolution().col_value)
    y = np.array(h.getSolution().row_dual)
    print(json.dumps({
        "scale_strategy": strategy,
        "status": h.modelStatusToString(h.getModelStatus()),
        "iterations": h.getInfo().simplex_iteration_count,
        "min_x": float(x.min()),
        "max_abs_Ax_minus_b": float(np.max(np.abs(A @ x - b))),
        "objective": float(c @ x),
        "min_dual_slack": float(np.min(c - A.T @ y)),
        "primal_minus_dual": float(c @ x - b @ y),
    }))

Trace of the cross-phase revisit

In a trace-only build, full basis index sets (not just hashes) match as follows:

Phase-2 candidate iteration	Candidate basis, sorted	Same set recorded after a phase-1 pivot at iteration_count_
50	[107, 108, 343, 587, 977]	48
51	[107, 343, 587, 977, 978]	47

Phase 2 starts at iteration 50. The first revisit is accepted. The next is rejected by HEkk::isBadBasisChange because the two revisit detections have consecutive iteration numbers. Thus the history responsible for the rejection predates restoration of the original bounds. The trace-only build returns exactly the same primal and dual arrays as the unmodified build.

Relevant original log excerpt:

dual-phase-1-optimal
dual-phase-2-start
DuPh2         50    -5.4795712499e-01 Pr: 2(0.201424) No reason
 basis change (108 out; 978 in) is bad
DuPh2         51    -5.4793251435e-01 Pr: 1(0.146115) Possibly optimal
 basis change (108 out; 978 in) is bad

For an additional feasibility check, exact rational elimination on the five positive columns returned by the successful control constructs a strictly positive solution of all five binary64 equalities. Rational dual-slack checks also bracket its objective with a gap below 3e-14. These checks were made both with the generated matrix and with entries at or below the default small-matrix threshold removed.

Candidate change for discussion

Resetting the visited-basis history only when phase 1 actually transitions to phase 2 restores this example without changing the LP, scaling strategy, or numerical tolerances:

--- a/highs/simplex/HEkkDual.cpp
+++ b/highs/simplex/HEkkDual.cpp
@@ -797,6 +797,11 @@
     ekk_instance_.initialiseBound(SimplexAlgorithm::kDual, kSolvePhase2);
     ekk_instance_.initialiseNonbasicValueAndMove();
     if (solve_phase == kSolvePhase2) {
+      // Phase 1 used auxiliary bounds: its visited bases do not establish a
+      // cycle under the restored phase 2 bounds.
+      ekk_instance_.visited_basis_.clear();
+      ekk_instance_.visited_basis_.insert(ekk_instance_.basis_.hash);
+      ekk_instance_.previous_iteration_cycling_detected = -kHighsIInf;
       // Moving to phase 2 so possibly reinstate cost perturbation
       if (ekk_instance_.dual_simplex_phase1_cleanup_level_ <
           ekk_instance_.options_->max_dual_simplex_phase1_cleanup_level) {

This keeps the phase-2 cycling guard active. With this candidate change, the existing unit_tests_all suite passed (365 test cases, 1,259,732 assertions), and capi_unit_tests passed. This is a tested suggestion, not a claim that this small check establishes a general anti-cycling guarantee.

Would a phase-local history reset be appropriate here, or would you prefer a retry mechanism along the lines of #1534?

AI assistance was used to construct and test this independent reproducer and investigate the code path.

[jajhall]

Thanks for this. I've reproduced the behaviour in a C++ unit test, and will study it further