Recovered (and deprecated) BFGSUpdater.

Author: maisonobe

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

@@ -16,7 +16,6 @@
  */
 package org.hipparchus.optim.nonlinear.vector.constrained;
 
-
 import org.hipparchus.linear.Array2DRowRealMatrix;
 import org.hipparchus.linear.ArrayRealVector;
 import org.hipparchus.linear.CholeskyDecomposition;
@@ -38,60 +37,37 @@
  *   <li>Automatic Hessian reset after configurable failures</li>
  * </ul>
  *
- * @since 4.1
+ * @deprecated as of 4.1, replaced by {@link BFGSUpdater2}
  */
+@Deprecated
 public class BFGSUpdater {
 
-    /**
-     * Damping factor.
-     */
+    /** Damping factor. */
     private static final double GAMMA = 0.2;
 
-    /**
-     * AutoScaling Flag.
-     */
-    private final boolean scale;
-
-    /**
-     * EPS.
-     */
-    private final double eps;
+    /** 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 eps treshold to apply auto scale sty<sqrt(eps)
-     * @param autoScale true apply auto hessian rescaling
-     * @param decompositionEpsilon tolerance for symmetric matrices
-     * decomposition
+     * @param initialHess          initial positive‐definite Hessian matrix
+     * @param regFactor            regularization factor to add on diagonal
+     * @param decompositionEpsilon tolerance for symmetric matrices decomposition
      */
-    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;
+    public BFGSUpdater(final RealMatrix initialHess, final double regFactor, final double decompositionEpsilon) {
+        this.initialH             = new Array2DRowRealMatrix(initialHess.getData());
+        this.regFactor            = regFactor;
         this.decompositionEpsilon = decompositionEpsilon;
         resetHessian();
     }
@@ -109,61 +85,20 @@ 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 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
+    public void update(RealVector s, RealVector y1) {
+        RealVector y = damp(s, y1);
+        if (y == null) {
+            return;
         }
-
-
-
-
-        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
-        if (rankOneUpdate(s, y, Hs, sHs))  return 0;
-
-        return 4;
+        rankOneUpdate(s, y);
     }
 
-
     /**
      * Resets the Hessian approximation to its initial value.
      */
@@ -172,80 +107,54 @@ 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, RealVector Hs, double sHs) {
+    public RealVector damp(RealVector s, RealVector y1) {
         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);
-            if (phi >= 1.0) return y1;
-             else if (phi <= 0.0) return Hs;
-//
-            y = (y1.mapMultiply(phi)).add(Hs.mapMultiply(1.0 - phi));
+            // 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));
             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, RealVector Hs, double sHs) {
+    private boolean rankOneUpdate(RealVector s, RealVector y) {
         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(sHs);
+        double theta = 1.0 / FastMath.sqrt(s.dotProduct(Hs));
         RealVector v = y.mapMultiply(rho);
         RealVector w = Hs.mapMultiply(theta);
-        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();
-
+        cholupdateLower(v, +1) ;
 
-            return false;
+        if (!cholupdateLower(w, -1)) {
+            //try to regularize
+            L.setSubMatrix(Lcopy.getData(), 0, 0);
         }
-
         return true;
     }
 
@@ -266,13 +175,9 @@ 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;
-
-            //skip or update
-            if (sigma < 0 && r2 < this.sqrtEPSmachine * lii * lii) {
-//                if (sigma < 0 && r2 < Precision.EPSILON) {
+            if (r2 < regFactor) {
                 return false;
             }
-
             double r = Math.sqrt(r2);
             double c = r / lii;
             double s = ui / lii;
@@ -289,8 +194,4 @@ private boolean cholupdateLower(RealVector u, int sigma) {
         return true;
     }
 
-
-
-
-
 }