Changeset 17294 for branches/2994-AutoDiffForIntervals

Timestamp:

10/02/19 17:06:47 (5 years ago)

Author:

gkronber

Message:

#2994: min() and max() are problematic in AutoDiff. We have used max(a,b) = (a+b + abs(b-a)) / 2 but this is problematic because of catastrophic cancellation. Min() and Max() are necessary for interval calculation. However, there we only need to determine min and max for double values (or duals). We do not need min() and max() in general.

File:

• branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs (modified) (13 diffs)

branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs

@@ -422 +422 @@
       // we assume that the for variable nodes ( v(x,w) = w * x ) the gradient is returned for parameter w
       instruction.value = new AlgebraicInterval(
-        low: new MultivariateDual<AlgebraicDouble>(v: lower, key: variable, dv: intervals[variable.VariableName].LowerBound),  
+        low: new MultivariateDual<AlgebraicDouble>(v: lower, key: variable, dv: intervals[variable.VariableName].LowerBound),
         high: new MultivariateDual<AlgebraicDouble>(v: upper, key: variable, dv: intervals[variable.VariableName].UpperBound)
         );
@@ -433 +433 @@
 
   public interface IAlgebraicType<T> {
-    T Zero { get; }
+    T Zero { get; } // Zero and One must create new objects
     T One { get; }
 
@@ -453 +453 @@
     T AssignIntRoot(T a, int r);
     T AssignSgn(T a); // set this to sign(a)
+    T AssignMin(T other); // set this min(this, other)
+    T AssignMax(T other); // set this max(this, other)
     T Clone();
   }
@@ -468 +470 @@
     public static T IntRoot<T>(this T a, int r) where T : IAlgebraicType<T> { a.AssignIntRoot(a.Clone(), r); return a; }
 
-    public static T Max<T>(T a, T b) where T : IAlgebraicType<T> {
-      // ((a + b) + abs(b - a)) / 2
-      return a.Clone().Add(b).Add(b.Clone().Sub(a).Abs()).Scale(1.0 / 2.0);
-    }
-    public static T Min<T>(T a, T b) where T : IAlgebraicType<T> {
-      // ((a + b) - abs(a - b)) / 2
-      return a.Clone().Add(b).Sub(a.Clone().Sub(b).Abs()).Scale(1.0 / 2.0);
-    }
+    internal static T Min<T>(T a, T b) where T : IAlgebraicType<T> { return a.Clone().AssignMin(b); }
+    internal static T Max<T>(T a, T b) where T : IAlgebraicType<T> { return a.Clone().AssignMax(b); }
+
+    // public static T Max<T>(T a, T b) where T : IAlgebraicType<T> {
+    //   // ((a + b) + abs(b - a)) / 2
+    //   return a.Clone().Add(b).Add(b.Clone().Sub(a).Abs()).Scale(1.0 / 2.0);
+    // }
+    // public static T Min<T>(T a, T b) where T : IAlgebraicType<T> {
+    //   // ((a + b) - abs(a - b)) / 2
+    //   return a.Clone().Add(b).Sub(a.Clone().Sub(b).Abs()).Scale(1.0 / 2.0);
+    // }
   }
 
@@ -511 +516 @@
     public AlgebraicDouble AssignIntPower(AlgebraicDouble a, int p) { Value = Math.Pow(a.Value, p); return this; }
     public AlgebraicDouble AssignIntRoot(AlgebraicDouble a, int r) { Value = IntRoot(a.Value, r); return this; }
+    public AlgebraicDouble AssignMin(AlgebraicDouble other) { Value = Math.Min(Value, other.Value); return this; }
+    public AlgebraicDouble AssignMax(AlgebraicDouble other) { Value = Math.Max(Value, other.Value); return this; }
 
     // helper
@@ -563 +570 @@
     public AlgebraicDoubleVector AssignIntPower(AlgebraicDoubleVector a, int p) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Pow(a.arr[i], p); } return this; }
     public AlgebraicDoubleVector AssignIntRoot(AlgebraicDoubleVector a, int r) { for (int i = 0; i < arr.Length; ++i) { arr[i] = IntRoot(a.arr[i], r); } return this; }
+    public AlgebraicDoubleVector AssignMin(AlgebraicDoubleVector other) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Min(arr[i], other.arr[i]); } return this; }
+    public AlgebraicDoubleVector AssignMax(AlgebraicDoubleVector other) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Max(arr[i], other.arr[i]); } return this; }
 
     // helper
@@ -778 +787 @@
     }
 
+    public AlgebraicSparseVector<K, T> AssignMin(AlgebraicSparseVector<K, T> other) {
+      // assumes that keys without a matching key in other are zero and vice versa
+      foreach (var kvp in elems) if (!other.elems.ContainsKey(kvp.Key)) kvp.Value.AssignMin(kvp.Value.Zero); // min(v, 0)
+      foreach (var kvp in other.elems) {
+        if (elems.TryGetValue(kvp.Key, out T value))
+          value.AssignMin(kvp.Value);
+        else
+          elems.Add(kvp.Key, kvp.Value.Zero.AssignMin(kvp.Value));
+      }
+      return this;
+    }
+
+    public AlgebraicSparseVector<K, T> AssignMax(AlgebraicSparseVector<K, T> other) {
+      // assumes that keys without a matching key in other are zero and vice versa
+      foreach (var kvp in elems) if (!other.elems.ContainsKey(kvp.Key)) kvp.Value.AssignMax(kvp.Value.Zero); // max(v, 0)
+      foreach (var kvp in other.elems) {
+        if (elems.TryGetValue(kvp.Key, out T value))
+          value.AssignMax(kvp.Value);
+        else
+          elems.Add(kvp.Key, kvp.Value.Zero.AssignMax(kvp.Value));
+      }
+      return this;
+    }
+
+
     public AlgebraicSparseVector<K, T> Clone() {
       return new AlgebraicSparseVector<K, T>(this);
@@ -836 +870 @@
     }
 
-    // algebraic min() / max() do not work for infinities 
-    // detect and handle infinite values explicitly
+
     private static MultivariateDual<AlgebraicDouble> Min(MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {
-      if (a.Value.IsInfinity || b.Value.IsInfinity) return a.Value < b.Value ? a : b; // NOTE: this is not differentiable but for infinite values we do not care
-      else return Algebraic.Min(a, b); // differentiable statement
+      return a.Value < b.Value ? a : b;
     }
     private static MultivariateDual<AlgebraicDouble> Max(MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {
-      if (a.Value.IsInfinity || b.Value.IsInfinity) return a.Value >= b.Value ? a : b; // NOTE: this is not differentiable but for infinite values we do not care
-      else return Algebraic.Max(a, b); // differentiable statement
+      return a.Value > b.Value ? a : b;
     }
 
@@ -919 +950 @@
         if (a.Contains(0.0)) {
           low = new MultivariateDual<AlgebraicDouble>(0.0);
-          high = Algebraic.Max(a.low.IntPower(p), a.high.IntPower(p));
+          high = a.low.IntPower(p).AssignMax(a.high.IntPower(p));
         } else {
           var lowPower = a.low.IntPower(p);
           var highPower = a.high.IntPower(p);
-          low = Algebraic.Min(lowPower, highPower);
-          high = Algebraic.Max(lowPower, highPower);
+          low = lowPower.AssignMin(highPower);
+          high = lowPower.AssignMax(highPower);
         }
       } else {
@@ -934 +965 @@
           var lowPower = a.low.IntPower(p);
           var highPower = a.high.IntPower(p);
-          low = Algebraic.Min(lowPower, highPower);
-          high = Algebraic.Max(lowPower, highPower);
+          low = lowPower.AssignMin(highPower);
+          high = lowPower.AssignMax(highPower);
         }
       }
@@ -1064 +1095 @@
       return this;
     }
+
+    public AlgebraicInterval AssignMin(AlgebraicInterval other) {
+      low.AssignMin(other.low);
+      high.AssignMin(other.high);
+      return this;
+    }
+
+    public AlgebraicInterval AssignMax(AlgebraicInterval other) {
+      low.AssignMax(other.low);
+      high.AssignMax(other.high);
+      return this;
+    }
   }
 
@@ -1116 +1159 @@
     public Dual<V> Clone() { return new Dual<V>(v.Clone(), dv.Clone()); }
 
+    public Dual<V> AssignMin(Dual<V> other) {
+      throw new NotImplementedException();
+    }
+
+    public Dual<V> AssignMax(Dual<V> other) {
+      throw new NotImplementedException();
+    }
   }
 
@@ -1199 +1249 @@
     public MultivariateDual<V> AssignSgn(MultivariateDual<V> a) { v.AssignSgn(a.v); dv = a.dv.Zero; return this; } // sign(f(x))' = 0;     
 
+    public MultivariateDual<V> AssignMin(MultivariateDual<V> other) {
+      XXX
+    }
+
+    public MultivariateDual<V> AssignMax(MultivariateDual<V> other) {
+      XXX
+    }
   }
 }