Fix bug in Trust Region subproblem implementation

Changed the ITrustRegionSubproblem interface to directly accept scaled gradient and Hessian matrices instead of accessing them from IObjectiveModel.

Author: diluculo

src/Numerics/Optimization/TrustRegion/ITrustRegionSubProblem.cs

@@ -20,10 +20,11 @@ public interface ITrustRegionSubproblem
         bool HitBoundary { get; }
 
         /// <summary>
-        /// Solves the trust region subproblem for the specified objective model and trust region radius.
+        /// Solves the trust region subproblem using the provided gradient and Hessian.
         /// </summary>
-        /// <param name="objective">The objective model containing function, gradient, and Hessian information.</param>
-        /// <param name="radius">The radius of the trust region.</param>
-        void Solve(IObjectiveModel objective, double radius);
+        /// <param name="gradient">The scaled gradient vector</param>
+        /// <param name="hessian">The scaled Hessian matrix</param>
+        /// <param name="delta">Trust region radius</param>
+        void Solve(Vector<double> gradient, Matrix<double> hessian, double delta);
     }
 }

src/Numerics/Optimization/TrustRegion/Subproblems/DogLegSubproblem.cs

@@ -8,17 +8,14 @@ internal class DogLegSubproblem : ITrustRegionSubproblem
 
         public bool HitBoundary { get; private set; }
 
-        public void Solve(IObjectiveModel objective, double delta)
+        public void Solve(Vector<double> gradient, Matrix<double> hessian, double delta)
         {
-            var Gradient = objective.Gradient;
-            var Hessian = objective.Hessian;
-
             // newton point, the Gauss–Newton step by solving the normal equations
-            var Pgn = -Hessian.PseudoInverse() * Gradient; // Hessian.Solve(Gradient) fails so many times...
+            var Pgn = -hessian.PseudoInverse() * gradient; // hessian.Solve(gradient) fails so many times...
 
             // cauchy point, steepest descent direction is given by
-            var alpha = Gradient.DotProduct(Gradient) / (Hessian * Gradient).DotProduct(Gradient);
-            var Psd = -alpha * Gradient;
+            var alpha = gradient.DotProduct(gradient) / (hessian * gradient).DotProduct(gradient);
+            var Psd = -alpha * gradient;
 
             // update step and prectted reduction
             if (Pgn.L2Norm() <= delta)

src/Numerics/Optimization/TrustRegion/Subproblems/NewtonCGSubproblem.cs

@@ -9,29 +9,36 @@ internal class NewtonCGSubproblem : ITrustRegionSubproblem
 
         public bool HitBoundary { get; private set; }
 
-        public void Solve(IObjectiveModel objective, double delta)
+        public void Solve(Vector<double> gradient, Matrix<double> hessian, double delta)
         {
-            var Gradient = objective.Gradient;
-            var Hessian = objective.Hessian;
-
             // define tolerance
-            var gnorm = Gradient.L2Norm();
+            var gnorm = gradient.L2Norm();
             var tolerance = Math.Min(0.5, Math.Sqrt(gnorm)) * gnorm;
 
             // initialize internal variables
-            var z = Vector<double>.Build.Dense(Hessian.RowCount);
-            var r = Gradient;
+            var z = Vector<double>.Build.Dense(hessian.RowCount);
+            var r = gradient;
             var d = -r;
 
             while (true)
             {
-                var Bd = Hessian * d;
+                var Bd = hessian * d;
                 var dBd = d.DotProduct(Bd);
 
                 if (dBd <= 0)
                 {
+                    // Direction of negative curvature found
+                    // Calculate two boundary points and choose the one with lower model value
                     var t = Util.FindBeta(1, z, d, delta);
-                    Pstep = z + t.Item1 * d;
+                    var pa = z + t.Item1 * d;
+                    var pb = z + t.Item2 * d;
+
+                    // Evaluate quadratic model at both points
+                    var valueA = Util.CalculateQuadraticModel(gradient, hessian, pa);
+                    var valueB = Util.CalculateQuadraticModel(gradient, hessian, pb);
+
+                    // Choose the point with the lower model value
+                    Pstep = valueA < valueB ? pa : pb;
                     HitBoundary = true;
                     return;
                 }

src/Numerics/Optimization/TrustRegion/Subproblems/Util.cs

@@ -5,6 +5,14 @@ namespace MathNet.Numerics.Optimization.TrustRegion.Subproblems
 {
     internal static class Util
     {
+        /// <summary>
+        /// Finds the two intersection points of a line with the trust region boundary.
+        /// </summary>
+        /// <param name="alpha">Scaling factor for steepest descent direction</param>
+        /// <param name="sd">Steepest descent direction vector</param>
+        /// <param name="gn">Gauss-Newton direction vector</param>
+        /// <param name="delta">Trust region radius</param>
+        /// <returns>A tuple containing two beta values, sorted from low to high</returns>
         public static (double, double) FindBeta(double alpha, Vector<double> sd, Vector<double> gn, double delta)
         {
             // Pstep is intersection of the trust region boundary
@@ -29,5 +37,24 @@ public static (double, double) FindBeta(double alpha, Vector<double> sd, Vector<
             // return sorted beta
             return beta1 < beta2 ? (beta1, beta2) : (beta2, beta1);
         }
+
+        /// <summary>
+        /// Calculates the value of the quadratic model at a given point.
+        /// The quadratic model is defined as:
+        ///     m(p) = g^T * p + 0.5 * p^T * H * p
+        /// where g is the gradient and H is the Hessian at the current point.
+        /// </summary>
+        /// <param name="gradient">The gradient vector</param>
+        /// <param name="hessian">The Hessian matrix</param>
+        /// <param name="p">The point at which to evaluate the quadratic model</param>
+        /// <returns>The value of the quadratic model at point p</returns>
+        public static double CalculateQuadraticModel(Vector<double> gradient, Matrix<double> hessian, Vector<double> p)
+        {
+            // Quadratic model: m(p) = g^T * p + 0.5 * p^T * H * p
+            var linearTerm = gradient.DotProduct(p);
+            var quadraticTerm = 0.5 * p.DotProduct(hessian * p);
+
+            return linearTerm + quadraticTerm;
+        }
     }
 }

src/Numerics/Optimization/TrustRegion/TrustRegionMinimizerBase.cs

@@ -88,8 +88,8 @@ public NonlinearMinimizationResult Minimum(ITrustRegionSubproblem subproblem, IO
             //    Residuals, R = L(y - f(x; p))
             //    Residual sum of squares, RSS = ||R||^2 = R.DotProduct(R)
             //    Jacobian J = df(x; p)/dp
-            //    Gradient g = -J'W(y − f(x; p)) = -J'LR
-            //    Approximated Hessian H = J'WJ
+            //    gradient g = -J'W(y − f(x; p)) = -J'LR
+            //    Approximated hessian H = J'WJ
             //
             // The trust region algorithm is summarized as follows:
             //    initially set trust-region radius, Δ
@@ -120,7 +120,7 @@ public NonlinearMinimizationResult Minimum(ITrustRegionSubproblem subproblem, IO
 
             objective.SetParameters(initialGuess, isFixed);
 
-            ExitCondition exitCondition = ExitCondition.None;
+            var exitCondition = ExitCondition.None;
 
             // First, calculate function values and setup variables
             var P = ProjectToInternalParameters(initialGuess); // current internal parameters
@@ -151,11 +151,11 @@ public NonlinearMinimizationResult Minimum(ITrustRegionSubproblem subproblem, IO
                 exitCondition = ExitCondition.Converged; // SmallRSS
             }
 
-            // evaluate projected gradient and Hessian
-            var (Gradient, Hessian) = EvaluateJacobian(objective, P);
+            // evaluate projected gradient and hessian
+            var (gradient, hessian) = EvaluateJacobian(objective, P);
 
             // if ||g||_oo <= gtol, found and stop
-            if (Gradient.InfinityNorm() <= gradientTolerance)
+            if (gradient.InfinityNorm() <= gradientTolerance)
             {
                 exitCondition = ExitCondition.RelativeGradient; // SmallGradient
             }
@@ -166,22 +166,22 @@ public NonlinearMinimizationResult Minimum(ITrustRegionSubproblem subproblem, IO
             }
 
             // initialize trust-region radius, Δ
-            double delta = Gradient.DotProduct(Gradient) / (Hessian * Gradient).DotProduct(Gradient);
+            var delta = gradient.DotProduct(gradient) / (hessian * gradient).DotProduct(gradient);
             delta = Math.Max(1.0, Math.Min(delta, maxDelta));
 
-            int iterations = 0;
+            var iterations = 0;
             bool hitBoundary;
             while (iterations < maximumIterations && exitCondition == ExitCondition.None)
             {
                 iterations++;
 
                 // solve the subproblem
-                subproblem.Solve(objective, delta);
+                subproblem.Solve(gradient, hessian, delta);
                 Pstep = subproblem.Pstep;
                 hitBoundary = subproblem.HitBoundary;
 
                 // predicted reduction = L(0) - L(Δp) = -Δp'g - 1/2 * Δp'HΔp
-                var predictedReduction = -Gradient.DotProduct(Pstep) - 0.5 * Pstep.DotProduct(Hessian * Pstep);
+                var predictedReduction = -gradient.DotProduct(Pstep) - 0.5 * Pstep.DotProduct(hessian * Pstep);
 
                 if (Pstep.L2Norm() <= stepTolerance * (stepTolerance + P.L2Norm()))
                 {
@@ -201,7 +201,7 @@ public NonlinearMinimizationResult Minimum(ITrustRegionSubproblem subproblem, IO
                 }
 
                 // calculate the ratio of the actual to the predicted reduction.
-                double rho = (predictedReduction != 0)
+                var rho = (predictedReduction != 0)
                         ? (RSS - RSSnew) / predictedReduction
                         : 0.0;
 
@@ -225,11 +225,11 @@ public NonlinearMinimizationResult Minimum(ITrustRegionSubproblem subproblem, IO
                     Pnew.CopyTo(P);
                     RSS = RSSnew;
 
-                    // evaluate projected gradient and Hessian
-                    (Gradient, Hessian) = EvaluateJacobian(objective, P);
+                    // evaluate projected gradient and hessian
+                    (gradient, hessian) = EvaluateJacobian(objective, P);
 
                     // if ||g||_oo <= gtol, found and stop
-                    if (Gradient.InfinityNorm() <= gradientTolerance)
+                    if (gradient.InfinityNorm() <= gradientTolerance)
                     {
                         exitCondition = ExitCondition.RelativeGradient;
                     }