Troubleshooting SciPy 1.15.2 MILP with HiGHS: O(N^2) postsolve packaging causes stalls; upgrade to 1.15.3/1.16.0

Seeing a long stall after HiGHS reports in SciPy 1.15.2 MILP? It is an O(N^2) wrapper slowdown, not time_limit. Upgrade to 1.15.3/1.16.0 to match timing.

When a MILP model in SciPy prints the HiGHS solving report and then the Python script stalls for minutes before returning, it looks like the time limit was ignored or the solver hung. In reality, this behavior can come from a specific slowdown in SciPy 1.15.2 that happens after the solver finishes: the result packaging step inside SciPy can dominate the wall time and create a long gap between the printed report and your next print statement.

Reproducible example that shows the delay

The snippet below mirrors a typical MILP setup using SciPy’s milp, but with different identifier names for clarity. It creates a random dataset, assembles the variables and linear constraints, solves with a time limit, and prints the total duration measured in Python. On SciPy 1.15.2 this can be much longer than the “Timing (total)” reported by HiGHS.

import pandas as pd
import numpy as np
from scipy.optimize import Bounds, milp, LinearConstraint
import time
import scipy
print(scipy.__version__)
cap_hi = 2
cap_lo = 1
fr_hi = 1
fr_lo = 1
N_total = 300
np.random.seed(42)
data_tbl = pd.DataFrame({'Col1': np.random.rand(N_total), 'Col2': np.random.rand(N_total)})
def run_and_collect(batch_df, max_wall):
    n_rows = len(batch_df)
    if n_rows < 10:
        return [], batch_df.copy(), False
    group_cnt = n_rows // 10
    cap_arr = batch_df['Col1'].astype(float).to_numpy()
    fr_arr  = batch_df['Col2'].astype(float).to_numpy()
    obj_vec = np.concatenate([np.zeros(n_rows*group_cnt), -np.ones(group_cnt), np.zeros(n_rows)])
    integrality_mask = np.ones(n_rows*group_cnt + group_cnt + n_rows, dtype=int)
    dom_limits = Bounds(0, 1)
    L1 = np.zeros((n_rows, n_rows*group_cnt + group_cnt + n_rows)); r1_up = np.zeros(n_rows)
    L2 = np.zeros((n_rows, n_rows*group_cnt + group_cnt + n_rows)); r2_up = np.zeros(n_rows)
    for i in range(n_rows):
        L1[i, i + np.arange(group_cnt)*n_rows] = 1
        L1[i, n_rows*group_cnt + group_cnt + i] = -1
        L2[i, i + np.arange(group_cnt)*n_rows] = -1
        L2[i, n_rows*group_cnt + group_cnt + i] =  1
    C2a = np.zeros((group_cnt, n_rows*group_cnt + group_cnt + n_rows)); rhs2a = np.zeros(group_cnt)
    C2b = np.zeros((group_cnt, n_rows*group_cnt + group_cnt + n_rows)); rhs2b = np.zeros(group_cnt)
    for j in range(group_cnt):
        C2a[j, j*n_rows:(j+1)*n_rows] = 1
        C2a[j, n_rows*group_cnt + j]  = -10
        C2b[j, j*n_rows:(j+1)*n_rows] = -1
        C2b[j, n_rows*group_cnt + j]  = 10
    C3 = np.zeros((group_cnt, n_rows*group_cnt + group_cnt + n_rows)); rhs3 = np.zeros(group_cnt)
    C4 = np.zeros((group_cnt, n_rows*group_cnt + group_cnt + n_rows)); rhs4 = np.zeros(group_cnt)
    C5 = np.zeros((group_cnt, n_rows*group_cnt + group_cnt + n_rows)); rhs5 = np.zeros(group_cnt)
    C6 = np.zeros((group_cnt, n_rows*group_cnt + group_cnt + n_rows)); rhs6 = np.zeros(group_cnt)
    for j in range(group_cnt):
        C3[j, j*n_rows:(j+1)*n_rows]   = cap_arr
        C3[j, n_rows*group_cnt + j]    = -cap_hi*10
        C4[j, j*n_rows:(j+1)*n_rows]   = -cap_arr
        C4[j, n_rows*group_cnt + j]    =  cap_lo*10
        C5[j, j*n_rows:(j+1)*n_rows]   = fr_arr
        C5[j, n_rows*group_cnt + j]    = -fr_hi*10
        C6[j, j*n_rows:(j+1)*n_rows]   = -fr_arr
        C6[j, n_rows*group_cnt + j]    =  fr_lo*10
    Csym = np.zeros((group_cnt-1, n_rows*group_cnt + group_cnt + n_rows)); rhs_sym = np.zeros(group_cnt-1)
    for j in range(group_cnt-1):
        Csym[j, n_rows*group_cnt + j]   = 1
        Csym[j, n_rows*group_cnt + j+1] = -1
    C7 = np.zeros((1, n_rows*group_cnt + group_cnt + n_rows))
    C7[0, n_rows*group_cnt + group_cnt:] = 1
    C7[0, n_rows*group_cnt:n_rows*group_cnt + group_cnt] = -10
    cons = [
        LinearConstraint(C2a, -np.inf, rhs2a),
        LinearConstraint(C2b, -np.inf, rhs2b),
        LinearConstraint(C3,  -np.inf, rhs3),
        LinearConstraint(C4,  -np.inf, rhs4),
        LinearConstraint(C5,  -np.inf, rhs5),
        LinearConstraint(C6,  -np.inf, rhs6),
        LinearConstraint(C7,  -np.inf, 0),
        LinearConstraint(L1,  -np.inf, r1_up),
        LinearConstraint(L2,  -np.inf, r2_up),
        LinearConstraint(Csym, rhs_sym, np.inf)
    ]
    t0 = time.time()
    ans = milp(
        c=obj_vec,
        integrality=integrality_mask,
        bounds=dom_limits,
        constraints=cons,
        options={'disp': True, 'time_limit': max_wall}
    )
    t1 = time.time()
    print(f'Duration: {t1 - t0:.2f}')
run_and_collect(data_tbl, 10)

What actually causes the stall

This behavior aligns with a known issue in SciPy 1.15.2 that affects HiGHS-based optimization. Inside SciPy’s _highs_wrapper(), reading col_status is implemented with an O(N^2) loop. That work happens after the solver stops, while SciPy is building the result object. Because it is postsolve bookkeeping, it is not governed by the time_limit option, which is why you see the solver’s report printed promptly, but Python returns much later. In practice, even small-to-medium problems can be impacted if the number of variables is large enough for the O(N^2) step to dominate.

Fix: upgrade SciPy

The practical resolution is to upgrade SciPy to a version that eliminates this slowdown. Moving from 1.15.2 to 1.15.3 or 1.16.0 avoids the O(N^2) col_status handling and brings the return time in line with the solver’s own timing. No changes to your model or code are required.

If you want to verify what is installed at runtime, you can print the version:

import scipy
print(scipy.__version__)

On the same benchmark setup, the observed durations shrink dramatically after the upgrade. For instance, with the example above the HiGHS “Timing (total)” is about 1.6 seconds, while the end-to-end Python Duration drops from roughly 48 seconds on 1.15.2 to about 1.7 seconds on 1.15.3.

Why this matters

MILP workflows often enforce wall-time contracts via time_limit and instrument code with timers around solver calls. When postsolve packaging becomes the bottleneck, those contracts silently break: jobs appear to exceed budgets despite the solver finishing on time. The mismatch leads to confusing telemetry, brittle pipelines, and wasted compute. Eliminating the overhead restores the expected relationship between the solver’s report and the application’s wall clock.

Takeaways

If you encounter large, variable delays after the HiGHS solving report in SciPy 1.15.2, the slowdown likely comes from result preparation in the wrapper rather than the MILP itself. Upgrading to SciPy 1.15.3 or 1.16.0 resolves the issue and requires no changes to your modeling code. Keep an eye on the printed SciPy version and compare your own wall-clock measurements with the solver’s “Timing (total)” to confirm the improvement.