diff --git a/burnman/optimize/nonlinear_solvers.py b/burnman/optimize/nonlinear_solvers.py index fb7e36ed..11cb7000 100644 --- a/burnman/optimize/nonlinear_solvers.py +++ b/burnman/optimize/nonlinear_solvers.py @@ -3,7 +3,6 @@ import numpy.typing as npt from typing import Any from scipy.linalg import lu_factor, lu_solve -from collections import namedtuple from types import SimpleNamespace @@ -200,6 +199,7 @@ def __init__( self.cond_lu_thresh = cond_lu_thresh self.cond_lstsq_thresh = cond_lstsq_thresh self.eps = 2.0 * np.finfo(float).eps + self.max_condition_number = 1.0 / np.finfo(float).eps assert np.all( self._constraints(self.guess) < self.eps @@ -222,8 +222,12 @@ def _update_lmda( h: npt.NDArray[np.float64], lmda_bounds: tuple[float, float], ) -> float: - assert lmda_bounds[1] < 1.0 + self.eps - assert lmda_bounds[0] > 1.0e-8 - self.eps + assert ( + lmda_bounds[1] < 1.0 + self.eps + ), f"max_lambda must be <= 1.0, got {lmda_bounds[1]}" + assert ( + lmda_bounds[0] > 1.0e-8 - self.eps + ), f"min_lambda must be > 1.0e-8, got {lmda_bounds[0]}" lmda_j = min(1.0 / (h + self.eps), lmda_bounds[1]) return max(lmda_j, lmda_bounds[0]) @@ -598,8 +602,10 @@ def _termination_info( n_it: int, max_iterations: int, violated_constraints: list[int], + condition_number: float, ) -> tuple[int, str, bool]: - if converged: + singular_jacobian = condition_number > self.max_condition_number + if converged and not singular_jacobian: return ( True, 0, @@ -611,7 +617,7 @@ def _termination_info( 1, f"The function is too non-linear for lower lambda bound ({lmda_bounds[0]})", ) - if persistent_bound_violation: + if persistent_bound_violation and not singular_jacobian: return ( False, 2, @@ -622,6 +628,30 @@ def _termination_info( ) if n_it == max_iterations: return False, 3, f"The solver reached max_iterations ({max_iterations})" + if singular_jacobian and not converged: + return ( + False, + 4, + ( + f"The system was effectively singular (cond={condition_number:.2e}) " + f"and tried to leave the feasible region at iteration {n_it} by " + "crossing the constraints with the following indices: " + f"{[i for i, lmda in violated_constraints]}.\n" + "You may want to try a different initial guess " + "or check the Jacobian implementation. " + "It may be that the problem is ill-posed." + ), + ) + if singular_jacobian and converged: + return ( + True, + 5, + ( + f"The solver successfully found a root after {n_it} iterations," + f"but the system is effectively singular (cond={condition_number:.2e}). " + "You should check the problem formulation before trusting this result." + ), + ) raise Exception("Unknown termination of solver") def solve(self) -> Solution: @@ -645,10 +675,12 @@ def solve(self) -> Solution: * 1 -- Failure: lambda reached its minimum bound * 2 -- Failure: descent direction violates constraints * 3 -- Failure: maximum iterations reached + * 4 -- Failure: singular system and violated constraints + * 5 -- Successful convergence but singular system * **text** (str) -- Human-readable termination description. * **n_it** (int) -- Number of Newton iterations performed. - * **iterates** (namedtuple, optional) -- Present if + * **iterates** (Iterates instance, optional) -- Present if ``store_iterates=True``. Contains iteration history: * **x** (list[np.ndarray]) -- Solution for each iteration. @@ -676,10 +708,16 @@ def solve(self) -> Solution: and not converged ): sol.J = self.J(sol.x) + condition_number = np.linalg.cond(sol.J) luJ = lu_factor(sol.J) dx = lu_solve(luJ, -sol.F) dx_norm = np.linalg.norm(dx, ord=2) + is_singular = condition_number > self.max_condition_number + if is_singular and any(np.isnan(dx)): + lmda_bounds = [0.0, 0.0] + break + lmda_bounds = self.lambda_bounds(dx, sol.x) h = ( lmda @@ -710,7 +748,6 @@ def solve(self) -> Solution: F_j = self.F(x_j) # Regularise ill-conditioned Jacobian - condition_number = np.linalg.cond(sol.J) if condition_number < self.cond_lu_thresh: dxbar_j = lu_solve(luJ, -F_j) else: @@ -751,14 +788,16 @@ def solve(self) -> Solution: sol.iterates.append(sol.x, sol.F, lmda) # Final adjustment for constraints - if converged and not persistent_bound_violation: - sol.x = x_j + dxbar_j - if not np.all(self._constraints(sol.x) <= 0.0): - sol.x -= dxbar_j - - sol.F = self.F(sol.x) - sol.F_norm = np.linalg.norm(sol.F, ord=2) - sol.J = self.J(sol.x) + if condition_number < self.max_condition_number: + if converged and not persistent_bound_violation: + sol.x = x_j + dxbar_j + if not np.all(self._constraints(sol.x) <= 0.0): + sol.x -= dxbar_j + + sol.F = self.F(sol.x) + sol.F_norm = np.linalg.norm(sol.F, ord=2) + sol.J = self.J(sol.x) + condition_number = np.linalg.cond(sol.J) sol.success, sol.code, sol.text = self._termination_info( converged, @@ -768,6 +807,7 @@ def solve(self) -> Solution: sol.n_it, self.max_iterations, locals().get("violated_constraints", []), + condition_number, ) return sol diff --git a/burnman/tools/equilibration.py b/burnman/tools/equilibration.py index 802cf0fa..a89da869 100644 --- a/burnman/tools/equilibration.py +++ b/burnman/tools/equilibration.py @@ -497,7 +497,13 @@ def lambda_bounds(dx, x, endmembers_per_phase): for i, step in enumerate(np.abs(dx)) ] ) - + if np.isnan(max_lmda): + raise ValueError( + "NaN in lambda bounds calculation.\n" + f"dx: {dx}\n" + f"x: {x}\n" + f"endmembers_per_phase: {endmembers_per_phase}." + ) return (1.0e-8, max_lmda) diff --git a/tests/test_equilibration.py b/tests/test_equilibration.py index f068ef3c..321d1fd3 100644 --- a/tests/test_equilibration.py +++ b/tests/test_equilibration.py @@ -244,6 +244,46 @@ def test_ol_wad_eqm_compositional_constraint(self): [1620.532183457096, 0.45, 0.6791743], ) + def test_convergence_with_singular_system(self): + + composition = {"Mg": 1.0, "Fe": 1.0, "Si": 1.0, "O": 4.0} + a = burnman.Composite( + [SLB_2011.mg_fe_olivine(), SLB_2011.mg_fe_olivine()], [0.5, 0.5] + ) + + a.phases[0].set_composition([0.1, 0.9]) + a.phases[1].set_composition([0.9, 0.1]) + equality_constraints = [("P", 1.0e5), ("T", 1000.0)] + sol, _ = equilibrate(composition, a, equality_constraints) + self.assertTrue(sol.success) # Expect successful convergence + self.assertTrue(sol.code == 5) # Expect singular system at solution + + def test_ill_posed_problem(self): + + composition = {"Mg": 1.8, "Fe": 0.2, "Si": 1.0, "O": 4.0} + a = burnman.Composite( + [SLB_2011.mg_fe_olivine(), SLB_2011.mg_fe_olivine()], [0.5, 0.5] + ) + + a.phases[0].set_composition([0.4, 0.6]) + a.phases[1].set_composition([0.5, 0.5]) + equality_constraints = [("P", 1.0e5), ("T", 300.0)] + sol, _ = equilibrate(composition, a, equality_constraints) + self.assertTrue(sol.code == 2) # Expect problem to leave the feasible region + + def test_ill_conditioned_jacobian(self): + + composition = {"Mg": 1.8, "Fe": 0.2, "Si": 1.0, "O": 4.0} + a = burnman.Composite( + [SLB_2011.mg_fe_olivine(), SLB_2011.mg_fe_olivine()], [0.5, 0.5] + ) + + a.phases[0].set_composition([0.4, 0.6]) + a.phases[1].set_composition([0.4, 0.6]) + equality_constraints = [("P", 1.0e5), ("T", 300.0)] + sol, _ = equilibrate(composition, a, equality_constraints) + self.assertTrue(sol.code == 4) # Expect singular system + if __name__ == "__main__": unittest.main()