milanofthe · Jul 26, 2026 · Jul 26, 2026
diff --git a/src/pathsim/optim/anderson.py b/src/pathsim/optim/anderson.py
@@ -8,10 +8,11 @@
 # IMPORTS ==============================================================================
 
 import numpy as np
+import warnings
 
 from collections import deque
 
-from scipy.linalg import lu_factor, lu_solve
+from scipy.linalg import lu_factor, lu_solve, LinAlgWarning
 
 from .numerical import num_jac
 
@@ -24,6 +25,86 @@
     )
 
 
+# HELPERS ==============================================================================
+
+def _factorize(A):
+    """LU factorization of the Newton matrix, or 'None' if it is singular.
+
+    The Newton matrix is structurally singular whenever a state's derivative
+    does not depend on that state. A 'PID' is the canonical case: its integrator
+    state has ``dx/dt = u``, so the corresponding row of the Jacobian is zero.
+    The same holds for any block whose dynamics are a pure quadrature.
+
+    'scipy.linalg.lu_factor' does not raise on a singular matrix, it emits a
+    'LinAlgWarning' and leaves a zero pivot behind, after which 'lu_solve'
+    returns 'inf'/'nan'. Those propagate silently into the block state and only
+    surface many iterations later, in a different place, as an SVD failure
+    inside the anderson least squares. Detecting it here keeps the failure where
+    its cause is.
+
+    Parameters
+    ----------
+    A : np.ndarray
+        Newton matrix
+
+    Returns
+    -------
+    lu : tuple, None
+        factorization for 'lu_solve', or 'None' if 'A' is singular
+    """
+    if A.size == 0:
+        return None
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", LinAlgWarning)
+        _lu, _piv = lu_factor(A)
+
+    #a zero pivot on the diagonal of 'U' means singular, scaled to the matrix
+    _scale = np.max(np.abs(A))
+    if _scale == 0.0:
+        return None
+    if np.min(np.abs(np.diag(_lu))) <= np.finfo(float).eps * _scale * len(A):
+        return None
+
+    return _lu, _piv
+
+
+
+def _singular_step(A, res):
+    """Newton correction for a singular Newton matrix.
+
+    Newton is only defined on the range of 'A'. There the pseudo-inverse gives
+    the usual correction; on the null space the residual carries no information
+    about the state, so there is no correction to make and the entry is left to
+    the outer iteration.
+
+    Substituting a plain fixed-point step 'x <- g' on the null space was tried
+    and is worse: on a closed loop it injects the full residual every sweep and
+    the iteration diverges rather than stalling.
+
+    Parameters
+    ----------
+    A : np.ndarray
+        singular Newton matrix
+    res : np.ndarray
+        fixed-point residual 'g - x'
+
+    Returns
+    -------
+    dx : np.ndarray
+        correction, subtracted from the current iterate
+    """
+    _U, _s, _Vt = np.linalg.svd(A)
+    _tol = np.finfo(float).eps * (_s[0] if _s.size else 0.0) * len(A)
+    _ok = _s > _tol
+
+    #pseudo-inverse on the range of 'A'
+    _inv = np.where(_ok, 1.0 / np.where(_ok, _s, 1.0), 0.0)
+    _dx = _Vt.T @ (_inv * (_U.T @ res))
+
+    return _dx
+
+
 # CLASS ================================================================================
 
 class Anderson:
@@ -343,22 +424,44 @@ def _newton(self, x, g, jac):
 
         #early exit for scalar or purely vectorial values
         if _res.size == 1 or np.ndim(_jac) == 1:
-
-            return _x - _res / (_jac - 1.0), np.linalg.norm(_res)
+
+            _den = np.broadcast_to(
+                np.atleast_1d(np.asarray(_jac, dtype=float).flatten() - 1.0),
+                _res.shape
+                )
+
+            #a vanishing denominator is the scalar face of a singular Newton
+            #matrix: the residual does not depend on the state in that
+            #direction, so there is no Newton correction to make. An 'Integrator'
+            #hits this exactly ('dx/dt = u', Jacobian 0), and so does any block
+            #whose Jacobian passes through zero, e.g. 'dx/dt = -x**2 + u' at
+            #'x = 0'. Dividing anyway yields inf/nan which propagates through the
+            #connections into other blocks and only surfaces much later as an SVD
+            #failure inside the anderson least squares. Leave those entries to
+            #the anderson step, which is what drives them through the system.
+            #a vanishing denominator is the scalar face of a singular Newton
+            #matrix, handled the same way: no correction where Newton is
+            #undefined, see '_singular_step'.
+            _ok = np.abs(_den) > np.finfo(float).eps
+            _dx = np.where(_ok, _res / np.where(_ok, _den, 1.0), 0.0)
+
+            return _x - _dx, np.linalg.norm(_res)
 
         #vectorial values (newton raphson), Newton matrix 'jac - I'
         _A = _jac - np.eye(len(_res))
 
         #reuse the cached factorization while the Newton matrix is unchanged
         #(constant jacobian across iterations / stages), refactor otherwise
-        if self._A is not None and _A.shape == self._A.shape \
-            and np.array_equal(_A, self._A):
-            _lu = self._lu
-        else:
-            _lu = lu_factor(_A)
-            self._A, self._lu = _A.copy(), _lu
+        if self._A is None or _A.shape != self._A.shape \
+            or not np.array_equal(_A, self._A):
+            self._A, self._lu = _A.copy(), _factorize(_A)
+
+        #singular Newton matrix -> pseudo-inverse where Newton is defined,
+        #plain fixed-point where it is not
+        if self._lu is None:
+            return _x - _singular_step(_A, _res), np.linalg.norm(_res)
 
-        return _x - lu_solve(_lu, _res), np.linalg.norm(_res)
+        return _x - lu_solve(self._lu, _res), np.linalg.norm(_res)
 
 
     def step(self, x, g, jac=None):

diff --git a/tests/pathsim/optim/test_anderson.py b/tests/pathsim/optim/test_anderson.py
@@ -15,7 +15,9 @@
 from pathsim.optim.anderson import (
     Anderson,
     NewtonAnderson,
-    solve_root
+    solve_root,
+    _factorize,
+    _singular_step
     )
 
 
@@ -238,6 +240,78 @@ def test_warmstart_converges_immediately(self):
         _, _, it = solve_root(NewtonAnderson(), func, x, tolerance=1e-12)
         self.assertEqual(it, 0)
 
+class TestSingularNewtonMatrix(unittest.TestCase):
+    """
+    Test the handling of a singular Newton matrix.
+
+    The Newton matrix 'jac - I' is structurally singular whenever a state's
+    derivative does not depend on that state: a 'PID's integrator has
+    'dx/dt = u', so its Jacobian row is zero. 'scipy.linalg.lu_factor' does not
+    raise on that, it warns and leaves a zero pivot, after which 'lu_solve'
+    returns inf/nan. Those propagate into the block state and only surface many
+    iterations later as an SVD failure inside the anderson least squares.
+    """
+
+    def test_factorize_rejects_singular(self):
+        """a singular matrix yields no factorization"""
+
+        self.assertIsNone(_factorize(np.array([[-50.0, 0.0], [0.0, 0.0]])))
+        self.assertIsNone(_factorize(np.zeros((2, 2))))
+        self.assertIsNone(_factorize(np.zeros((0, 0))))
+
+
+    def test_factorize_accepts_regular(self):
+        """a regular matrix factorizes and solves correctly"""
+
+        from scipy.linalg import lu_solve
+
+        A = np.array([[-50.0, 0.0], [0.0, -2.0]])
+        lu = _factorize(A)
+
+        self.assertIsNotNone(lu)
+        self.assertTrue(np.allclose(lu_solve(lu, np.array([1.0, 4.0])), [-0.02, -2.0]))
+
+
+    def test_singular_step_is_finite(self):
+        """the correction stays finite where 'lu_solve' would give inf/nan"""
+
+        A = np.array([[-50.0, 0.0], [0.0, 0.0]])
+        dx = _singular_step(A, np.array([1.0, 0.382]))
+
+        self.assertTrue(np.all(np.isfinite(dx)))
+
+        #newton on the range, no correction on the null space
+        self.assertAlmostEqual(dx[0], 1.0 / -50.0, places=12)
+        self.assertAlmostEqual(dx[1], 0.0, places=12)
+
+
+    def test_newton_step_survives_singular_jacobian(self):
+        """the vector path does not produce inf/nan"""
+
+        opt = NewtonAnderson()
+        x = np.array([0.0, 0.0])
+        g = np.array([1.0, 0.382])
+        jac = np.eye(2) + np.array([[-50.0, 0.0], [0.0, 0.0]])
+
+        y, res = opt._newton(x, g, jac)
+
+        self.assertTrue(np.all(np.isfinite(y)), f"got {y}")
+        self.assertTrue(np.isfinite(res))
+
+
+    def test_newton_step_survives_zero_scalar_jacobian(self):
+        """the scalar path does not divide by zero
+
+        'jac - 1 == 0' is the scalar face of the same singularity, and an
+        'Integrator' ('dx/dt = u') hits it exactly.
+        """
+
+        opt = NewtonAnderson()
+        y, res = opt._newton(np.array([0.0]), np.array([1.0]), np.array([[1.0]]))
+
+        self.assertTrue(np.all(np.isfinite(y)), f"got {y}")
+        self.assertTrue(np.isfinite(res))
+
 
 # RUN TESTS LOCALLY ====================================================================
 

diff --git a/tests/pathsim/test_simulation.py b/tests/pathsim/test_simulation.py
@@ -1183,6 +1183,100 @@ def test_assemble_graph_resolves_all_connections(self):
         self.assertGreaterEqual(len(B.outputs), 1)
 
 
+class TestSteadyStateSingularJacobian(unittest.TestCase):
+    """
+    Regression test for a crash in 'steadystate' on systems whose blocks have a
+    singular Newton matrix.
+
+    A 'PID' has 'dx/dt = u' for its integrator state and a 'DynamicalSystem'
+    with 'dx/dt = -x**2 + u' has a zero Jacobian at 'x = 0'. Both make the
+    Newton matrix 'jac - I' singular, which used to produce inf/nan silently and
+    surface many iterations later, in an unrelated place, as
+    'LinAlgError: SVD did not converge in Linear Least Squares'.
+    """
+
+    def _closed_loop(self):
+        from pathsim.blocks import DynamicalSystem
+        from pathsim.blocks.ctrl import PID
+
+        ref, err = Constant(1.0), Adder("+-")
+        ctrl = PID(Kp=1.0, Ki=2.0, Kd=0.0, f_max=50)
+        plant = DynamicalSystem(
+            func_dyn=lambda x, u, t: -x**2 + u,
+            func_alg=lambda x, u, t: x,
+            initial_value=1.0,
+            jac_dyn=lambda x, u, t: -2*x
+            )
+        Sim = Simulation(
+            blocks=[ref, err, ctrl, plant, Scope()],
+            connections=[
+                Connection(ref, err[0]),
+                Connection(plant, err[1]),
+                Connection(err, ctrl),
+                Connection(ctrl, plant),
+                Connection(plant, Scope())
+                ],
+            dt=0.01, log=False
+            )
+        return Sim, ctrl, plant
+
+
+    def test_no_linalg_error(self):
+        """the failure mode is a truthful convergence report, not an SVD crash
+
+        This system is not solvable by the block-local steady-state scheme, so
+        it is expected to report non-convergence. What it must not do is crash
+        somewhere unrelated.
+        """
+        Sim, _, _ = self._closed_loop()
+
+        try:
+            Sim.steadystate()
+        except np.linalg.LinAlgError as exc:
+            self.fail(f"singular Newton matrix crashed the solve: {exc}")
+        except RuntimeError:
+            pass
+
+
+    def test_states_stay_finite(self):
+        """no inf/nan leaks into the block states"""
+        Sim, ctrl, plant = self._closed_loop()
+
+        try:
+            Sim.steadystate()
+        except RuntimeError:
+            pass
+
+        self.assertTrue(np.all(np.isfinite(np.atleast_1d(ctrl.state))))
+        self.assertTrue(np.all(np.isfinite(np.atleast_1d(plant.state))))
+
+
+    def test_integrator_loop_still_converges(self):
+        """a solvable system with the same singularity is unaffected
+
+        An 'Integrator' has a zero Jacobian too, so it takes the same code path;
+        closing a stable loop around it must still reach its equilibrium.
+        """
+        src, add, integ, gain = (
+            Constant(2.0), Adder("+-"), Integrator(0.0), Amplifier(1.0)
+            )
+        Sim = Simulation(
+            blocks=[src, add, integ, gain],
+            connections=[
+                Connection(src, add[0]),
+                Connection(integ, gain),
+                Connection(gain, add[1]),
+                Connection(add, integ)
+                ],
+            dt=0.01, log=False
+            )
+
+        Sim.steadystate()
+
+        #dx/dt = u - x -> x = u = 2
+        self.assertAlmostEqual(float(np.atleast_1d(integ.state)[0]), 2.0, places=4)
+
+
 # RUN TESTS LOCALLY ====================================================================
 
 if __name__ == '__main__':