NumPy 2 changed np.linalg.lstsq's default rcond: rank and solution shift
np.linalg.lstsq(A, b) without rcond used the old cutoff eps*s_max in NumPy 1.x, with a FutureWarning (1.26.4 signature rcond='warn'). NumPy 2.x defaults to rcond=None, i.e. eps*max(M,N)*s_max, so the warning is gone and the result can change. Measured with np.vander(linspace(0,1,50), 20) plus 1e-3 noise: rank 20 on 1.26.4 vs 19 on 2.5.4, coefficient norm 1.5e10 vs 8.2e8, and the fitted curves differ by 9.3e-4 (noise-level). At 26 columns the ranks were 23 vs 21 and the fits differed by 5.5e-3. It only matters when s_min/s_max lies between eps and eps*max(M,N), but then a NumPy upgrade silently changes regression output. Pass rcond explicitly (-1 reproduces 1.x on 2.x), or better, use a well-conditioned basis (np.polynomial.Polynomial.fit or Chebyshev on scaled x) so the cutoff never bites. np.polyfit is unaffected: it already defaulted to len(x)*eps.
Confirmed with NumPy 1.26.4, 2.0.0 and 2.5.4 on Python 3.12.3. The lstsq signature default is rcond='warn' on 1.26.4 (FutureWarning on every call without rcond) and None on both 2.0.0 and 2.5.4 (no warning), so the change is already in 2.0.0. With np.vander(np.linspace(0,1,50), k) and b = sin(3x) + 1e-3 * default_rng(seed).standard_normal(50), seeds 0, 1 and 2 all gave rank 20 vs 19 at k = 20 and 23 vs 21 at k = 26. Seed 0 also gave your coefficient norms (1.5e10 vs 8.2e8). On 2.5.4, rcond=-1 gives back 20 and 23. The polyfit remark also holds. polyfit(x, b, 19, full=True) reports rank 19 and rcond 1.11e-14 (50*eps) on all three versions. The fitted values differ by at most 6.2e-8 between 1.26.4 and 2.5.4, and by the same amount between 2.0.0 and 2.5.4, so that difference does not come from the lstsq default.