Changeset 17204

Timestamp:

08/13/19 09:29:11 (5 years ago)

Author:

gkronber

Message:

#2994: added parameter for gradient checks and experimented with preconditioning

Location:

branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions

Files:

• ConstrainedNLS.cs (modified) (7 diffs)
• ConstrainedNLSInternal.cs (modified) (17 diffs)
• NLOpt.cs (modified) (1 diff)

branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions/ConstrainedNLS.cs

@@ -38 +38 @@
       get { return (IFixedValueParameter<DoubleValue>)Parameters["MaxTime"]; }
     }
+    public IFixedValueParameter<BoolValue> CheckGradientParameter {
+      get { return (IFixedValueParameter<BoolValue>)Parameters["CheckGradient"]; }
+    }
     public int Iterations { get { return IterationsParameter.Value.Value; } set { IterationsParameter.Value.Value = value; } }
 
@@ -48 +51 @@
       get { return ModelStructureParameter.Value.Value; }
       set { ModelStructureParameter.Value.Value = value; }
+    }
+    public bool CheckGradient {
+      get { return CheckGradientParameter.Value.Value; }
+      set { CheckGradientParameter.Value.Value = value; }
     }
 
@@ -64 +71 @@
       Parameters.Add(new FixedValueParameter<DoubleValue>("FuncToleranceAbs", new DoubleValue(0)));
       Parameters.Add(new FixedValueParameter<DoubleValue>("MaxTime", new DoubleValue(10)));
+      Parameters.Add(new FixedValueParameter<BoolValue>("CheckGradient", "Flag to indicate whether the gradient should be checked using numeric approximation", new BoolValue(false)));
+
+      CheckGradientParameter.Hidden = true;
     }
 
     public ConstrainedNLS(ConstrainedNLS original, Cloner cloner) : base(original, cloner) {
+    }
+
+    [StorableHook(HookType.AfterDeserialization)]
+    public void AfterDeserializationHook() {
+      if (!Parameters.ContainsKey("CheckGradient")) {
+        Parameters.Add(new FixedValueParameter<BoolValue>("CheckGradient", "Flag to indicate whether the gradient should be checked using numeric approximation", new BoolValue(false)));
+
+        CheckGradientParameter.Hidden = true;
+      }
     }
 
@@ -89 +108 @@
       Results.AddOrUpdateResult("Evaluations", functionEvaluations);
       var bestError = new DoubleValue(double.MaxValue);
+      var curError = new DoubleValue(double.MaxValue);
       Results.AddOrUpdateResult("Best error", bestError);
+      Results.AddOrUpdateResult("Current error", curError);
       Results.AddOrUpdateResult("Tree", tree);
       var qualitiesTable = new DataTable("Qualities");
@@ -112 +133 @@
 
       var state = new ConstrainedNLSInternal(Solver.Value, tree, Iterations, problem.ProblemData, FuncToleranceRel, FuncToleranceAbs, MaxTime);
+      if (CheckGradient) state.CheckGradient = true;
+      int idx = 0;
+      foreach(var constraintTree in state.constraintTrees) {
+        Results.AddOrUpdateResult($"Constraint {idx++}", constraintTree);
+      }
 
       // we use a listener model here to get state from the solver      
@@ -123 +149 @@
       bestQualityRow.Values.Add(bestError.Value);
 
-
-      Results.AddOrUpdateResult("Best solution", CreateSolution(state.BestTree, problem.ProblemData));
+      Results.AddOrUpdateResult("Best solution", CreateSolution((ISymbolicExpressionTree)state.BestTree.Clone(), problem.ProblemData));
+      Results.AddOrUpdateResult("Best solution constraint values", new DoubleArray(state.BestConstraintValues));
 
 
@@ -132 +158 @@
         functionEvaluations.Value++;
         bestError.Value = state.BestError;
+        curError.Value = state.CurError;
         curQualityRow.Values.Add(state.CurError);
         bestQualityRow.Values.Add(bestError.Value);

branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions/ConstrainedNLSInternal.cs

@@ -38 +38 @@
     public ISymbolicExpressionTree BestTree => bestTree;
 
+    private double[] bestConstraintValues;
+    public double[] BestConstraintValues => bestConstraintValues;
+
+    public bool CheckGradient { get; internal set; }
+
     // begin internal state
     private IntPtr nlopt;
     private SymbolicDataAnalysisExpressionTreeLinearInterpreter interpreter;
     private readonly NLOpt.nlopt_func calculateObjectiveDelegate; // must hold the delegate to prevent GC
+    private readonly NLOpt.nlopt_precond preconditionDelegate;
     private readonly IntPtr[] constraintDataPtr; // must hold the objects to prevent GC
     private readonly NLOpt.nlopt_func[] calculateConstraintDelegates; // must hold the delegates to prevent GC
@@ -51 +57 @@
     private readonly ISymbolicExpressionTreeNode[] preparedTreeParameterNodes;
     private readonly List<ConstantTreeNode>[] allThetaNodes;
+    public List<ISymbolicExpressionTree> constraintTrees;    // TODO make local in ctor (public for debugging)
+    
     private readonly double[] fi_eval;
     private readonly double[,] jac_eval;
     private readonly ISymbolicExpressionTree scaledTree;
+    private readonly VectorAutoDiffEvaluator autoDiffEval;
+    private readonly VectorEvaluator eval;
 
     // end internal state
@@ -72 +82 @@
       this.problemData = problemData;
       this.interpreter = new SymbolicDataAnalysisExpressionTreeLinearInterpreter();
-
+      this.autoDiffEval = new VectorAutoDiffEvaluator();
+      this.eval = new VectorEvaluator();
+
+      CheckGradient = false;
 
       var intervalConstraints = problemData.IntervalConstraints;
@@ -103 +116 @@
       Dictionary<string, ISymbolicExpressionTree> derivatives = new Dictionary<string, ISymbolicExpressionTree>();
       allThetaNodes = thetaNames.Select(_ => new List<ConstantTreeNode>()).ToArray();
-      var constraintTrees = new List<ISymbolicExpressionTree>();
+      constraintTrees = new List<ISymbolicExpressionTree>();
       foreach (var constraint in intervalConstraints.Constraints) {
         if (constraint.IsDerivation) {
@@ -117 +130 @@
             // convert variables named theta back to constants
             var df_prepared = ReplaceVarWithConst(df_smaller_upper, thetaNames, thetaValues, allThetaNodes);
-            constraintTrees.Add(df_prepared);
+            constraintTrees.Add((ISymbolicExpressionTree)df_prepared.Clone());
           }
           if (constraint.Interval.LowerBound > double.NegativeInfinity) {
             var df_larger_lower = Subtract(CreateConstant(constraint.Interval.LowerBound), (ISymbolicExpressionTree)df.Clone());
             // convert variables named theta back to constants
-            var df_prepared = ReplaceVarWithConst(df_larger_lower, thetaNames, thetaValues, allThetaNodes);
-            constraintTrees.Add(df_prepared);
+            var df_prepared = ReplaceVarWithConst(df_larger_lower, thetaNames, thetaValues, allThetaNodes);            
+            constraintTrees.Add((ISymbolicExpressionTree)df_prepared.Clone());
           }
         } else {
@@ -130 +143 @@
             // convert variables named theta back to constants
             var df_prepared = ReplaceVarWithConst(f_smaller_upper, thetaNames, thetaValues, allThetaNodes);
-            constraintTrees.Add(df_prepared);
+            constraintTrees.Add((ISymbolicExpressionTree)df_prepared.Clone());
           }
           if (constraint.Interval.LowerBound > double.NegativeInfinity) {
@@ -136 +149 @@
             // convert variables named theta back to constants
             var df_prepared = ReplaceVarWithConst(f_larger_lower, thetaNames, thetaValues, allThetaNodes);
-            constraintTrees.Add(df_prepared);
+            constraintTrees.Add((ISymbolicExpressionTree)df_prepared.Clone());
           }
         }
@@ -149 +162 @@
 
 
-      var minVal = Math.Min(-1000.0, thetaValues.Min());
-      var maxVal = Math.Max(1000.0, thetaValues.Max());
+      var minVal = Math.Min(-100.0, thetaValues.Min());
+      var maxVal = Math.Max(100.0, thetaValues.Max());
       var lb = Enumerable.Repeat(minVal, thetaValues.Count).ToArray();
       var up = Enumerable.Repeat(maxVal, thetaValues.Count).ToArray();
@@ -158 +171 @@
       NLOpt.nlopt_set_upper_bounds(nlopt, up);
       calculateObjectiveDelegate = new NLOpt.nlopt_func(CalculateObjective); // keep a reference to the delegate (see below)
-      NLOpt.nlopt_set_min_objective(nlopt, calculateObjectiveDelegate, IntPtr.Zero);
+      NLOpt.nlopt_set_min_objective(nlopt, calculateObjectiveDelegate, IntPtr.Zero); // --> without preconditioning
+
+      // preconditionDelegate = new NLOpt.nlopt_precond(PreconditionObjective);
+      // NLOpt.nlopt_set_precond_min_objective(nlopt, calculateObjectiveDelegate, preconditionDelegate, IntPtr.Zero);
 
 
@@ -169 +185 @@
         calculateConstraintDelegates[i] = new NLOpt.nlopt_func(CalculateConstraint);
         NLOpt.nlopt_add_inequality_constraint(nlopt, calculateConstraintDelegates[i], constraintDataPtr[i], 1e-8);
+        //NLOpt.nlopt_add_precond_inequality_constraint(nlopt, calculateConstraintDelegates[i], preconditionDelegate, constraintDataPtr[i], 1e-8);
       }
 
@@ -176 +193 @@
       NLOpt.nlopt_set_maxeval(nlopt, maxIterations);
     }
+
 
     ~ConstrainedNLSInternal() {
@@ -195 +213 @@
       bestError = minf;
 
-      if (res < 0) {
+      if (res < 0 && res != NLOpt.nlopt_result.NLOPT_FORCED_STOP) {
         throw new InvalidOperationException($"NLOpt failed {res} {NLOpt.nlopt_get_errmsg(nlopt)}");
       } else {
+        // calculate constraints of final solution
+        double[] _ = new double[x.Length];
+        bestConstraintValues = new double[calculateConstraintDelegates.Length];
+        for(int i=0;i<calculateConstraintDelegates.Length;i++) {
+          bestConstraintValues[i] = calculateConstraintDelegates[i].Invoke((uint)x.Length, x, _, constraintDataPtr[i]);
+        }
+
         // update parameters in tree
         var pIdx = 0;
@@ -218 +243 @@
 
       if (grad != null) {
-        var autoDiffEval = new VectorAutoDiffEvaluator();
         autoDiffEval.Evaluate(preparedTree, problemData.Dataset, trainingRows,
           preparedTreeParameterNodes, fi_eval, jac_eval);
@@ -235 +259 @@
 
         #region check gradient
-        var eval = new VectorEvaluator();
-        for (int i = 0; i < dim; i++) {
-          // make two additional evaluations
-          var xForNumDiff = (double[])curX.Clone();
-          double delta = Math.Abs(xForNumDiff[i] * 1e-5);
-          xForNumDiff[i] += delta;
-          UpdateThetaValues(xForNumDiff);
-          var evalHigh = eval.Evaluate(preparedTree, problemData.Dataset, trainingRows);
-          var mseHigh = MSE(target, evalHigh);
-          xForNumDiff[i] = curX[i] - delta;
-          UpdateThetaValues(xForNumDiff);
-          var evalLow = eval.Evaluate(preparedTree, problemData.Dataset, trainingRows);
-          var mseLow = MSE(target, evalLow);
-
-          var numericDiff = (mseHigh - mseLow) / (2 * delta);
-          var autoDiff = grad[i];
-          if ((Math.Abs(autoDiff) < 1e-10 && Math.Abs(numericDiff) > 1e-2)
-            || (Math.Abs(autoDiff) >= 1e-10 && Math.Abs((numericDiff - autoDiff) / numericDiff) > 1e-2))
-            throw new InvalidProgramException();
+        if (grad != null && CheckGradient) {
+          for (int i = 0; i < dim; i++) {
+            // make two additional evaluations
+            var xForNumDiff = (double[])curX.Clone();
+            double delta = Math.Abs(xForNumDiff[i] * 1e-5);
+            xForNumDiff[i] += delta;
+            UpdateThetaValues(xForNumDiff);
+            var evalHigh = eval.Evaluate(preparedTree, problemData.Dataset, trainingRows);
+            var mseHigh = MSE(target, evalHigh);
+            xForNumDiff[i] = curX[i] - delta;
+            UpdateThetaValues(xForNumDiff);
+            var evalLow = eval.Evaluate(preparedTree, problemData.Dataset, trainingRows);
+            var mseLow = MSE(target, evalLow);
+
+            var numericDiff = (mseHigh - mseLow) / (2 * delta);
+            var autoDiff = grad[i];
+            if ((Math.Abs(autoDiff) < 1e-10 && Math.Abs(numericDiff) > 1e-2)
+              || (Math.Abs(autoDiff) >= 1e-10 && Math.Abs((numericDiff - autoDiff) / numericDiff) > 1e-2))
+              throw new InvalidProgramException();
+          }
         }
         #endregion
@@ -275 +300 @@
     }
 
+    // // NOT WORKING YET? WHY is det(H) always zero?!
+    // // TODO
+    // private void PreconditionObjective(uint n, double[] x, double[] v, double[] vpre, IntPtr data) {
+    //   UpdateThetaValues(x); // calc H(x)
+    //   
+    //   autoDiffEval.Evaluate(preparedTree, problemData.Dataset, trainingRows,
+    //     preparedTreeParameterNodes, fi_eval, jac_eval);
+    //   var k = jac_eval.GetLength(0);
+    //   var h = new double[n, n];
+    //   
+    //   // calc residuals and scale jac_eval
+    //   var f = -2.0 / k;
+    //   for (int i = 0;i<k;i++) {
+    //     var r = target[i] - fi_eval[i];
+    //     for(int j = 0;j<n;j++) {
+    //       jac_eval[i, j] *= f * r;
+    //     }
+    //   }
+    //   
+    //   // approximate hessian H(x) = J(x)^T * J(x)
+    //   alglib.rmatrixgemm((int)n, (int)n, k,
+    //     1.0, jac_eval, 0, 0, 1,  // transposed
+    //     jac_eval, 0, 0, 0, 
+    //     0.0, ref h, 0, 0,
+    //     null
+    //     );
+    //   
+    //   // scale v
+    //   alglib.rmatrixmv((int)n, (int)n, h, 0, 0, 0, v, 0, ref vpre, 0, alglib.serial);
+    //   var det = alglib.matdet.rmatrixdet(h, (int)n, alglib.serial);
+    // }
+
+
     private double MSE(double[] a, double[] b) {
       Trace.Assert(a.Length == b.Length);
@@ -288 +346 @@
         bestSolution = (double[])curX.Clone();
       }
+      curError = curF;
     }
 
@@ -313 +372 @@
 
       #region check gradient
-      for (int i = 0; i < dim; i++) {
-        // make two additional evaluations
-        var xForNumDiff = (double[])curX.Clone();
-        double delta = Math.Abs(xForNumDiff[i] * 1e-5);
-        xForNumDiff[i] += delta;
-        UpdateThetaValues(xForNumDiff);
-        var evalHigh = intervalEvaluator.Evaluate(constraintData.Tree, dataIntervals, constraintData.ParameterNodes,
-          out double[] unusedLowerGradientHigh, out double[] unusedUpperGradientHigh);
-
-        xForNumDiff[i] = curX[i] - delta;
-        UpdateThetaValues(xForNumDiff);
-        var evalLow = intervalEvaluator.Evaluate(constraintData.Tree, dataIntervals, constraintData.ParameterNodes,
-          out double[] unusedLowerGradientLow, out double[] unusedUpperGradientLow);
-        
-        var numericDiff = (evalHigh.UpperBound - evalLow.UpperBound) / (2 * delta);
-        var autoDiff = grad[i];
-        // XXX GRADIENT FOR UPPER BOUND FOR FUNCTION IS ZERO FOR THE FIRST SET OF VARIABLES? GRADIENT FOR ADDITION INCORRECT?!
-        if ((Math.Abs(autoDiff) < 1e-10 && Math.Abs(numericDiff) > 1e-2)
-          || (Math.Abs(autoDiff) >= 1e-10 && Math.Abs((numericDiff - autoDiff) / numericDiff) > 1e-2))
-          throw new InvalidProgramException();
-      }
+      if (grad != null && CheckGradient)
+        for (int i = 0; i < dim; i++) {
+          // make two additional evaluations
+          var xForNumDiff = (double[])curX.Clone();
+          double delta = Math.Abs(xForNumDiff[i] * 1e-5);
+          xForNumDiff[i] += delta;
+          UpdateThetaValues(xForNumDiff);
+          var evalHigh = intervalEvaluator.Evaluate(constraintData.Tree, dataIntervals, constraintData.ParameterNodes,
+            out double[] unusedLowerGradientHigh, out double[] unusedUpperGradientHigh);
+
+          xForNumDiff[i] = curX[i] - delta;
+          UpdateThetaValues(xForNumDiff);
+          var evalLow = intervalEvaluator.Evaluate(constraintData.Tree, dataIntervals, constraintData.ParameterNodes,
+            out double[] unusedLowerGradientLow, out double[] unusedUpperGradientLow);
+
+          var numericDiff = (evalHigh.UpperBound - evalLow.UpperBound) / (2 * delta);
+          var autoDiff = grad[i];
+
+          if ((Math.Abs(autoDiff) < 1e-10 && Math.Abs(numericDiff) > 1e-2)
+            || (Math.Abs(autoDiff) >= 1e-10 && Math.Abs((numericDiff - autoDiff) / numericDiff) > 1e-2))
+            throw new InvalidProgramException();
+        }
       #endregion
 

branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions/NLOpt.cs

@@ -18 +18 @@
     /* A preconditioner, which preconditions v at x to return vpre.
        (The meaning of "preconditioning" is algorithm-dependent.) */
-    public delegate IntPtr nlopt_precond(uint n,
+    public delegate void nlopt_precond(uint n,
       [In] [MarshalAs(UnmanagedType.LPArray, SizeParamIndex = 0)] double[] x, 
       [In] [MarshalAs(UnmanagedType.LPArray, SizeParamIndex = 0)] double[] v,