ot.solve(..., reg_type="L2") reports a value that is not the optimum of the solved problem #896

Summary

For balanced quadratic-regularized OT, ot.solve returns the optimal plan of

$$ \min_{P\ge 0}\ \langle M,P\rangle+\frac{\lambda}{2}|P|_F^2 \quad\text{s.t.}\quad P\mathbf 1=a,\ P^\top \mathbf 1=b, $$

but reports

$$ \texttt{res.value}=\langle M,P\rangle+\lambda|P|_F^2 =J(P)+\frac{\lambda}{2}|P|_F^2 . $$

res.plan and res.value_linear are correct; res.value is not the optimal value of the problem that was solved. In the exact example below, res.value = 3.0 while the true optimum is 1.7.

User impact: any downstream use of res.value as a loss, a benchmark number, or a model-selection criterion is affected; solve_sample and solve_bary_sample inherit the same scalar. Transport plans are not affected.

Minimal reproduction

import numpy as np
import ot

print(ot.__version__, ot.__file__)  # 0.9.8.dev / checkout

M = np.array([[0.0, 1.0], [1.0, 0.0]])
a = np.array([0.5, 0.5])
b = np.array([0.5, 0.5])
reg = 10.0

res = ot.solve(M, a, b, reg=reg, reg_type="L2")
P = np.asarray(res.plan)
J = np.sum(M * P) + 0.5 * reg * np.sum(P ** 2)

print(P)                 # [[0.3 0.2]
                         #  [0.2 0.3]]
print(res.value_linear)  # 0.4
print(res.value)         # 3.0
print(J)                 # 1.7
print(res.value - J)     # 1.3 == 0.5 * reg * ||P||^2

Actual: res.value = 3.0, res.value_linear = 0.4, res.plan = [[0.3, 0.2], [0.2, 0.3]].

Expected: the plan is correct; res.value should be the optimal value 17/10 = 1.7 (or, if the maintainers intend the alternative convention discussed below, the plan should be a different matrix and the value should be 2.975). It cannot be 3.0 in either case.

The plan is the exact KKT solution for this instance (the solver returns [[3/10, 1/5], [1/5, 3/10]] to machine precision), so this is not a convergence artifact.

Mathematical basis

• ot/smooth.py:SquaredL2.Omega defines $\Omega(T)=\frac{\gamma}{2}|T|^2$; delta_Omega(X)=\sum[\max(X,0)]^2/(2\gamma), and get_plan_from_dual returns $P=[\alpha+\beta-M]_+/\gamma$. With gamma=reg, the solver called by ot.solve minimizes the half-squared primal.
• On the support of the returned plan, $\alpha_i+\beta_j=M_{ij}+\lambda P_{ij}$, hence

$$ a^\top\alpha+b^\top\beta=\langle M,P\rangle+\lambda|P|^2=\texttt{res.value}. $$

In other words, the reported value is the affine part of the L2 dual; the conjugate term $\sum_j\delta_\Omega(X_{:j})=\frac{\lambda}{2}|P|^2$ is omitted. The dual objective itself is 1.7.

• This matches the cited paper (Blondel, Seguy & Rolet, AISTATS 2018, Table 1) and POT's other L2 paths (reg_div="l2" is documented and implemented as half-squared), and it is bit-identical to calling ot.smooth.smooth_ot_dual(a, b, M, reg=reg) directly.

Why this is not only a documentation/convention question

Even if one redefines the intended objective as $\langle M,P\rangle+\lambda|P|^2$ (the no-half formula currently in docs/source/user_guide.rst), the current outputs are still not consistent: the returned plan is the half-squared minimizer, not the no-half minimizer, and the reported value is the no-half objective at that suboptimal plan. For the same instance:

• half-squared optimum: 1.7, at P=[[3/10,1/5],[1/5,3/10]] (returned plan);
• no-half optimum: 2.975 = 119/40, at P=[[11/40,9/40],[9/40,11/40]];
• reported value: 3.0, which is neither optimum (and is 0.025 above the no-half optimum).

So no convention makes both plan and value correct. The plan/value pair has to be made consistent one way or the other.

A note on visibility: when the two problems happen to share a minimizer (sparse/vertex plans, small regularization, or cost structures that force a vertex), the reported value can coincide numerically with the no-half optimum. That does not make the half-squared value correct: it is still larger than the half-squared optimum by exactly λ/2‖P‖², because sharing a plan does not mean the two objectives have the same value. In that regime the bug is hidden only if one interprets the value under the no-half convention; when regularization genuinely moves the plan, the no-half value is also visibly non-optimal.

A second exact instance isolates the value error from the convention question (there the two conventions happen to share the minimizer):

M=[[1,2],[3,4]], a=b=(1/2,1/2), reg=1
plan = 1/4 * ones(2,2)          (correct)
true J_half = 21/8 = 2.625
reported value = 2.7500008      (+0.125 = reg/2*||plan||^2)

A third exact instance ($M=\begin{pmatrix}0&1\1&0\end{pmatrix}$, uniform marginals, reg=4) gives plan [[3/8,1/8],[1/8,3/8]], true optimum 7/8, reported 3/2, and the no-half optimum 23/16 = 1.4375; again 3/2 is neither.

Scope

• ot.solve(..., reg_type="L2") and delegates: ot.solve_sample(..., reg_type="L2") and ot.solve_bary_sample(..., reg_type="L2") (the latter also uses the value for stopping_criterion="loss"). The solve_sample and solve_bary_sample propagation was verified in small runs; the barycenter support itself was unchanged because the BCD update uses the plans, while the reported loss is inflated.
• Affected releases: 0.9.0 through 0.9.7.post1, plus master 0.9.8.dev (the line was introduced in PR [MRG] New API for OT solver (with pre-computed ground cost matrix) #388, 2022-12-15).
• Backends: verified on NumPy, PyTorch and JAX; the expression is backend-generic (TensorFlow/CuPy not tested here).
• The transport plans are not affected; only the reported value (and the barycenter loss criterion).

Documentation issue

docs/source/user_guide.rst defines the quadratic regularizer as $\Omega(\gamma)=\sum\gamma^2$ and points the reader to ot.smooth.smooth_ot_dual, whose Omega() is $\frac{\gamma}{2}\sum\gamma^2$. The ot.solve docstring writes $\lambda R(T)$ without defining $R$. This contradiction is likely why the value line survived four years: the solver and the value line encode two different conventions.

Suggested fix

Minimal, non-breaking:

-                value = value_linear + reg * nx.sum(plan**2)
+                value = value_linear + 0.5 * reg * nx.sum(plan**2)

plus a statement in the ot.solve/ot.solve_sample/ot.solve_bary_sample docstrings that $R(T)=\frac12|T|_F^2$, and a correction of the user-guide formula to $\Omega(\gamma)=\frac12\sum\gamma^2$. This keeps every plan unchanged and makes value the true optimal value; it matches ot.smooth, the cited paper, and POT's "half-squared" unbalanced convention. I can open a PR with this change and a regression test (exact rational instances, NumPy/PyTorch/JAX) if maintainers agree.

The alternative — treating $\Omega=\sum\gamma^2$ as authoritative and passing 2*reg to smooth_ot_dual — would change every L2 plan returned since 0.9.0 and break the current bit-identical agreement with ot.smooth.smooth_ot_dual. That looks like a versioned behavior change rather than a bug fix, so I would prefer to confirm the intended convention before touching the plan.

Relation to existing work

PR #895 (open) notes this as an “independent pre-existing issue” in its open questions but deliberately does not fix it. This issue is meant to track it separately and provide the exact counterexample and convention analysis. PR #874 (merged) is the precedent for fixing the reported value of an OT solver without changing plans.

Environment

Python 3.13.1, numpy 2.4.6, scipy 1.16.3, POT master 5cb82c9 (0.9.8.dev)
Plan and value independently checked with exact rational KKT enumeration (SymPy),
CVXPY (Clarabel), and Wolfram Language Minimize.