Native heap corruption in v1.14.0 when passHessian receives a square Hessian with epsilon-level asymmetry (via CVXPY/highspy)

[raphaelsaavedra]

Summary

When a square Hessian with floating-point-level asymmetry (on the order of machine epsilon, ~1e-18) is passed to HiGHS 1.14.0 via passHessian from the C API / highspy, HiGHS first emits ERROR: Square Hessian contains N non-symmetries and then the Python process crashes with native heap corruption (malloc: Incorrect checksum for freed object, free(): invalid size, double free or corruption, munmap_chunk(): invalid pointer, SIGABRT / SIGSEGV). The crash is not catchable from Python since the entire interpreter is terminated.

From #2821 / #2857: passHessian now returns kError on any asymmetry instead of using (Q+Q^T)/2. The check at highs/model/HighsHessianUtils.cpp#L422 uses strict != equality with no tolerance, so two values that differ by 1e-18 are counted as asymmetric.

The practical impact is that CVXPY's HiGHS backend fails to solve certain QPs whose Hessian has dense off-diagonal entries. In our case, the Hessian originates from cp.quad_form(w, Sigma) where Sigma is the result of a PSD projection (V @ diag(λ) @ V.T) of a covariance matrix. Although Sigma is mathematically symmetric, the matrix that ultimately reaches HiGHS is not bit-exactly symmetric after floating-point reconstruction and CVXPY's canonicalization, and the new strict equality check flags all such asymmetries as errors.

Environment

• HiGHS: 1.14.0 (via highspy 1.14.0, PyPI wheels)
• CVXPY: 1.7.5
• Python: 3.12
• NumPy 1.26.4, SciPy 1.15.3
• Reproduced on both macOS arm64 (local) and Linux x86_64 (containers). Reproduced sequentially, so not a threading issue.

Reproduction

import numpy as np
import cvxpy as cp
import scipy.sparse as sp

rng = np.random.default_rng(0)

N_PATHS = 150
N_NODES = 50

# ---- Sparse "incidence" matrix mapping paths to nodes ----
rows, cols, vals = [], [], []
for p in range(N_PATHS):
    src, snk = rng.choice(N_NODES, size=2, replace=False)
    rows += [p, p]
    cols += [src, snk]
    vals += [1.0, -1.0]

incidence = sp.csr_matrix((vals, (rows, cols)), shape=(N_PATHS, N_NODES))

# ---- Dense PSD covariance on nodes, built by eigenvalue clamping ----
# (mirrors a typical "project onto PSD cone" step used after estimating a
#  covariance matrix from noisy / partially-missing data)
A = rng.standard_normal((N_NODES, N_NODES))
Sigma_raw = A @ A.T + rng.standard_normal((N_NODES, N_NODES)) * 0.01
Sigma_sym = 0.5 * (Sigma_raw + Sigma_raw.T)
w_eig, V_eig = np.linalg.eigh(Sigma_sym)
w_eig = np.maximum(w_eig, 0.0)
Sigma = V_eig @ np.diag(w_eig) @ V_eig.T
# Sigma is mathematically symmetric and PSD, but not bit-exact symmetric:
print("max |Sigma - Sigma.T|:", np.max(np.abs(Sigma - Sigma.T)))
print("min eigenvalue:", float(np.linalg.eigvalsh(Sigma).min()))

# ---- Build the CVXPY problem ----
w_path = cp.Variable(N_PATHS, nonneg=True)

# Linear transformation to "nodal" space via the sparse incidence matrix.
w_nodal = w_path @ incidence  # shape: (N_NODES,)

# Some realistic constraints so the problem is bounded and not trivial.
path_coef = rng.standard_normal(N_PATHS)
path_upper = 50.0
constraints = [
    w_path <= path_upper,
    cp.sum(w_path) <= 500.0,
]

cov_penalty = 0.1
diag_penalty = 0.01
diag_weights = rng.uniform(0.1, 2.0, size=N_NODES)

obj_terms = [
    cp.sum(cp.multiply(path_coef, w_path)),
    -diag_penalty * cp.sum_squares(cp.multiply(diag_weights, w_nodal)),
    -cov_penalty * cp.quad_form(w_nodal, Sigma),
]
prob = cp.Problem(cp.Maximize(cp.sum(obj_terms)), constraints)

prob.solve(solver=cp.HIGHS, verbose=True)

results in

Running HiGHS 1.14.0 (git hash: 7df0786): Copyright (c) 2026 under MIT licence terms
ERROR:   Square Hessian contains 743 non-symmetries
QP has 401 rows; 250 cols; 1150 matrix nonzeros; 1325 Hessian nonzeros
Coefficient ranges:
  Matrix  [1e-01, 2e+00]
  Cost    [3e-03, 3e+00]
  Hessian [3e-03, 1e+01]
  Bound   [0e+00, 0e+00]
  RHS     [5e+01, 5e+02]
  Iteration        Objective     NullspaceDim
python(76500,0x1f2066140) malloc: Incorrect checksum for freed object 0x132c4f600: probably modified after being freed.
Corrupt value: 0xc004115f9637ddc9
python(76500,0x1f2066140) malloc: *** set a breakpoint in malloc_error_break to debug
[1]    76500 abort      python

[jajhall]

I could argue that it's for users to ensure that the redundant information in a square Hessian is consistent, but it's not hard to replace both values with their average if the difference is less than a small tolerance

[raphaelsaavedra]

Thanks, the tolerance-averaging would fix this.

I could argue that it's for users to ensure that the redundant information in a square Hessian is consistent

Agreed if you're using HiGHS directly, but many HiGHS users reach the solver through a modeling layer (CVXPY, JuMP) rather than constructing the Hessian directly, and those users have no natural place to enforce bit-exact symmetry. Even starting from a symmetric matrix, sparse-matrix canonicalization (e.g. A^T Σ A expansions) doesn't generally preserve bit-exact symmetry.

Separately, the crash itself (heap corruption and SIGABRT) looks like a bug worth tracking.

[jajhall]

Thanks, the tolerance-averaging would fix this.

I could argue that it's for users to ensure that the redundant information in a square Hessian is consistent

Agreed if you're using HiGHS directly, but many HiGHS users reach the solver through a modeling layer (CVXPY, JuMP) rather than constructing the Hessian directly, and those users have no natural place to enforce bit-exact symmetry. Even starting from a symmetric matrix, sparse-matrix canonicalization (e.g. A^T Σ A expansions) doesn't generally preserve bit-exact symmetry.

I can well believe this, and it's interesting to know that some interfaces provide a square Hessian (ie with redundant information that needs to be consistent - modulo tolerance discussed).

Separately, the crash itself (heap corruption and SIGABRT) looks like a bug worth tracking.

Yes, very odd. I guess I can simulate in C++ what you've observed.

[odow]

CVXPY could also fix this by passing the triangular Hessian. It's what we do in HiGHS.jl.

The issue with a square Hessian is "what exactly does symmetric mean".

Unrelated to HiGHS, but this also comes up in the context of PSD matrices, where we need to add constraints that x[i,j]==x[j,i] if they are not equal. JuMP has used 1e-10 as the approximate comparison for 10 years, and I don't think anyone has ever noticed or complained: jump-dev/JuMP.jl#976

It shows up here: https://github.com/jump-dev/MathOptInterface.jl/blob/692f8f2de733351eaf4336f8252361e27d76e29a/src/Bridges/Constraint/bridges/SquareBridge.jl#L107-L112

[jajhall]

In #2984 I've added a hard-coded tolerance (of 1e-10) on the difference between lower/upper triangular entries, but couldn't reproduce the crash that you observed. Note that if the Hessian is rejected due to asymmetry, then there will be no model in HiGHS if it was passed with its Hessian, or no Hessian if the model was passed as an LP and then a Hessian passed with passHessian. Return codes indicate this. Perhaps the lack of a model or Hessian leads to you crash if the return codes are not acted on appropriately.