Changeset 8420 for branches/NCA/HeuristicLab.Algorithms.NCA/3.3/NeighborhoodComponentsAnalysis.cs

Timestamp:

08/07/12 00:42:33 (12 years ago)

Author:

abeham

Message:

#1913:

• Worked on NCA
• Added scatter plot view for the model to show training data when it is reduced to two dimensions

It works, but I don't think it works correctly yet. I have randomized the initial matrix, because the starting point influences the achievable quality quite a bit.

Location:

branches/NCA

Files:

• . (modified) (1 prop)
• HeuristicLab.Algorithms.NCA/3.3 (modified) (1 prop)
• HeuristicLab.Algorithms.NCA/3.3/NeighborhoodComponentsAnalysis.cs (modified) (4 diffs)

branches/NCA/HeuristicLab.Algorithms.NCA/3.3/NeighborhoodComponentsAnalysis.cs

@@ -34 +34 @@
       var matrix = new double[attributes * reduceDimensions];
       // TODO: make some more clever initialization of matrix
+      var rand = new Random();
+      var counter = 0;
       for (var i = 0; i < attributes; i++)
         for (var j = 0; j < reduceDimensions; j++)
-          matrix[j * attributes + i] = 1;
+          matrix[counter++] = rand.NextDouble();
 
       alglib.mincgstate state;
@@ -48 +50 @@
 
       var transformationMatrix = new double[attributes, reduceDimensions];
+      counter = 0;
       for (var i = 0; i < attributes; i++)
         for (var j = 0; j < reduceDimensions; j++)
-          transformationMatrix[i, j] = matrix[j * attributes + i];
+          transformationMatrix[i, j] = matrix[counter++];
 
       var transformedTrainingset = new double[instances, reduceDimensions];
       var rowCount = 0;
-      for (var r = data.TrainingPartition.Start; r < data.TrainingPartition.End; r++) {
-        if (r >= data.TestPartition.Start && r < data.TestPartition.End) continue;
+      foreach (var r in data.TrainingIndices) {
         var j = 0;
         foreach (var v in data.AllowedInputVariables) {
@@ -77 +79 @@
       var info = (OptimizationInfo)obj;
       var instances = info.ProblemData.TrainingIndices.ToArray();
+      var AMatrix = new Matrix(A, A.Length / info.ReduceDimensions, info.ReduceDimensions);
       var distances = new double[instances.Length, instances.Length];
+
       for (int i = 0; i < instances.Length - 1; i++) {
+        var iVector = new Matrix(GetRow(info.ProblemData, instances[i]));
         for (int j = i + 1; j < instances.Length; j++) {
-          distances[i, j] = VectorLength(VectorSubtract(
-            Ax(A, GetRow(info.ProblemData, instances[i])),
-            Ax(A, GetRow(info.ProblemData, instances[j]))));
+          var jVector = new Matrix(GetRow(info.ProblemData, instances[j]));
+          distances[i, j] = iVector.Multiply(AMatrix).Subtract(jVector.Multiply(AMatrix)).QuadraticSum();
           distances[j, i] = distances[i, j];
         }
       }
-      var probabilities = new double[instances.Length, instances.Length];
-      for (int i = 0; i < instances.Length - 1; i++) {
+      alglib.sparsematrix probabilities;
+      alglib.sparsecreate(instances.Length, instances.Length, out probabilities);
+      for (int i = 0; i < instances.Length; i++) {
         var denom = 0.0;
         for (int k = 0; k < instances.Length; k++) {
@@ -93 +98 @@
           denom += Math.Exp(-(distances[i, k] * distances[i, k]));
         }
-        for (int j = i + 1; j < instances.Length; j++) {
-          probabilities[i, j] = Math.Exp(-(distances[i, j] * distances[i, j])) / denom;
-        }
-      }
-
-      var target = info.ProblemData.TargetVariable;
-      var ds = info.ProblemData.Dataset;
-      var pi = new double[instances.Length];
-      func = 0.0;
-      for (int i = 0; i < instances.Length; i++) {
-        var classI = ds.GetDoubleValue(target, instances[i]);
-        pi[i] = 0.0;
-        for (int j = 0; j < instances.Length; j++) {
-          if (i == j) continue;
-          var classJ = ds.GetDoubleValue(target, instances[j]);
-          if (classI.IsAlmost(classJ)) {
-            func += probabilities[i, j];
-            pi[i] += probabilities[i, j];
+        if (denom > 0) {
+          for (int j = 0; j < instances.Length; j++) {
+            if (i == j) continue;
+            var v = Math.Exp(-(distances[i, j] * distances[i, j])) / denom;
+            alglib.sparseset(probabilities, i, j, v);
           }
         }
       }
+      alglib.sparseconverttocrs(probabilities); // needed to enumerate in order
 
-      var sum = GetSum(probabilities, pi, info.ProblemData, instances);
-      for (int i = 0; i < A.Length; i++) {
-        grad[i] = 2.0 * A[i] * sum;
+      int t0 = 0, t1 = 0, r, c;
+      double val;
+      var classes = info.ProblemData.Dataset.GetDoubleValues(info.ProblemData.TargetVariable, instances).ToArray();
+      var pi = new double[instances.Length];
+      while (alglib.sparseenumerate(probabilities, ref t0, ref t1, out r, out c, out val)) {
+        if (classes[r].IsAlmost(classes[c]))
+          pi[r] += val;
       }
+      func = pi.Sum();
+
+      grad = AMatrix.Multiply(2.0).Transpose().Multiply(CalculateInnerSum(probabilities, pi, info.ProblemData, instances, classes)).Transpose().ToArray();
     }
 
-    private static double[,] GetSum(double[,] p, double[] pi, IClassificationProblemData data, int[] instances) {
+    private static Matrix CalculateInnerSum(alglib.sparsematrix p, double[] pi, IClassificationProblemData data, int[] instances, double[] classes) {
+      var attributes = data.AllowedInputVariables.Count();
       var target = data.TargetVariable;
-      var sum = new double[pi.Length, pi.Length];
-      for (int i = 0; i < pi.Length; i++) {
-        var classI = data.Dataset.GetDoubleValue(target, instances[i]);
-        var sumK = 0.0;
-        for (int k = 0; k < pi.Length; k++) {
-          if (i == k) continue;
-          sumK += p[i, k] * MatrixMultiply(VectorSubtract(GetRow(data, i), GetRow(data, k)), VectorSubtract(GetRow(data, i), GetRow(data, k)));
-        }
-        var sumJ = 0.0;
-        for (int j = 0; j < pi.Length; j++) {
-          if (i == j) continue;
-          var classJ = data.Dataset.GetDoubleValue(target, instances[j]);
-          if (classI.IsAlmost(classJ)) sumJ += p[i, j] * ScalarProduct(VectorSubtract(GetRow(data, i), GetRow(data, j)), VectorSubtract(GetRow(data, i), GetRow(data, j)));
-        }
-        sum += pi[i] * sumK - sumJ;
-      }
-      return sum;
-    }
+      int t0 = 0, t1 = 0, r, c;
+      double v;
+      var result = new double[attributes, attributes];
+      while (alglib.sparseenumerate(p, ref t0, ref t1, out r, out c, out v)) {
+        var vector = new Matrix(GetRow(data, instances[r])).Subtract(new Matrix(GetRow(data, instances[c]))).Apply();
+        vector.OuterProduct(vector).Multiply(v * pi[r]).AddTo(result);
 
-    private static double[,] MatrixMultiply(IEnumerable<double> vec1, IEnumerable<double> vec2) {
-      var enumVec1 = vec1.GetEnumerator();
-      var enumVec2 = vec2.GetEnumerator();
-      bool f1, f2;
-      while (true) {
-        f1 = enumVec1.MoveNext();
-        f2 = enumVec2.MoveNext();
-        if (!f1 && f2 || f1 && !f2) throw new ArgumentException("vectors are of unequal length.");
-        if (!f1 && !f2) yield break;
-        yield return enumVec1.Current - enumVec2.Current;
-      }
-    }
-
-    private static IEnumerable<double> Ax(IEnumerable<double> A, IEnumerable<double> x) {
-      var enumA = A.GetEnumerator();
-      var enumX = x.GetEnumerator();
-      var sum = 0.0;
-      while (true) {
-        if (!enumA.MoveNext()) {
-          yield return sum;
-          yield break;
-        } else if (!enumX.MoveNext()) {
-          yield return sum;
-          enumX = x.GetEnumerator();
-          sum = 0.0;
-        } else {
-          sum += enumA.Current * enumX.Current;
+        if (classes[r].IsAlmost(classes[c])) {
+          vector.OuterProduct(vector).Multiply(-v).AddTo(result);
         }
       }
+      return new Matrix(result);
     }
 
     private static IEnumerable<double> GetRow(IClassificationProblemData data, int row) {
       return data.AllowedInputVariables.Select(v => data.Dataset.GetDoubleValue(v, row));
-    }
-
-    private static IEnumerable<double> VectorSubtract(IEnumerable<double> vec1, IEnumerable<double> vec2) {
-      var enumVec1 = vec1.GetEnumerator();
-      var enumVec2 = vec2.GetEnumerator();
-      bool f1, f2;
-      while (true) {
-        f1 = enumVec1.MoveNext();
-        f2 = enumVec2.MoveNext();
-        if (!f1 && f2 || f1 && !f2) throw new ArgumentException("vectors are of unequal length.");
-        if (!f1 && !f2) yield break;
-        yield return enumVec1.Current - enumVec2.Current;
-      }
-    }
-
-    private static double ScalarProduct(IEnumerable<double> vec1, IEnumerable<double> vec2) {
-      var enumVec1 = vec1.GetEnumerator();
-      var enumVec2 = vec2.GetEnumerator();
-      bool f1, f2;
-      double product = 0.0;
-      while (true) {
-        f1 = enumVec1.MoveNext();
-        f2 = enumVec2.MoveNext();
-        if (!f1 && f2 || f1 && !f2) throw new ArgumentException("vectors are of unequal length.");
-        if (!f1 && !f2) return product;
-        product += enumVec1.Current * enumVec2.Current;
-      }
-    }
-
-    private static double VectorLength(IEnumerable<double> vector) {
-      return Math.Sqrt(vector.Sum(x => x * x));
     }