Potential regression in the QP solver of HiGHS 1.14.0

[stephane-caron]

In the continuous integration pipeline of qpsolvers, we have noticed a regression from HiGHS 1.14.0 on solving a small dense QP problem.

The test fails e.g. in this run (I reproduced it locally on my computer as well):

======================================================================
FAIL: test_highs_tolerances (tests.test_highs.TestHiGHS)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/runner/work/qpsolvers/qpsolvers/tests/test_highs.py", line 37, in test_highs_tolerances
    self.assertIsNotNone(x)
AssertionError: unexpectedly None

The corresponding test only checks whether the solver returns a (primal vector in its) solution:

def test_highs_tolerances(self):
    problem = get_sd3310_problem()
    x = highs_solve_qp(
        problem.P, problem.q, problem.G, problem.h, problem.A, problem.b,
        time_limit=0.1,
        primal_feasibility_tolerance=1e-1,
        dual_feasibility_tolerance=1e-1,
    )
    self.assertIsNotNone(x)

The corresponding problem is:

def get_sd3310_problem() -> Problem:
    """
    Get a small dense problem with 3 optimization variables, 3 inequality
    constraints, 1 equality constraint and 0 box constraint.
    """
    M = np.array([[1.0, 2.0, 0.0], [-8.0, 3.0, 2.0], [0.0, 1.0, 1.0]])
    P = np.dot(M.T, M)  # this is a positive definite matrix
    q = np.dot(np.array([3.0, 2.0, 3.0]), M).reshape((3,))
    G = np.array([[1.0, 2.0, 1.0], [2.0, 0.0, 1.0], [-1.0, 2.0, -1.0]])
    h = np.array([3.0, 2.0, -2.0]).reshape((3,))
    A = np.array([1.0, 1.0, 1.0])
    b = np.array([1.0])
    return Problem(P, q, G, h, A, b)

The test does not fail when we the version range of the dependency to highspy = ">=1.5.3,<1.14.0", which suggests highspy 1.14.0 may have a regression there.

Hoping this report helps!

[jajhall]

Thanks, I'll try to reproduce this tomorrow

[jajhall]

I can't reproduce this in C++, getting the same objective value with v1.3.0 and 1.14.0, and same size (3) of x

I can't reproduce this in Python because I don't have qpsolvers installed - and I'm wary of doing so since I'm away from support if I mess up Python on my laptop

Do you get the same issue if you save your model as an MPS (or LP) file, and run from that?

What optimal objective value are you expecting?

[odow]

I also can't reproduce this in JuMP:

julia> using JuMP, HiGHS

julia> begin
           M = [1 2 0; -8 3 2; 0 1 1]
           P = M' * M
           q = M' * [3, 2, 3]
           G = [1 2 1; 2 0 1; -1 2 -1]
           h = [3, 2, -2]
           A = [1 1 1]
           b = [1]
           model = Model(HiGHS.Optimizer)
           @variable(model, x[1:3])
           @objective(model, Min, 0.5 * x' * P * x + q' * x)
           @constraint(model, G * x .<= h)
           @constraint(model, A * x .== b)
           optimize!(model)
           solution_summary(model)
       end
Running HiGHS 1.14.0 (git hash: 7df0786de3): Copyright (c) 2026 under Apache 2.0 license terms
Using BLAS: blastrampoline 
QP has 4 rows; 3 cols; 11 matrix nonzeros; 6 Hessian nonzeros
Coefficient ranges:
  Matrix  [1e+00, 2e+00]
  Cost    [7e+00, 2e+01]
  Hessian [5e+00, 6e+01]
  Bound   [0e+00, 0e+00]
  RHS     [1e+00, 3e+00]
  Iteration        Objective     NullspaceDim
          0        3.5555556                0      0.00s
          6       -2.3076922                1      0.00s
Model status        : Optimal
QP ASM    iterations: 6
Objective value     : -2.3076923077e+00
P-D objective error :  4.4362488438e-08
HiGHS run time      :          0.00
solution_summary(; result = 1, verbose = false)
├ solver_name          : HiGHS
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 1
│ ├ raw_status         : kHighsModelStatusOptimal
│ └ objective_bound    : 7.69234e-02
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : FEASIBLE_POINT
│ └ objective_value      : -2.30769e+00
└ Work counters
  ├ solve_time (sec)   : 2.34042e-04
  ├ simplex_iterations : 0
  ├ barrier_iterations : 0
  └ node_count         : -1

[stephane-caron]

Thank you for taking a look at this. While reproducing the issue I realized there is a difference in the Hessian warning between the two versions of HiGHS: 1.13.1 silently ignores the off-triangle entries of the Hessian matrix of the example, while 1.14.0 symmetrises them ("added to lower triangle, requiring 3 non-trivial summations"). I see now that this was documented in the release notes of v1.14.0.

This change in behavior was the culprit. By changing:

-model.hessian_.dim_ = P_csc.shape[0]
-model.hessian_.start_ = P_csc.indptr
-model.hessian_.index_ = P_csc.indices
-model.hessian_.value_ = P_csc.data
+P_lower = spa.tril(P_csc, format="csc")
+model.hessian_.dim_ = P_lower.shape[0]
+model.hessian_.start_ = P_lower.indptr
+model.hessian_.index_ = P_lower.indices
+model.hessian_.value_ = P_lower.data

The solver returns the expected primal solution as before.

I'm writing down more details below. Feel free to share your thoughts if you have feedback on this fix, and close the issue of course. I will merge a corresponding update in qpsolvers (not updating its highspy compatibility versions since this change should be backward-compatible).

Reproduction details

I reproduced the issue locally as follows, with only highspy, numpy, and scipy as dependencies:

#!/usr/bin/env python3

import numpy as np
import scipy.sparse as spa
import highspy

# Problem definition (sd3310: small dense, 3-3-1-0)
M = np.array([[1.0, 2.0, 0.0], [-8.0, 3.0, 2.0], [0.0, 1.0, 1.0]])
P = M.T @ M          # positive definite Hessian
q = (np.array([3.0, 2.0, 3.0]) @ M).reshape((3,))
G = np.array([[1.0, 2.0, 1.0], [2.0, 0.0, 1.0], [-1.0, 2.0, -1.0]])
h = np.array([3.0, 2.0, -2.0])
A = np.array([[1.0, 1.0, 1.0]])
b = np.array([1.0])

# Convert to CSC sparse matrices (as qpsolvers does internally)
P_csc = spa.csc_matrix(P)
G_csc = spa.csc_matrix(G)
A_csc = spa.csc_matrix(A)
n = q.shape[0]

model = highspy.HighsModel()

# Hessian
model.hessian_.dim_ = P_csc.shape[0]
model.hessian_.start_ = P_csc.indptr
model.hessian_.index_ = P_csc.indices
model.hessian_.value_ = P_csc.data

# Columns
lp = model.lp_
lp.num_col_ = n
lp.col_cost_ = q
lp.col_lower_ = np.full((n,), -highspy.kHighsInf)
lp.col_upper_ = np.full((n,), highspy.kHighsInf)

# Rows: stack [G; A] with appropriate lower/upper bounds
row_matrix = spa.vstack([G_csc, A_csc], format="csc")
lp.num_row_ = row_matrix.shape[0]
lp.a_matrix_.format_ = highspy.MatrixFormat.kColwise
lp.a_matrix_.start_ = row_matrix.indptr
lp.a_matrix_.index_ = row_matrix.indices
lp.a_matrix_.value_ = row_matrix.data
lp.a_matrix_.num_row_ = row_matrix.shape[0]
lp.a_matrix_.num_col_ = row_matrix.shape[1]
lp.row_lower_ = np.hstack([np.full((G.shape[0],), -highspy.kHighsInf), b])
lp.row_upper_ = np.hstack([h, b])

solver = highspy.Highs()

# NB: solver options for tolerance and time limit
solver.setOptionValue("dual_feasibility_tolerance", 1e-1)
solver.setOptionValue("log_to_console", True)
solver.setOptionValue("primal_feasibility_tolerance", 1e-1)
solver.setOptionValue("time_limit", 0.1)

solver.passModel(model)
solver.run()
model_status = solver.getModelStatus()
result = solver.getSolution()
print(f"\nHiGHS version : {solver.versionMajor()}.{solver.versionMinor()}.{solver.versionPatch()}")
print(f"Model status  : {model_status}")
print(f"col_value     : {list(result.col_value)}")

# Print solution
found = model_status == highspy.HighsModelStatus.kOptimal
x = np.array(result.col_value) if found else None
print(f"\nx = {x}")
if x is None:
    print("\n[FAIL] Solver returned no solution — regression reproduced.")
else:
    print("\n[PASS] Solver returned a primal solution.")

I checked the regression locally using a small pixi.toml:

[workspace]
name = "highs-reprod"
version = "0.1.0"
channels = ["conda-forge"]
platforms = ["linux-64"]

[dependencies]
python = "3.11.*"

[pypi-dependencies]
numpy = "*"
scipy = "*"

[feature.highs1131.pypi-dependencies]
highspy = "==1.13.1"

[feature.highs1131.tasks]
test-highspy-1-13-1 = "python repro.py"

[feature.highs1140.pypi-dependencies]
highspy = "==1.14.0"

[feature.highs1140.tasks]
test-highspy-1-14-0 = "python repro.py"

[environments]
highs1131 = ["highs1131"]
highs1140 = ["highs1140"]

The output of pixi run -e highs1131 test-highspy-1-13-1:

Running HiGHS 1.13.1 (git hash: 1d267d9): Copyright (c) 2026 under MIT licence terms
WARNING: Ignored 3 entries of Hessian in opposite triangle
QP has 4 rows; 3 cols; 11 matrix nonzeros; 6 Hessian nonzeros
Coefficient ranges:
  Matrix  [1e+00, 2e+00]
  Cost    [7e+00, 2e+01]
  Hessian [5e+00, 6e+01]
  Bound   [0e+00, 0e+00]
  RHS     [1e+00, 3e+00]
  Iteration        Objective     NullspaceDim
          0        3.5555556                0      0.00s
          6       -2.3076922                1      0.00s
Model status        : Optimal
QP ASM    iterations: 6
Objective value     : -2.3076923077e+00
P-D objective error :  4.4362487172e-08
HiGHS run time      :          0.00

HiGHS version : 1.13.1
Model status  : HighsModelStatus.kOptimal
col_value     : [0.3076923111580727, -0.6923076888419275, 1.3846153776838548]

x = [ 0.30769231 -0.69230769  1.38461538]

[PASS] Solver returned a primal solution.

And here is the output of pixi run -e highs1140 test-highspy-1-14-0:

Running HiGHS 1.14.0 (git hash: 7df0786): Copyright (c) 2026 under MIT licence terms
WARNING: Triangular Hessian contains 3 entries in upper triangle: added to lower triangle, requiring 3 non-trivial summations
QP has 4 rows; 3 cols; 11 matrix nonzeros; 6 Hessian nonzeros
Coefficient ranges:
  Matrix  [1e+00, 2e+00]
  Cost    [7e+00, 2e+01]
  Hessian [5e+00, 6e+01]
  Bound   [0e+00, 0e+00]
  RHS     [1e+00, 3e+00]
  Iteration        Objective     NullspaceDim
          0      -0.22222217                0      0.00s
QP solver model status: Undetermined
Model status        : Solve error
HiGHS run time      :          0.00

HiGHS version : 1.14.0
Model status  : HighsModelStatus.kSolveError
col_value     : [0.0, 0.0, 0.0]

x = None

[FAIL] Solver returned no solution — regression reproduced.

There is a difference in the Hessian warning between the two versions: 1.13.1 silently ignores the off-triangle entries, while 1.14.0 symmetrises them ("added to lower triangle, requiring 3 non-trivial summations"). I see now that this was documented in the release notes of v1.14.0.

This change in behavior was the culprit here. If we change:

-model.hessian_.dim_ = P_csc.shape[0]
-model.hessian_.start_ = P_csc.indptr
-model.hessian_.index_ = P_csc.indices
-model.hessian_.value_ = P_csc.data
+P_lower = spa.tril(P_csc, format="csc")
+model.hessian_.dim_ = P_lower.shape[0]
+model.hessian_.start_ = P_lower.indptr
+model.hessian_.index_ = P_lower.indices
+model.hessian_.value_ = P_lower.data

Then the solver returns the expected primal solution as before.