diff --git a/src/Numerics/Optimization/ITrustRegionSubProblem.cs b/src/Numerics/Optimization/ITrustRegionSubProblem.cs index 7e2d7e07..5305fc44 100644 --- a/src/Numerics/Optimization/ITrustRegionSubProblem.cs +++ b/src/Numerics/Optimization/ITrustRegionSubProblem.cs @@ -5,7 +5,6 @@ namespace MathNet.Numerics.Optimization public interface ITrustRegionSubproblem { Vector Pstep { get; } - double PredictedReduction { get; } bool HitBoundary { get; } void Solve(IObjectiveModel objective, double radius); diff --git a/src/Numerics/Optimization/Subproblems/QuadraticSubproblem.cs b/src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs similarity index 80% rename from src/Numerics/Optimization/Subproblems/QuadraticSubproblem.cs rename to src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs index 499982eb..e3f06d15 100644 --- a/src/Numerics/Optimization/Subproblems/QuadraticSubproblem.cs +++ b/src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs @@ -2,12 +2,10 @@ namespace MathNet.Numerics.Optimization.Subproblems { - internal class QuadraticSubproblem : ITrustRegionSubproblem + internal class DogLegSubproblem : ITrustRegionSubproblem { public Vector Pstep { get; private set; } - public double PredictedReduction { get; private set; } - public bool HitBoundary { get; private set; } public void Solve(IObjectiveModel objective, double delta) @@ -32,14 +30,12 @@ namespace MathNet.Numerics.Optimization.Subproblems // Pgn is inside trust region radius HitBoundary = false; Pstep = Pgn; - PredictedReduction = RSS; } else if (alpha * Psd.L2Norm() >= delta) { // Psd is outside trust region radius HitBoundary = true; Pstep = delta / Psd.L2Norm() * Psd; - PredictedReduction = delta * (2.0 * (alpha * Gradient).L2Norm() - delta) / 2.0 / alpha; } else { @@ -47,7 +43,6 @@ namespace MathNet.Numerics.Optimization.Subproblems HitBoundary = true; var beta = Util.FindBeta(alpha, Psd, Pgn, delta).Item2; Pstep = alpha * Psd + beta * (Pgn - alpha * Psd); - PredictedReduction = 0.5 * alpha * (1 - beta) * (1 - beta) * Gradient.DotProduct(Gradient) + beta * (2 - beta) * RSS; } } } diff --git a/src/Numerics/Optimization/Subproblems/Util.cs b/src/Numerics/Optimization/Subproblems/Util.cs index c15f636b..4ea986f8 100644 --- a/src/Numerics/Optimization/Subproblems/Util.cs +++ b/src/Numerics/Optimization/Subproblems/Util.cs @@ -26,7 +26,10 @@ namespace MathNet.Numerics.Optimization.Subproblems var beta1 = -aux / 2.0 / a; var beta2 = -2.0 * c / aux; - return new Tuple(beta1, beta2); + // return sorted beta + return (beta1 < beta2) + ? new Tuple(beta1, beta2) + : new Tuple(beta2, beta1); } } } diff --git a/src/Numerics/Optimization/TrustRegionDogLegMinimizer.cs b/src/Numerics/Optimization/TrustRegionDogLegMinimizer.cs index af4efc4d..596f7637 100644 --- a/src/Numerics/Optimization/TrustRegionDogLegMinimizer.cs +++ b/src/Numerics/Optimization/TrustRegionDogLegMinimizer.cs @@ -6,7 +6,7 @@ namespace MathNet.Numerics.Optimization public sealed class TrustRegionDogLegMinimizer : TrustRegionMinimizerBase { public TrustRegionDogLegMinimizer(double gradientTolerance = 1E-8, double stepTolerance = 1E-8, double functionTolerance = 1E-8, double radiusTolerance = 1E-8, int maximumIterations = -1) - : base(TrustRegionSubproblem.Quadratic(), gradientTolerance, stepTolerance, functionTolerance, radiusTolerance, maximumIterations) + : base(TrustRegionSubproblem.DogLeg(), gradientTolerance, stepTolerance, functionTolerance, radiusTolerance, maximumIterations) { } } } diff --git a/src/Numerics/Optimization/TrustRegionMinimizerBase.cs b/src/Numerics/Optimization/TrustRegionMinimizerBase.cs index 4bfc315d..f4125e8f 100644 --- a/src/Numerics/Optimization/TrustRegionMinimizerBase.cs +++ b/src/Numerics/Optimization/TrustRegionMinimizerBase.cs @@ -189,8 +189,9 @@ namespace MathNet.Numerics.Optimization // solve the subproblem subproblem.Solve(objective, delta); var Pstep = subproblem.Pstep; - var predictedReduction = subproblem.PredictedReduction; var hitBoundary = subproblem.HitBoundary; + // predicted reduction = L(0) - L(Δp) = Δp'g - 1/2 * Δp'HΔp + var predictedReduction = objective.Gradient.DotProduct(Pstep) - 0.5 * Pstep.DotProduct(objective.Hessian * Pstep); if (Pstep.L2Norm() <= stepTolerance * (stepTolerance + P.L2Norm())) { diff --git a/src/Numerics/Optimization/TrustRegionSubProblem.cs b/src/Numerics/Optimization/TrustRegionSubProblem.cs index aed6886b..cde8ecf8 100644 --- a/src/Numerics/Optimization/TrustRegionSubProblem.cs +++ b/src/Numerics/Optimization/TrustRegionSubProblem.cs @@ -4,9 +4,9 @@ namespace MathNet.Numerics.Optimization { public static class TrustRegionSubproblem { - public static ITrustRegionSubproblem Quadratic() + public static ITrustRegionSubproblem DogLeg() { - return new QuadraticSubproblem(); + return new DogLegSubproblem(); } } }