ENH: Normalize the payoff scale in _initialize_tableaux so lemke_howson's tolerances are effectively relative

Follow-up deferred from #949, where @oyamad wrote: "Payoff normalization in _initialize_tableaux was considered here and deferred, as it changes what the tolerances mean; the 1e-10 row of the sweep is a good argument for revisiting it."

The problem

_lex_min_ratio_test compares tableau entries against absolute tolerances — TOL_PIV = 1e-10 decides whether a pivot-column entry counts as positive, and TOL_RATIO_DIFF = 1e-15 decides whether two ratios are tied (pivoting.py:10-11). But _initialize_tableaux writes raw payoff values into the tableau, shifted by abs(min) + 1 when the minimum is non-positive but never rescaled (lemke_howson.py:295-306). So the tolerances mean something completely different depending on the units the payoffs happen to be written in, even though the equilibria themselves are invariant to a positive affine transformation of each player's payoffs.

Since #949 landed, most of the damage shows up as an honest converged=False rather than a silent wrong answer, which is a big improvement. It is still a case where lemke_howson declines to answer a question it could answer, and at the extreme end (1e15) wrong answers with converged=True survive.

Evidence

30 random 4x4 games per payoff scale (default_rng(0)), each answer validated against the equivalent unit-scale game via is_nash, so the absolute tolerance inside is_nash doesn't flatter tiny payoffs. Run against main at e504433 (i.e. with #949 included). "Normalized" applies a per-player positive affine transform — shift to positive, then divide by the maximum — before calling lemke_howson, which is exactly what doing the work inside _initialize_tableaux would achieve:

Normalization recovers a correct equilibrium in every case, including the twelve games at 1e15 that currently come back wrong with converged=True.

Proposed change

_initialize_tableaux already shifts each player's payoffs by abs(min) + 1 to keep the tableau non-negative with no identically-zero column. The suggestion is to divide by the maximum of the shifted values as well, so that the payoff block of each tableau lands in a fixed range and the absolute tolerances become effectively relative to the payoff scale. The right-hand side is already 1 for every row, so the two sides of each ratio would then be on comparable scales.

Caveats and open questions

• This changes what the tolerances mean, which is precisely why it was deferred from #949 rather than folded into it. It should be a PR of its own, with the sweep above (or a wider version of it) as evidence.
• Every user-visible answer for well-scaled games ought to be unchanged, but that needs demonstrating rather than assuming — the full game_theory and optimize suites plus the Netlib benchmark from #949 are the obvious checks.
• Worth deciding whether TOL_PIV / TOL_RATIO_DIFF should be revisited once the tableau scale is pinned down; the current values were chosen against unnormalized tableaux.
• _initialize_tableaux_ig in _compute_fp.py:286 builds its tableau from squared distances with the same shift-only treatment, so compute_fixed_point inherits the same sensitivity to the coordinate scale of X and Y. Whether it gets the same treatment in the same PR or a later one is worth a decision.
• There is a mirror in QuantEcon.jl, which took the #949 changes as QuantEcon/QuantEcon.jl#406; a port would keep the two implementations in step.

Reproducing the sweep

import numpy as np
from quantecon.game_theory import NormalFormGame, lemke_howson

def normalize(payoffs):
    out = payoffs.copy()
    for pl in range(2):
        P = out[:, :, pl]
        m = P.min()
        if m <= 0:
            P = P - m + 1
        out[:, :, pl] = P / np.abs(P).max()
    return out

for s in [1e-12, 1e-10, 1.0, 1e12, 1e15]:
    rng = np.random.default_rng(0)
    counts = {'as-is': [0, 0], 'normalized': [0, 0]}
    for _ in range(30):
        base = rng.random((4, 4, 2))
        g_unit = NormalFormGame(base)
        for key, payoffs in (('as-is', base * s), ('normalized', normalize(base * s))):
            NE, res = lemke_howson(NormalFormGame(payoffs), full_output=True)
            counts[key][0] += res.converged
            counts[key][1] += res.converged and g_unit.is_nash(tuple(NE))
    print(f"{s:8.0e}  as-is {counts['as-is']}  normalized {counts['normalized']}")

Other follow-ups from #949

Two further ideas were raised there and are deliberately not part of this issue, since they are independent and the second needs a design discussion of its own:

• A RuntimeWarning on non-convergence in lemke_howson / polym_lcp_solver, so that the breakdown reporting added in #949 reaches callers using the default full_output=False.
• Zeroing x and lambd on a Phase 1 failure in linprog_simplex, where they currently come back as allocated by np.empty.