Candidate fix for enumpoly numerical stability

gambitproject · Jan 10, 2025 · 247ab0e · 247ab0e
1 parent c1a222c
commit 247ab0e
Showing 1 changed file with 7 additions and 3 deletions.
diff --git a/src/solvers/enumpoly/polysolver.cc b/src/solvers/enumpoly/polysolver.cc
@@ -76,7 +76,7 @@ bool PolynomialSystemSolver::SystemHasNoRootsIn(const Rectangle<double> &r,
 
 namespace {
 
-static bool fuzzy_equals(double x, double y)
+bool fuzzy_equals(double x, double y)
 {
   const double epsilon = 0.000000001;
 
@@ -220,7 +220,12 @@ Vector<double> PolynomialSystemSolver::NewtonPolishOnce(const Vector<double> &po
 {
   Vector<double> oldevals = m_derivatives.ValuesOfRootPolys(point, NumEquations());
   Matrix<double> Df = m_derivatives.DerivativeMatrix(point, NumEquations());
-  return point - (Df.Transpose() * SquareMatrix<double>(Df * Df.Transpose()).Inverse()) * oldevals;
+  auto transpose = Df.Transpose();
+  auto sqmat = SquareMatrix<double>(Df * transpose);
+  if (std::abs(sqmat.Determinant()) <= 1.0e-9) {
+    throw SingularMatrixException();
+  }
+  return point - (transpose * sqmat.Inverse()) * oldevals;
 }
 
 Vector<double> PolynomialSystemSolver::SlowNewtonPolishOnce(const Vector<double> &point) const
@@ -297,7 +302,6 @@ void PolynomialSystemSolver::FindRoots(List<Vector<double>> &rootlist, const Rec
     }
     return;
   }
-
   for (const auto &cell : r.Orthants()) {
     if (max_roots == 0 || rootlist.size() < max_roots) {
       if (iterations >= max_iterations || depth == MAX_DEPTH) {