refactor: Enhance Newton's method with early termination, improved line search, and refined candidate processing in extremum evaluations

Author: dpasukhi

src/FoundationClasses/TKMath/MathSys/MathSys_Newton2D.hxx

@@ -450,14 +450,38 @@ Newton2DResult Newton2DSymmetric(const Function& theFunc,
     aNewU = std::max(aExtUMin, std::min(aExtUMax, aNewU));
     aNewV = std::max(aExtVMin, std::min(aExtVMax, aNewV));
 
+    // Step size relative to domain - used for early termination and line search skip
+    const double aStepNormSq = aDU * aDU + aDV * aDV;
+    const double aDomainSizeSq = aDomainU * aDomainU + aDomainV * aDomainV;
+    const double aRelStepSq = aStepNormSq / std::max(aDomainSizeSq, THE_HESSIAN_DEGENERACY_ABS);
+
+    // Early convergence: if step is very small relative to domain, we're converged
+    if (aRelStepSq < THE_NEWTON_STEP_TOL_FACTOR * THE_NEWTON_STEP_TOL_FACTOR)
+    {
+      aRes.FNorm  = std::sqrt(aFNormSq);
+      aRes.Status = Status::OK;
+      return aRes;
+    }
+
+    // Skip line search for small steps - Newton is working well
+    // This avoids expensive Value() calls when not needed
+    const bool aSkipLineSearch = (aRelStepSq < THE_NEWTON2D_SKIP_LINESEARCH_SQ);
+
     // Line search if step increases function value significantly
     double aNewF1, aNewF2;
-    if (theFunc.Value(aNewU, aNewV, aNewF1, aNewF2))
+    if (!aSkipLineSearch && theFunc.Value(aNewU, aNewV, aNewF1, aNewF2))
     {
       double aNewFNormSq = aNewF1 * aNewF1 + aNewF2 * aNewF2;
 
-      // Backtrack if function increased too much
-      if (aNewFNormSq > aFNormSq * THE_NEWTON2D_BACKTRACK_TRIGGER)
+      // Compute directional derivative for Armijo condition: grad_f . d
+      // For f = 0.5*(F1^2 + F2^2), grad_f = [F1, F2] (simplified)
+      // d = [dU, dV], so grad_f . d = F1*dU + F2*dV
+      const double aDirectionalDeriv = aF1 * aDU + aF2 * aDV;
+
+      // Backtrack if function increased too much OR Armijo not satisfied
+      // Armijo: f(x + d) <= f(x) + c1 * grad_f . d
+      const bool aArmijoFailed = (aNewFNormSq > aFNormSq + THE_ARMIJO_C1 * aDirectionalDeriv);
+      if (aNewFNormSq > aFNormSq * THE_NEWTON2D_BACKTRACK_TRIGGER || aArmijoFailed)
       {
         // Use quadratic interpolation for first estimate
         // Model: f(α) ≈ f(0) + f'(0)*α + c*α², where c = (f(1) - f(0) - f'(0))/1
@@ -492,7 +516,9 @@ Newton2DResult Newton2DSymmetric(const Function& theFunc,
           if (theFunc.Value(aTryU, aTryV, aTryF1, aTryF2))
           {
             double aTryFNormSq = aTryF1 * aTryF1 + aTryF2 * aTryF2;
-            if (aTryFNormSq < aFNormSq * THE_NEWTON2D_BACKTRACK_ACCEPT)
+            // Accept if Armijo condition satisfied: f(x + α*d) <= f(x) + c1*α*(grad·d)
+            const double aArmijoThreshold = aFNormSq + THE_ARMIJO_C1 * aAlpha * aDirectionalDeriv;
+            if (aTryFNormSq <= aArmijoThreshold || aTryFNormSq < aFNormSq * THE_NEWTON2D_BACKTRACK_ACCEPT)
             {
               aNewU = aTryU;
               aNewV = aTryV;

src/FoundationClasses/TKMath/MathUtils/MathUtils_Config.hxx

@@ -47,7 +47,8 @@ constexpr double THE_NEWTON2D_MAX_STEP_RATIO = 0.5;
 
 //! Domain extension factor for soft boundary handling.
 //! Solutions slightly outside domain (by this fraction) are allowed if converged.
-constexpr double THE_NEWTON2D_DOMAIN_EXT = 0.01;
+//! Reduced to 1e-4 to closely match old algorithm behavior.
+constexpr double THE_NEWTON2D_DOMAIN_EXT = 1.0e-4;
 
 //! Stagnation progress ratio - improvement required each iteration to avoid stagnation.
 //! Value of 0.999 means at least 0.1% improvement required.
@@ -57,7 +58,11 @@ constexpr double THE_NEWTON2D_STAGNATION_RATIO = 0.999;
 constexpr int THE_NEWTON2D_STAGNATION_COUNT = 3;
 
 //! Maximum line search backtracking iterations.
-constexpr int THE_NEWTON2D_LINE_SEARCH_MAX = 12;
+constexpr int THE_NEWTON2D_LINE_SEARCH_MAX = 8;
+
+//! Relative step threshold below which line search is skipped.
+//! If step^2 / domain^2 < threshold, Newton is converging well.
+constexpr double THE_NEWTON2D_SKIP_LINESEARCH_SQ = 0.01;
 
 //! Tolerance relaxation factor for accepting stagnated solutions.
 //! Stagnated solution accepted if |F| < tolerance * factor.
@@ -74,8 +79,27 @@ constexpr double THE_NEWTON2D_BACKTRACK_ACCEPT = 1.2;
 //! Maximum relaxation factor for solutions at max iterations.
 constexpr double THE_NEWTON2D_MAXITER_RELAX = 10.0;
 
-//! High precision tolerance target (5e-8, better than Precision::Confusion).
-constexpr double THE_HIGH_PRECISION_TOL = 5.0e-8;
+//==================================================================================================
+//! @name Hessian Classification Constants
+//! Constants for extremum classification using second derivatives.
+//==================================================================================================
+
+//! Relative tolerance for Hessian degeneracy detection.
+//! If |det(H)| < threshold * |H_ii| * |H_jj|, Hessian is considered degenerate.
+constexpr double THE_HESSIAN_DEGENERACY_REL = 1.0e-8;
+
+//! Absolute tolerance for Hessian degeneracy detection.
+//! Minimum threshold for determinant comparison.
+constexpr double THE_HESSIAN_DEGENERACY_ABS = 1.0e-20;
+
+//==================================================================================================
+//! @name Line Search Constants
+//! Constants for line search acceptance conditions.
+//==================================================================================================
+
+//! Armijo condition constant (sufficient decrease).
+//! Step accepted if: f(x + alpha*d) <= f(x) + c1 * alpha * grad_f . d
+constexpr double THE_ARMIJO_C1 = 1.0e-4;
 
 //! Configuration for iterative solvers.
 //! Provides common settings for convergence criteria and iteration limits.

src/ModelingData/TKGeomBase/ExtremaPS/ExtremaPS.hxx

@@ -75,11 +75,27 @@ constexpr double THE_CELL_EXPAND_RATIO = 0.1;
 //! Ratio of distance for scaling tolerance in cell search.
 constexpr double THE_DISTANCE_SCALE_RATIO = 0.1;
 
+//==================================================================================================
+//! @name Early Termination and Candidate Processing Constants
+//! Constants controlling candidate filtering, early termination, and iteration limits.
+//==================================================================================================
+
 //! Threshold for skipping max candidates that are clearly worse.
+//! For Max mode: skip if estDist < bestDist * threshold (i.e., 90% of best).
 constexpr double THE_MAX_SKIP_THRESHOLD = 0.9;
 
-//! Threshold for refined distance comparison (slightly less than 1.0).
-constexpr double THE_REFINED_DIST_THRESHOLD = 0.99;
+//! Threshold margin for Min mode early termination.
+//! For Min mode: skip if estDist > bestDist * (1 + margin).
+//! Value of 0.2 means skip candidates 20% worse than current best.
+constexpr double THE_MIN_SKIP_MARGIN = 0.2;
+
+//! Maximum number of candidates to process before stopping.
+//! Limits Newton iterations per grid scan for performance.
+constexpr int THE_MAX_CANDIDATES_TO_PROCESS = 50;
+
+//! Threshold for refined distance comparison.
+//! Grid fallback accepted if refinedDist < currentMin * threshold.
+constexpr double THE_REFINED_DIST_THRESHOLD = 1.0;
 
 //! Relaxation factor for gradient-based zero detection.
 //! Multiplied with tolerance for more robust extremum detection.
@@ -89,24 +105,46 @@ constexpr double THE_GRADIENT_TOL_FACTOR = 100.0;
 //! Multiplier for Newton retry tolerance when initial Newton fails.
 constexpr double THE_NEWTON_RETRY_TOL_FACTOR = 10.0;
 
-//! High precision tolerance target (5e-8, better than Precision::Confusion ~1e-7).
-//! Used for multi-start Newton and two-phase refinement.
-constexpr double THE_HIGH_PRECISION_TOL = 5.0e-8;
+//==================================================================================================
+//! @name Cache and Coherent Scanning Constants
+//! Constants for spatial coherence optimization and trajectory prediction.
+//==================================================================================================
 
 //! Squared threshold for spatial coherence optimization.
 //! Query points within this squared distance use cached solution.
+//! Value of 100.0 corresponds to distance of 10.0 units.
 constexpr double THE_COHERENCE_THRESHOLD_SQ = 100.0;
 
 //! Minimum ratio for trajectory prediction (step magnitude ratio).
+//! Ratios below this indicate irregular step sizes.
 constexpr double THE_TRAJECTORY_MIN_RATIO = 0.5;
 
 //! Maximum ratio for trajectory prediction (step magnitude ratio).
+//! Ratios above this indicate acceleration/deceleration too large.
 constexpr double THE_TRAJECTORY_MAX_RATIO = 2.0;
 
 //! Minimum cosine for trajectory prediction (direction alignment).
 //! Values above this indicate roughly same direction (< ~45 degrees).
 constexpr double THE_TRAJECTORY_MIN_COS = 0.7;
 
+//! Minimum cosine for quadratic extrapolation (high coherence).
+//! Values above this indicate very smooth trajectory (< ~25 degrees).
+constexpr double THE_TRAJECTORY_QUADRATIC_COS = 0.9;
+
+//! Cache size for coherent scanning (number of solutions stored).
+//! Increased from 3 to 5 to enable quadratic extrapolation.
+constexpr int THE_CACHE_SIZE = 5;
+
+//! Sanity factor for cache result validation.
+//! If cache-based Newton result has squared distance > (theTol^2 * factor),
+//! it's rejected as likely converging to wrong local minimum.
+//! Value of 1e6 allows reasonable tolerance range while catching gross errors.
+constexpr double THE_CACHE_SANITY_FACTOR = 1.0e6;
+
+//! Velocity adjustment limit for trajectory prediction.
+//! Limits acceleration factor to prevent overshoot.
+constexpr double THE_VELOCITY_ADJUST_LIMIT = 1.5;
+
 //! Maximum number of Newton iterations for grid optimization.
 constexpr int THE_MAX_GOLDEN_ITERATIONS = 50;