Damp the Newton step in the solver when the full step makes things worse. Fixes #1247

src/system.cpp

@@ -316,32 +316,62 @@ bool System::NewtonSolve() {
     for(i = 0; i < mat.m; i++) {
         mat.B.num[i] = (mat.B.sym[i])->Eval();
     }
+
+    // Where we were before we took the step, so that we can take a shorter
+    // one from the same operating point if the first one is no good.
+    std::vector<double> prevVal(mat.n);
+
+    // Take the Newton step,
+    //      J(x_n) (x_{n+1} - x_n) = 0 - F(x_n)
+    // scaled by relax, and re-evaluate the functions, since the params have
+    // just changed. Returns false if that took us somewhere clearly useless.
+    auto takeStep = [&](double relax) {
+        for(int k = 0; k < mat.n; k++) {
+            Param *p = param.FindById(mat.param[k]);
+            p->val = prevVal[k] - relax*mat.X[k];
+            if(IsReasonable(p->val)) return false;
+        }
+        for(int k = 0; k < mat.m; k++) {
+            mat.B.num[k] = (mat.B.sym[k])->Eval();
+            if(IsReasonable(mat.B.num[k])) return false;
+        }
+        return true;
+    };
+
     do {
         // And evaluate the Jacobian at our initial operating point.
         EvalJacobian();
 
         if(!SolveLeastSquares()) break;
 
-        // Take the Newton step;
-        //      J(x_n) (x_{n+1} - x_n) = 0 - F(x_n)
         for(i = 0; i < mat.n; i++) {
-            Param *p = param.FindById(mat.param[i]);
-            p->val -= mat.X[i];
-            if(IsReasonable(p->val)) {
-                // Very bad, and clearly not convergent
-                return false;
-            }
+            prevVal[i] = param.FindById(mat.param[i])->val;
         }
-
-        // Re-evalute the functions, since the params have just changed.
-        for(i = 0; i < mat.m; i++) {
-            mat.B.num[i] = (mat.B.sym[i])->Eval();
-            if(IsReasonable(mat.B.num[i])) {
+        const double err = mat.B.num.squaredNorm();
+
+        // A Newton step is only as good as the linearization of F about our
+        // operating point, so where F is strongly curved the full step can
+        // land much further from the solution than it started; and if it
+        // lands near a critical point of F, then the step after that one is
+        // enormous, and the geometry runs away to nowhere. But the step is a
+        // descent direction for |F|^2, so if the full step makes things
+        // worse, then a short enough step along the same direction makes
+        // them better. So backtrack until the residual actually decreases.
+        bool accepted = false;
+        double relax = 1.0;
+        for(int tries = 0; tries < 8 && !accepted; tries++, relax /= 2) {
+            accepted = takeStep(relax) && mat.B.num.squaredNorm() < err;
+        }
+        if(!accepted) {
+            // No step along this direction is an improvement, so just take
+            // the whole thing, like we always used to; if that was a bad
+            // idea, then the iteration limit below will catch it.
+            if(!takeStep(1.0)) {
                 // Very bad, and clearly not convergent
                 return false;
             }
         }
-        
+
         // Check for convergence
         converged = true;
         for(i = 0; i < mat.m; i++) {

test/CMakeLists.txt

@@ -19,6 +19,7 @@ set(testsuite_SOURCES
     core/expr/test.cpp
     core/locale/test.cpp
     core/path/test.cpp
+    core/solver/test.cpp
     constraint/points_coincident/test.cpp
     constraint/pt_pt_distance/test.cpp
     constraint/pt_plane_distance/test.cpp