Avoid normal equations in inverse least-squares updates

[nicamo]

Summary

Replace the unconstrained inverse solver's regularised normal equations,

(A.T @ A + R) delta = A.T @ b,

with the equivalent augmented least-squares problem,

[A; sqrt(R)] delta ~= [b; 0],

solved using numpy.linalg.lstsq.

Motivation

Forming A.T @ A squares the nonzero condition number of the response matrix and obscures rank deficiency. This matters for isoflux constraints because N points generate N(N-1)/2 pairwise residual rows, although those rows contain at most N-1 independent flux differences.

Representative FreeGSNKE matrices confirmed the redundancy:

• unconstrained test machine: A has shape (32, 10) and rank 6;
• MAST-U-like case: A has shape (32, 11) and rank 7;
• the nonzero condition number of A is only about 30, but the unregularised A.T @ A has a condition number of order 10^17.

The current default regularisation protects these two cases, but low or zero regularisation can make the normal-equation update inaccurate without an explicit failure. In a controlled pairwise-isoflux system with cond(A) = 10^8, the relative solution errors were:

Solver	Relative error
Normal equations	5.6e-1
Augmented least squares	6.2e-10

Uniform complete pairwise differences do not themselves increase the nonzero condition number; they expose rank redundancy. The avoidable numerical loss comes from subsequently forming the normal equations.

Changes

• Add a small private helper that solves diagonally regularised least-squares systems in augmented form.
• Validate that regularisation values are finite, non-negative, and have one entry per variable.
• Use the helper for unconstrained frozen-plasma coil-current updates.
• Use the helper for unconstrained full-plasma finite-difference Jacobian updates.
• Use the helper for the two-parameter plasma-profile update.
• Leave all constrained CVXPY paths unchanged.
• Document the augmented formulation and regularisation convention in solver docstrings.
• Update Example 1a, which feeds the generated user guide, to distinguish unconstrained augmented least squares from constrained quadratic optimisation.

For diagonal R, appending sqrt(R) preserves the existing objective exactly:

||A delta - b||^2 + delta.T R delta.

With zero regularisation and a rank-deficient matrix, lstsq returns the well-defined minimum-norm solution.

Validation

Focused tests cover:

• equivalence to normal equations on a well-conditioned regularised problem;
• accuracy for a pairwise system with cond(A) > 10^7;
• minimum-norm behavior for an exactly rank-deficient system;
• rejection of negative, infinite, and NaN regularisation.

Results:

6 focused tests passed
16 repository tests passed, 4 skipped

The MAST-U-like diverted inverse case was also run without current limits to force the unconstrained paths through the augmented solver. It converged successfully in 28 iterations with finite coil currents and flux.

The updated Example 1a notebook was also parsed and converted to HTML successfully.

Integration with current main

Current main was merged in 2dadc31. This includes the up-down symmetry correction from #101. The two changes combine without a behavioural conflict: #101 constrains the plasma states used by the inverse and full-Jacobian forward solves, while this PR changes only how their unconstrained least-squares update is calculated. The combined branch retains both the symmetry propagation and the augmented least-squares solve.