Fix Hessian matrix in HiGHS

We should sort the elements in the same column by its row number

Author: metab0t

include/pyoptinterface/solver_common.hpp

@@ -40,7 +40,8 @@ class OnesideLinearConstraintMixin
 	{
 		if (function.degree() >= 2)
 		{
-			throw std::runtime_error("add_linear_constraint expects linear expression but receives a quadratic expression.");
+			throw std::runtime_error("add_linear_constraint expects linear expression but receives "
+			                         "a quadratic expression.");
 		}
 		ScalarAffineFunction f(function);
 		return get_base()->add_linear_constraint(f, sense, rhs, name);
@@ -617,7 +618,16 @@ struct CSCMatrix
 		// Sorting based on column indices
 		std::vector<IDXT> idx(numnz);
 		std::iota(idx.begin(), idx.end(), 0);
-		std::sort(idx.begin(), idx.end(), [&](int i, int j) { return cols[i] < cols[j]; });
+		std::sort(idx.begin(), idx.end(), [&](int i, int j) {
+			if (cols[i] == cols[j])
+			{
+				return rows[i] < rows[j];
+			}
+			else
+			{
+				return cols[i] < cols[j];
+			}
+		});
 
 		// Creating CSC arrays
 		values_CSC.reserve(numnz);

tests/test_qp.py

@@ -1,6 +1,8 @@
 import pyoptinterface as poi
 from pytest import approx
 
+import numpy as np
+
 
 def test_simple_qp(model_interface):
     model = model_interface
@@ -21,3 +23,46 @@ def test_simple_qp(model_interface):
 
     obj_val = model.get_model_attribute(poi.ModelAttribute.ObjectiveValue)
     assert obj_val == approx(N**2)
+
+
+# reported by https://github.com/metab0t/PyOptInterface/issues/59
+def test_shuffle_qp_objective(model_interface):
+    model = model_interface
+
+    N = 3
+    weights = model.add_m_variables(N, lb=0)
+
+    expected_returns = [0.05, 0.07, 0.12]
+    min_return = 0.06
+    cov = [
+        (0, 0, 0.1),
+        (0, 1, 0.1),
+        (0, 2, 0.04),
+        (1, 1, 0.4),
+        (1, 2, 0.2),
+        (2, 2, 0.9),
+    ]
+    cov_objs = [weights[i] * weights[j] * v for i, j, v in cov]
+
+    model.add_linear_constraint(poi.quicksum(weights) == 1)
+    model.add_linear_constraint(
+        poi.quicksum(expected_returns[i] * weights[i] for i in range(N)) >= min_return
+    )
+
+    trial = 120
+    obj_values = []
+    for _ in range(trial):
+        import random
+
+        random.shuffle(cov_objs)
+        obj = poi.quicksum(cov_objs)
+        model.set_objective(obj, poi.ObjectiveSense.Minimize)
+
+        model.optimize()
+
+        obj_value = model.get_model_attribute(poi.ModelAttribute.ObjectiveValue)
+        obj_values.append(obj_value)
+
+    obj_values = np.array(obj_values)
+    # test all values are the same
+    assert np.all(np.abs(obj_values - obj_values[0]) < 1e-8)