Resolve trust-constr vs linear_sum_assignment mismatches by relaxing strict equality constraints in SciPy

Learn why SciPy trust-constr diverges from linear_sum_assignment in the assignment problem, and how relaxing equality constraints restores accuracy reliably.

When you compare the optimal assignment from linear_sum_assignment with a solution from minimize on the same linear cost and constraints, you expect identical answers. In practice, adding an analytical Jacobian and Hessian to trust-constr can make the two diverge. Without derivatives the error stays around 1.5e-7, with only a Hessian it grows to roughly 8e-4, and with both Jacobian and Hessian it can exceed 0.07. The twist is that the derivatives are fine; the strict equality constraints are the part that needs a small adjustment.

Reproducing the setup

The example below sets up the linear assignment as a vectorized linear program, solves it once with linear_sum_assignment for a reference solution, and then calls minimize with trust-constr, linear constraints, bounds, and supplied derivatives.

import numpy as np
from scipy.optimize import linear_sum_assignment, minimize, LinearConstraint
# weight matrix for the assignment
wmat = np.array(((54, 54, 51, 53), (51, 57, 52, 52), (50, 53, 54, 56), (56, 54, 55, 53)))
# problem size
m = 4
m2 = m * m
# reference optimal assignment via Hungarian algorithm
r_idx, c_idx = linear_sum_assignment(wmat)
Z = np.zeros_like(wmat)
Z[r_idx, c_idx] = 1
# linear objective and its derivatives in vector form
def obj(v, q):
    return q.dot(v)
def d_obj(v, q):
    return q
def h_obj(v, q):
    return np.zeros((len(q), len(q)))
qvec = wmat.flatten().astype(float)
# equality constraints: row sums == 1 and column sums == 1
G = np.zeros((2 * m, m2))
for j in range(m):
    G[j, j * m:(j + 1) * m] = 1
    G[m + j, j::m] = 1
unit = np.ones(2 * m)
aff = LinearConstraint(G, lb=unit, ub=unit)
# bounds and starting point
box = [(0, 1)] * m2
v0 = np.eye(m).flatten()
# trust-constr with analytical Jacobian and Hessian
res1 = minimize(
    obj, v0, args=(qvec,), constraints=[aff], bounds=box,
    method='trust-constr', jac=d_obj, hess=h_obj
)
Y = res1.x.reshape((m, m))
print(np.abs(Y - Z).max())

What is actually going on

The objective is linear, the constraints are linear, and the provided Jacobian and Hessian are consistent with that formulation. The discrepancy does not stem from incorrect derivatives. The key observation is that making the equality constraints slightly elastic resolves the mismatch. There is no definitive explanation offered beyond that; it is enough to know that the strict equalities are sensitive in this setup and that a small tolerance around them helps trust-constr behave well.

The fix: relax equalities into tight inequalities

Instead of enforcing row and column sums to be exactly one, allow them to deviate by an epsilon and tighten solver tolerances. This change aligns the trust-constr solution with the Hungarian result to very high accuracy, around 1e-12.

import numpy as np
from scipy.optimize import linear_sum_assignment, minimize, LinearConstraint
wmat = np.array(((54, 54, 51, 53), (51, 57, 52, 52), (50, 53, 54, 56), (56, 54, 55, 53)))
m = 4
m2 = m * m
r_idx, c_idx = linear_sum_assignment(wmat)
Z = np.zeros_like(wmat)
Z[r_idx, c_idx] = 1
def obj(v, q):
    return q.dot(v)
def d_obj(v, q):
    return q
def h_obj(v, q):
    return np.zeros((len(q), len(q)))
qvec = wmat.flatten().astype(float)
G = np.zeros((2 * m, m2))
for j in range(m):
    G[j, j * m:(j + 1) * m] = 1
    G[m + j, j::m] = 1
unit = np.ones(2 * m)
# relaxed equalities
delta = 1e-15
aff_loose = LinearConstraint(G, lb=unit - delta, ub=unit + delta)
box = [(0, 1)] * m2
v0 = np.eye(m).flatten()
res2 = minimize(
    obj, v0, args=(qvec,), constraints=[aff_loose], bounds=box,
    method='trust-constr', jac=d_obj, hess=h_obj,
    options=dict(gtol=1e-13, xtol=1e-13)
)
Y2 = res2.x.reshape((m, m))
print(np.abs(Y2 - Z).max())

Why this matters

Even for a fully linear assignment with a linear objective, switching from an exact equality to a near-equality can be the difference between visibly wrong and near machine-precision results when using trust-constr with supplied derivatives. If you are validating a nonlinear pipeline on a linear problem first, this is an easy and effective adjustment that removes an unnecessary source of error from the comparison.

SciPy has a solver for linear programming problems and the assignment problem can be formulated as a linear program. This is usually faster and more accurate than non-linear methods. But perhaps your real problem has nonlinearity.

Takeaways

When trust-constr with an analytical Jacobian and Hessian disagrees with linear_sum_assignment on a linear assignment formulation, keep the same linear objective and constraints but relax the equalities into tight inequalities using a very small epsilon and set stricter tolerances. If your production problem remains linear, prefer linear_sum_assignment or a linear programming approach. If you are preparing for nonlinear extensions, this small relaxation is a practical way to get minimize to match the known optimum.