[xihaian251]

What this fixes

ot.solve documents that unbalanced=float("inf") recovers balanced OT and
that (inf, rho) / (rho, inf) request the semi-relaxed problem. Only the
Sinkhorn KL path implemented that sentinel. The exact (MM) and L-BFGSB paths
treated inf as a number and failed silently:

• ot.solve(M, a, b, reg=0, unbalanced=float("inf")) returned
  value=nan, plan=nan, status=None, with no warning.
• ot.solve(M, a, b, reg=0.5, reg_type="l2", unbalanced=(inf, 1.0)) returned
  value=nan together with a plan satisfying both marginals to 2.8e-17,
  i.e. it silently solved the balanced problem instead of the semi-relaxed one.

Root cause: mm_unbalanced builds r1 = reg_m1 / (reg + reg_m1 + reg_m2)
(inf / inf = nan), and lbfgsb_unbalanced evaluates reg_m1 * kl_div(...)
(inf * 0 = nan).

Why a plain "treat inf as balanced" is not enough

A two-sided infinite relaxation is the hard-constraint limit of the same
objective. It equals the problem solved by unbalanced=None only when the
balanced branch minimizes the same objective. That is not the case for every
regularizer. On the 2x2 fixture a = b = (0.5, 0.5), M = [[0, 1], [1, 0]]:

reg	reg_type	hard-limit value	unbalanced=None value
1	kl, default c = a b^T	0.3798855	0.3798855 (same)
1	kl, custom c (entropic path)	0.0652983 (x*=0.4804)	0.3798855 (custom c is ignored)
1	l2	0.125	0.5
2	l2	0.25	1.0
1	entropy	-2.0064	-1.0064 (off by exactly reg)

The custom reference measure used in that row is
c = [[0.45, 0.05], [0.05, 0.45]], which is feasible for these marginals;
the hard-constraint limit is therefore not c once M != 0
(x* = 0.4804, while c corresponds to x = 0.45).

Two independent reasons behind those differences:

• l2: the balanced branch solves reg/2 * ||P||^2 with no reference
  measure and reports <P,M> + reg * ||P||^2, while the UOT limit is
  <P,M> + reg/2 * ||P - c||^2.
• entropy: the balanced value omits -sum(P), so it differs by exactly reg.

The plans are not guaranteed to match either. With the default c = a b^T,
<P,c> is constant on the feasible set only when a or b is uniform. For
a = (0.7, 0.3), b = (0.6, 0.4), M = 0, reg = 1 the UOT limit is
P00 = 0.42 while the balanced L2 solver returns P00 = 0.40.

Dispatching those combinations would therefore return the solution of a
different optimization problem, so this PR rejects them instead.

What this changes

• reg == 0 with both relaxations infinite → the exact balanced solver
  (unbalanced=None), for any c (the regularizer vanishes).
• reg > 0 with reg_type="KL" and the default reference measure c = a b^T
  → the balanced KL solver; its objective is exactly the limit.
• A user supplied (f, df) pair → the balanced CG solver, which minimizes the
  same f.
• The entropic KL path is left untouched: it already implements the sentinel,
  honours a custom c when unbalanced is infinite, and is the only path that
  supports potentials_init and method. Note that the plain balanced KL path
  (unbalanced=None) ignores c and always uses a b^T.
• One-sided infinite relaxations on the exact / L-BFGSB paths, and two-sided
  ones for l2 / entropy / KL with a custom c outside the entropic KL
  path, raise an explicit NotImplementedError before any optimization.
  No silent NaN, no silently-wrong plan.
• mm_unbalanced / lbfgsb_unbalanced docstrings state that the wrapper
  resolves or rejects the sentinel while the low-level solvers still require a
  finite reg_m.
• RELEASES.md entry added.

Tests

test/unbalanced/test_inf_marginal.py, 99 cases over NumPy / PyTorch / JAX.
The 2x2 fixtures make the hard-constraint feasible set one dimensional
(P = [[x, .5-x], [.5-x, x]]), so every accepted combination is checked
against an independent analytic oracle (bounded scalar minimization of the
original objective), including the value - a value-only guard would miss the
entropy and l2 cases. The rejected combinations are covered by
exception-only tests; their numeric divergence is documented above, not
asserted in CI.

• new file: 97 passed, 2 skipped
• test/unbalanced + test/test_solvers.py: 1593 passed, 131 skipped
• the same test file against the previous (rejected) implementation: 30 failed
• patch applied to a pristine 5cb82c9: test/unbalanced 681 passed, 33 skipped

All numbers above are local runs (Windows, Python 3.13, NumPy 2.4.6,
SciPy 1.16.3, PyTorch 2.10 CPU, JAX 0.12); GitHub Actions results will
appear separately and are not claimed here.

Compatibility

Inputs that previously returned NaN or a wrong plan now either return the
correct limit or raise. Already-correct calls are unchanged, in particular the
entropic KL path with potentials_init / method. The only behaviour change
for a previously "successful" call is that l2 / entropy / KL + custom c
with unbalanced=inf now raise instead of returning the balanced solution.

Not verified: CUDA / CuPy / TensorFlow (no environment available).

Open questions for maintainers

1. Should entropy be supported instead of rejected? The plan coincides with
   the balanced one; only the reported value differs, by exactly reg. This
   could be handled by compensating the value, if you prefer that convention.
2. Independent pre-existing issue, not touched here: the balanced l2 branch
   reports reg * ||P||^2 while its solver minimizes reg/2 * ||P||^2.

[codecov]

Codecov Report

❌ Patch coverage is 96.49123% with 6 lines in your changes missing coverage. Please review.
✅ Project coverage is 97.00%. Comparing base (5cb82c9) to head (07ee2b3).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #895      +/-   ##
==========================================
- Coverage   97.00%   97.00%   -0.01%     
==========================================
  Files         130      131       +1     
  Lines       26842    27012     +170     
==========================================
+ Hits        26039    26203     +164     
- Misses        803      809       +6     

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.