Document augmented inverse least squares

Author: nicamo

examples/example01a - static_inverse_solve_MASTU.ipynb

@@ -387,19 +387,23 @@
                 "\n",
                 "During an inverse solve, a minimisation problem involving the responses, changes in coil currents, and constraints, is repeatedly solved:\n",
                 "\n",
-                "$$ \\min_{x} \\| A\\vec{x} - \\vec{b}\\|^2 + \\| \\gamma \\vec{x} \\|^2, $$\n",
+                "$$ \\min_{x} \\| A\\vec{x} - \\vec{b}\\|^2 + \\vec{x}^{T} \\Gamma \\vec{x}, $$\n",
                 "\n",
                 "where\n",
                 "- $A$ is the fixed response matrix (that determines how a change in coil currents $x$, affects constraint values). \n",
                 "- $\\vec{x} = \\Delta \\vec{I}^c$ is the step change in active coil currents required to match the constraints. \n",
                 "- $\\vec{b}$ is the vector of constraint values being enforced. \n",
-                "- $\\gamma > 0$ is the regularisation parameter/vector.\n",
+                "- $\\Gamma$ is the non-negative diagonal regularisation matrix: a scalar `l2_reg` gives $\\Gamma = \\mathtt{l2\\_reg} I$, while a vector supplies its diagonal entries.\n",
                 "\n",
-                "We solve for $\\vec{x}$ using a gradient-based optimiser. \n",
+                "Without inequality constraints, FreeGSNKE solves the equivalent augmented system\n",
+                "\n",
+                "$$ \\min_x \\left\\| \\begin{bmatrix} A \\\\ \\Gamma^{1/2} \\end{bmatrix} \\vec{x} - \\begin{bmatrix} \\vec{b} \\\\ 0 \\end{bmatrix} \\right\\|^2, $$\n",
+                "\n",
+                "using a direct least-squares solve. This avoids forming the normal equations $(A^T A + \\Gamma)\\vec{x} = A^T\\vec{b}$. When coil-current or normalised-flux limits are active, the same regularised objective is handled by the constrained quadratic optimiser instead.\n",
                 "\n",
                 "[Song et al. (2024)](https://www.mdpi.com/2571-6182/7/4/45) provide a nice overview of the inverse problem.\n",
                 "\n",
-                "The coil limits are handled in a slightly different manner because they are box constraints that cannot be encoded in the minimisation problem. Instead, FreeGSNKE adds additional constraints to the solver that bound the solution $\\Delta \\vec{I}^c$. These bounds have a tolerance (slack) which allows the coils to violate their limits, but penalises violations; this allows the coil limits to be violated 'on the path' to a solution that respects the prescribed limits.\n"
+                "The coil limits are handled in a slightly different manner because they are box constraints that cannot be represented by the unconstrained augmented least-squares system. Instead, FreeGSNKE adds constraints to the quadratic optimiser that bound the solution $\\Delta \\vec{I}^c$. These bounds have a tolerance (slack) which allows the coils to violate their limits, but penalises violations; this allows the coil limits to be violated 'on the path' to a solution that respects the prescribed limits.\n"
             ]
         },
         {
@@ -414,7 +418,7 @@
                 "\n",
                 "The `verbose=True` option will provide an indication of the progression of the solve. \n",
                 "\n",
-                "The `l2_reg` parameter defines the Tikhonov regularisation (i.e. $\\gamma$ in the cell above) used by the optimiser. This can be set as a scalar or as a vector (equal to the number of coil currents being solved for). Larger values force coil currents to stay closer to their original values while lower values will encourage more freedom. This can be useful in particular for vertical control coils (e.g. the P6 coil in MAST-U), in which we don't want the coil current to \"jump around\" during optimisation and cause vertical instability issues. \n",
+                "The `l2_reg` parameter defines the diagonal of the Tikhonov regularisation matrix $\\Gamma$ in the cell above. It can be set as a scalar or as a vector (equal to the number of coil currents being solved for). Larger values force coil currents to stay closer to their original values while lower values encourage more freedom. This can be useful in particular for vertical control coils (e.g. the P6 coil in MAST-U), in which we don't want the coil current to \"jump around\" during optimisation and cause vertical instability issues. \n",
                 "\n",
                 "\n",
                 "The solver steps are (roughly):\n",

freegsnke/GSstaticsolver.py

@@ -1232,14 +1232,15 @@ def optimize_currents(
         Once A is constructed, the Newton step is computed by solving the
         Tikhonov-regularised least-squares problem:
 
-            min || A ΔI + b0 ||² + ||R ΔI||²
+            min || A ΔI + b0 ||² + ΔIᵀ R ΔI
 
         where:
             b0 = current constraint residual
             R  = regularisation matrix
 
-        If current or flux limits are active, a quadratic optimisation
-        routine is used instead of the closed-form normal equations.
+        Without inequality limits, this is solved as an augmented
+        least-squares system. If current or flux limits are active, a
+        constrained quadratic optimisation routine is used instead.
 
         Parameters
         ----------

freegsnke/inverse.py

@@ -27,7 +27,7 @@
 
 
 def _solve_regularized_lstsq(A, b, regularization):
-    """Solve a diagonally regularized least-squares problem without normal equations."""
+    """Solve ``min ||A x - b||² + xᵀ diag(regularization) x`` in augmented form."""
     regularization = np.asarray(regularization)
     if regularization.shape != (A.shape[1],):
         raise ValueError(
@@ -1213,13 +1213,13 @@ def optimize_currents(
 
         This method computes optimal coil current corrections by solving:
 
-            min_I || A I − b ||² + λ || I ||²
+            min_I || A I − b ||² + Iᵀ R I
 
         where:
 
             A = combined constraint Jacobian matrix
             b = combined constraint residual vector
-            λ = Tikhonov (L2) regularisation parameter
+            R = diagonal Tikhonov regularisation matrix
 
         The optimisation accounts for:
 
@@ -1252,10 +1252,10 @@ def optimize_currents(
             Tikhonov regularisation parameter.
 
             If float:
-                λ I² penalty is applied uniformly.
+                R is the scalar value times the identity matrix.
 
             If array:
-                Allows coil-wise regularisation weighting.
+                R contains the supplied coil-wise values on its diagonal.
 
         Returns
         -------