SQP improvements with new Test Classes

hipparchus-optim/nbactions.xml

@@ -0,0 +1,29 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<actions>
+        <action>
+            <actionName>test</actionName>
+            <packagings>
+                <packaging>*</packaging>
+            </packagings>
+            <goals>
+                <goal>test</goal>
+            </goals>
+            <properties>
+                <hipparchus.debug.sqp>false</hipparchus.debug.sqp>
+            </properties>
+        </action>
+        <action>
+            <actionName>test.single</actionName>
+            <packagings>
+                <packaging>*</packaging>
+            </packagings>
+            <goals>
+                <goal>process-test-classes</goal>
+                <goal>surefire:test</goal>
+            </goals>
+            <properties>
+                <test>${packageClassName}</test>
+                <hipparchus.debug.sqp>true</hipparchus.debug.sqp>
+            </properties>
+        </action>
+    </actions>

hipparchus-optim/src/main/java/org/hipparchus/optim/BaseMultivariateOptimizer.java

@@ -37,7 +37,7 @@
 public abstract class BaseMultivariateOptimizer<P>
     extends BaseOptimizer<P> {
     /** Initial guess. */
-    private double[] start;
+    public double[] start;
     /** Lower bounds. */
     private double[] lowerBound;
     /** Upper bounds. */

hipparchus-optim/src/main/java/org/hipparchus/optim/nonlinear/vector/constrained/AbstractSQPOptimizer2.java

@@ -21,6 +21,7 @@
 import org.hipparchus.linear.EigenDecompositionSymmetric;
 import org.hipparchus.optim.LocalizedOptimFormats;
 import org.hipparchus.optim.OptimizationData;
+import org.hipparchus.optim.SimpleBounds;
 import org.hipparchus.optim.nonlinear.scalar.ObjectiveFunction;
 import org.hipparchus.util.MathUtils;
 
@@ -49,7 +50,10 @@ public abstract class AbstractSQPOptimizer2 extends ConstraintOptimizer {
 
     /** Inequality constraint (may be null). */
     private BoundedConstraint boxConstraint;
-
+
+    /** Simple Bounds (may be null). */
+    private SimpleBounds simpleBounds;
+
     /** Default QPSolver. */
     private QPOptimizer QPSolver = new QPDualActiveSolver();
 
@@ -102,6 +106,13 @@ public InequalityConstraint getIqConstraint() {
     public BoundedConstraint getBoxConstraint() {
         return boxConstraint;
     }
+
+     /** Getter for simple bounds.
+     * @return simple bounds
+     */
+    public SimpleBounds getSimpleBounds() {
+        return simpleBounds;
+    }
 
     /** Getter for QP Solver.
      * @return QP Solver
@@ -120,6 +131,21 @@ public void parseOptimizationData(OptimizationData... optData) {
         super.parseOptimizationData(optData);
         for (OptimizationData data : optData) {
 
+            if (data instanceof SQPProblem) {
+                SQPProblem problem = (SQPProblem) data;
+                obj = new SQPObj(problem);
+
+                eqConstraint= (problem.hasEquality())?new SQPEq((SQPProblem) data):null;
+                iqConstraint= (problem.hasInequality())?new SQPIneq((SQPProblem) data):null;
+
+                double[] lb = ((SQPProblem) data).getBoxConstraintLB();
+                double[]ub = ((SQPProblem) data).getBoxConstraintUB();
+                simpleBounds = (problem.hasBounds())?new SimpleBounds(lb,ub):null;
+
+                start=(problem.hasInitialGuess())?problem.getInitialGuess():null;
+                continue;
+            }
+
             if (data instanceof ObjectiveFunction) {
                 obj = (TwiceDifferentiableFunction) ((ObjectiveFunction) data).getObjectiveFunction();
                 continue;
@@ -138,6 +164,11 @@ public void parseOptimizationData(OptimizationData... optData) {
                 boxConstraint = (BoundedConstraint) data;
                 continue;
             }
+
+            if (data instanceof SimpleBounds) {
+                simpleBounds = (SimpleBounds) data;
+                continue;
+            }
 
             if (data instanceof SQPOption) {
                 settings = (SQPOption) data;

hipparchus-optim/src/main/java/org/hipparchus/optim/nonlinear/vector/constrained/BFGSUpdater.java

@@ -16,6 +16,7 @@
  */
 package org.hipparchus.optim.nonlinear.vector.constrained;
 
+
 import org.hipparchus.linear.Array2DRowRealMatrix;
 import org.hipparchus.linear.ArrayRealVector;
 import org.hipparchus.linear.CholeskyDecomposition;
@@ -30,43 +31,67 @@
  * Manages Hessian updates for SQP solvers by:
  * </p>
  * <ul>
- *   <li>Checking curvature condition</li>
- *   <li>Applying dynamic damping if necessary</li>
- *   <li>Skipping update if curvature still fails after damping</li>
- *   <li>Soft regularization of diagonal entries on repeated failures</li>
- *   <li>Automatic Hessian reset after configurable failures</li>
+ * <li>Checking curvature condition</li>
+ * <li>Applying dynamic damping if necessary</li>
+ * <li>Skipping update if curvature still fails after damping</li>
+ * <li>Soft regularization of diagonal entries on repeated failures</li>
+ * <li>Automatic Hessian reset after configurable failures</li>
  * </ul>
  *
  * @since 4.1
  */
 public class BFGSUpdater {
 
-    /** Damping factor. */
+    /**
+     * AutoScaling Flag.
+     */
+    private final boolean SCALE;
+
+    /**
+     * EPS.
+     */
+    private final double EPS;
+
+    /**
+     * Damping factor.
+     */
     private static final double GAMMA = 0.2;
 
-    /** Regularization factor for diagonal of Hessian. */
-    private final double regFactor;
+    /**
+     * trigger skip update for diagonal of Hessian.
+     */
+    private final double sqrtEPSmachine = FastMath.sqrt(Precision.EPSILON);
 
-    /** Tolerance for symmetric matrices decomposition.
+    /**
+     * Tolerance for symmetric matrices decomposition.
+     *
      * @since 4.1
      */
     private final double decompositionEpsilon;
 
-    /** Stored initial Hessian for resets. */
+    /**
+     * Stored initial Hessian for resets.
+     */
     private final RealMatrix initialH;
 
-    /** Current Cholesky factor L such that H = L·Lᵀ. */
+    /**
+     * Current Cholesky factor L such that H = L·Lᵀ.
+     */
     private RealMatrix L;
 
     /**
      * Creates a new updater.
-     * @param initialHess          initial positive‐definite Hessian matrix
-     * @param regFactor            regularization factor to add on diagonal
-     * @param decompositionEpsilon tolerance for symmetric matrices decomposition
+     *
+     * @param initialHess initial positive‐definite Hessian matrix
+     * @param eps treshold to apply auto scale sty<sqrt(eps)
+     * @param autoSCale true apply auto hessain rescaling
+     * @param decompositionEpsilon tolerance for symmetric matrices
+     * decomposition
      */
-    public BFGSUpdater(final RealMatrix initialHess, final double regFactor, final double decompositionEpsilon) {
-        this.initialH             = new Array2DRowRealMatrix(initialHess.getData());
-        this.regFactor            = regFactor;
+    public BFGSUpdater(final RealMatrix initialHess, final double eps, final boolean autoScale, final double decompositionEpsilon) {
+        this.initialH = new Array2DRowRealMatrix(initialHess.getData());
+        this.EPS = eps;
+        this.SCALE = autoScale;
         this.decompositionEpsilon = decompositionEpsilon;
         resetHessian();
     }
@@ -84,20 +109,61 @@ public RealMatrix getHessian() {
      * Updates the Hessian approximation using the BFGS formula.
      * <p>
      * If curvature condition fails, applies damping or regularization.
-     *</p>
+     * </p>
      *
      * @param s displacement vector (x_{k+1} − x_k)
      * @param y1 gradient difference (∇f_{k+1} − ∇f_k)
+     * @return type of update
+     * 0: update is done
+     * 1:sHs too small skip update
+     * 2:sty<gamma*sHs skipped update
+     * 3:sty<sqrt(eps)auto scale gamma*Hini
+     * 4:singulatity detected during downdate reset to gamma*Hinit
      */
-    public void update(RealVector s, RealVector y1) {
-        RealVector y = damp(s, y1);
-        if (y == null) {
-            return;
+    public int update(RealVector s, RealVector y1) {
+        RealVector Hs = L.operate(L.preMultiply(s)); // al posto di getHessian().operate(s)
+        double sHs = s.dotProduct(Hs);
+        if (sHs <= 0.0) {
+
+//            double sty = s.dotProduct(y1);
+//            double yy = y1.dotProduct(y1);
+//            double gamma =yy /sty;
+//            if (sty>0) {
+//                gamma = FastMath.max(sqrtEPSmachine, FastMath.min(1/sqrtEPSmachine, gamma));
+//                this.resetHessian(gamma);
+//                return 1;
+//            }
+//            this.resetHessian();
+            return 1;
+        }
+
+        RealVector y = damp(s, y1, Hs, sHs);
+        if (y == null ) {
+
+            return 2;//damp failed
         }
+
+
+
+
+        if (SCALE ) {
+            double sty = s.dotProduct(y);
+            double yy = y.dotProduct(y);
+            double gamma =yy /sty;
+            if (gamma < FastMath.sqrt(EPS)) {
+                gamma = FastMath.max(sqrtEPSmachine, FastMath.min(1/sqrtEPSmachine, gamma));
+                this.resetHessian(gamma);
+                return 3;
+            }
+        }
+
         // Attempt rank‐one BFGS update; regularize on failure
-        rankOneUpdate(s, y);
+        if (rankOneUpdate(s, y, Hs, sHs))  return 0;
+
+        return 4;
     }
-
+
+
     /**
      * Resets the Hessian approximation to its initial value.
      */
@@ -106,54 +172,80 @@ public void resetHessian() {
         L = ch.getL();
     }
 
+    /**
+     * Resets the Hessian approximation with information on the curvature.
+     *
+     * @param gamma scale factor
+     */
+    public void resetHessian(double gamma) {
+        final CholeskyDecomposition ch = new CholeskyDecomposition(initialH.scalarMultiply(gamma), decompositionEpsilon, decompositionEpsilon);
+        L = ch.getL();
+    }
+
     /**
      * Applies dynamic damping to y to satisfy curvature condition sᵀy ≥ γ sᵀHs.
      *
      * @param s search direction
      * @param y1 raw gradient difference
+     * @param Hs vector
+     * @param sHs value
      * @return damped y, or null if update should be skipped
      */
-    public RealVector damp(RealVector s, RealVector y1) {
+    public RealVector damp(RealVector s, RealVector y1, RealVector Hs, double sHs) {
         RealVector y = new ArrayRealVector(y1);
         double sty = s.dotProduct(y1);
-        RealVector Hs = getHessian().operate(s);
-        double sHs = s.dotProduct(Hs);
-        if (sty <= Precision.EPSILON) {
-            return null;
-        }
+
         if (sty < GAMMA * sHs) {
             double phi = (1.0 - GAMMA) * sHs / (sHs - sty);
-            // clamp phi to [0,1]
-            phi = FastMath.max(0.0, FastMath.min(1.0, phi));
-            y = y1.mapMultiply(phi).add(Hs.mapMultiply(1.0 - phi));
+            if (phi >= 1.0) return y1;
+             else if (phi <= 0.0) return Hs;
+//            
+            y = (y1.mapMultiply(phi)).add(Hs.mapMultiply(1.0 - phi));
             sty = s.dotProduct(y);
+
+
             if (sty < GAMMA * sHs) {
+
                 return null;
             }
+
         }
         return y;
     }
+
+
 
     /**
      * Performs a BFGS rank‐one update on L.
      *
      * @param s displacement vector
      * @param y gradient difference vector
+     * @param Hs vector
+     * @param sHs value
      * @return true if update succeeded, false otherwise
      */
-    private boolean rankOneUpdate(RealVector s, RealVector y) {
+    private boolean rankOneUpdate(RealVector s, RealVector y, RealVector Hs, double sHs) {
         RealMatrix Lcopy = new Array2DRowRealMatrix(L.getData());
-        RealVector Hs = L.operate(L.preMultiply(s));
         double rho = 1.0 / FastMath.sqrt(s.dotProduct(y));
-        double theta = 1.0 / FastMath.sqrt(s.dotProduct(Hs));
+        double theta = 1.0 / FastMath.sqrt(sHs);
         RealVector v = y.mapMultiply(rho);
         RealVector w = Hs.mapMultiply(theta);
-        cholupdateLower(v, +1) ;
-
-        if (!cholupdateLower(w, -1)) {
-            //try to regularize
-            L.setSubMatrix(Lcopy.getData(), 0, 0);
+        cholupdateLower(v,+1);//upgrade
+        if (!cholupdateLower(w,-1)) { //downdate    
+            double gamma = 1.0;
+            double sty = s.dotProduct(y);
+            double yy = y.dotProduct(y);
+            if (sty> Precision.SAFE_MIN) {
+                gamma =yy /sty;
+                gamma = FastMath.max(sqrtEPSmachine, FastMath.min(1/sqrtEPSmachine, gamma));
+                this.resetHessian(gamma);
+            }
+            else this.resetHessian();
+
+
+            return false;
         }
+
         return true;
     }
 
@@ -174,9 +266,13 @@ private boolean cholupdateLower(RealVector u, int sigma) {
             double lii = L.getEntry(i, i);
             double ui = temp.getEntry(i);
             double r2 = lii * lii + sigma * ui * ui;
-            if (r2 < regFactor) {
+
+            //skip or update 
+            if (sigma < 0 && r2 < this.sqrtEPSmachine * lii * lii) {
+//                if (sigma < 0 && r2 < Precision.EPSILON) {
                 return false;
             }
+
             double r = Math.sqrt(r2);
             double c = r / lii;
             double s = ui / lii;
@@ -192,5 +288,9 @@ private boolean cholupdateLower(RealVector u, int sigma) {
         }
         return true;
     }
+
+
+
+
 
 }