Changeset 8681

Timestamp:

09/23/12 00:31:49 (12 years ago)

Author:

abeham

Message:

#1913: Added regularization term introduced by Yang and Laaksonen 2007

Location:

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/Nca

Files:

• Matrix.cs (modified) (3 diffs)
• NcaAlgorithm.cs (modified) (11 diffs)

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/Nca/Matrix.cs

@@ -70 +70 @@
 
     public Matrix Transpose() {
-      var result = new Matrix(Transpose(values, Columns, Rows), Columns, Rows);
-      return result;
+      return new Matrix(Transpose(values, Columns, Rows), Columns, Rows);
     }
 
@@ -125 +124 @@
     }
 
-    public double VectorLength() {
-      return Math.Sqrt(SquaredVectorLength());
-    }
-
-    public double SquaredVectorLength() {
-      if (Rows != 1) throw new ArgumentException("Length only works on vectors.");
+    public double EuclideanNorm() {
+      return Math.Sqrt(SumOfSquares());
+    }
+
+    public double SumOfSquares() {
       return values.Sum(x => x * x);
     }
@@ -137 +135 @@
       if (Rows != 1 || other.Rows != 1) throw new ArgumentException("OuterProduct can only be applied to vectors.");
       return Transpose().Multiply(other);
+    }
+
+    public IEnumerable<double> ColumnSums() {
+      return Transpose().RowSums();
+    }
+
+    public IEnumerable<double> RowSums() {
+      var sum = 0.0;
+      int counter = 0;
+      foreach (var v in values) {
+        sum += v;
+        counter++;
+        if (counter == Rows) {
+          yield return sum;
+          sum = 0.0;
+          counter = 0;
+        }
+      }
     }
 

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/Nca/NcaAlgorithm.cs

@@ -35 +35 @@
 
 namespace HeuristicLab.Algorithms.DataAnalysis {
-  internal delegate void Reporter(double quality, double[] coefficients);
+  internal delegate void Reporter(double quality, double[] coefficients, double[] gradients);
   /// <summary>
   /// Neighborhood Components Analysis
   /// </summary>
-  [Item("Neighborhood Components Analysis (NCA)", "Implementation of Neighborhood Components Analysis based on the description of J. Goldberger, S. Roweis, G. Hinton, R. Salakhutdinov. 2005. Neighbourhood Component Analysis. Advances in Neural Information Processing Systems, 17. pp. 513-520.")]
+  [Item("Neighborhood Components Analysis (NCA)", @"Implementation of Neighborhood Components Analysis
+based on the description of J. Goldberger, S. Roweis, G. Hinton, R. Salakhutdinov. 2005.
+Neighbourhood Component Analysis. Advances in Neural Information Processing Systems, 17. pp. 513-520
+with additional regularizations described in Z. Yang, J. Laaksonen. 2007.
+Regularized Neighborhood Component Analysis. Lecture Notes in Computer Science, 4522. pp. 253-262.")]
   [Creatable("Data Analysis")]
   [StorableClass]
@@ -59 +63 @@
       get { return (IFixedValueParameter<IntValue>)Parameters["Iterations"]; }
     }
+    public IFixedValueParameter<DoubleValue> RegularizationParameter {
+      get { return (IFixedValueParameter<DoubleValue>)Parameters["Regularization"]; }
+    }
     #endregion
 
     #region Properties
-    private int K {
+    public int K {
       get { return KParameter.Value.Value; }
       set { KParameter.Value.Value = value; }
     }
-    private int Dimensions {
+    public int Dimensions {
       get { return DimensionsParameter.Value.Value; }
       set { DimensionsParameter.Value.Value = value; }
     }
-    private int NeighborSamples {
+    public int NeighborSamples {
       get { return NeighborSamplesParameter.Value.Value; }
       set { NeighborSamplesParameter.Value.Value = value; }
     }
-    private int Iterations {
+    public int Iterations {
       get { return IterationsParameter.Value.Value; }
       set { IterationsParameter.Value.Value = value; }
+    }
+    public double Regularization {
+      get { return RegularizationParameter.Value.Value; }
+      set { RegularizationParameter.Value.Value = value; }
     }
     #endregion
@@ -85 +96 @@
     public NcaAlgorithm()
       : base() {
-      Parameters.Add(new FixedValueParameter<IntValue>("K", "The K for the nearest neighbor.", new IntValue(1)));
+      Parameters.Add(new FixedValueParameter<IntValue>("K", "The K for the nearest neighbor.", new IntValue(3)));
       Parameters.Add(new FixedValueParameter<IntValue>("Dimensions", "The number of dimensions that NCA should reduce the data to.", new IntValue(2)));
       Parameters.Add(new ConstrainedValueParameter<INCAInitializer>("Initialization", "Which method should be used to initialize the matrix. Typically LDA (linear discriminant analysis) should provide a good estimate."));
-      Parameters.Add(new FixedValueParameter<IntValue>("NeighborSamples", "How many of the neighbors should be sampled in order to speed up the calculation. This should be at least the value of k and at most the number of training instances minus one.", new IntValue(50)));
-      Parameters.Add(new FixedValueParameter<IntValue>("Iterations", "How many iterations the conjugate gradient (CG) method should be allowed to perform. The method might still terminate earlier if a local optima has already been reached.", new IntValue(20)));
+      Parameters.Add(new FixedValueParameter<IntValue>("NeighborSamples", "How many of the neighbors should be sampled in order to speed up the calculation. This should be at least the value of k and at most the number of training instances minus one.", new IntValue(60)));
+      Parameters.Add(new FixedValueParameter<IntValue>("Iterations", "How many iterations the conjugate gradient (CG) method should be allowed to perform. The method might still terminate earlier if a local optima has already been reached.", new IntValue(50)));
+      Parameters.Add(new FixedValueParameter<DoubleValue>("Regularization", "A non-negative paramter which can be set to increase generalization and avoid overfitting. If set to 0 the algorithm is similar to NCA as proposed by Goldberger et al.", new DoubleValue(0)));
 
       INCAInitializer defaultInitializer = null;
@@ -105 +117 @@
     }
 
+    [StorableHook(HookType.AfterDeserialization)]
+    private void AfterDeserialization() {
+      if (!Parameters.ContainsKey("Regularization")) {
+        Parameters.Add(new FixedValueParameter<DoubleValue>("Regularization", "A non-negative paramter which can be set to increase generalization and avoid overfitting. If set to 0 the algorithm is similar to NCA as proposed by Goldberger et al.", new DoubleValue(0)));
+      }
+    }
+
     public override void Prepare() {
       if (Problem != null) base.Prepare();
@@ -113 +132 @@
 
       var clonedProblem = (IClassificationProblemData)Problem.ProblemData.Clone();
-      var model = Train(clonedProblem, K, Dimensions, NeighborSamples, Iterations, initializer.Initialize(clonedProblem, Dimensions), ReportQuality, CancellationToken.None);
-      Results.Add(new Result("ClassificationSolution", "The classification solution.", model.CreateClassificationSolution(clonedProblem)));
-    }
-
-    public static INcaClassificationSolution CreateClassificationSolution(IClassificationProblemData data, int k, int dimensions, int neighborSamples, int iterations, INCAInitializer initializer) {
+      var model = Train(clonedProblem, K, Dimensions, NeighborSamples, Regularization, Iterations, initializer.Initialize(clonedProblem, Dimensions), ReportQuality, CancellationToken.None);
+      var solution = model.CreateClassificationSolution(clonedProblem);
+      if (!Results.ContainsKey("ClassificationSolution"))
+        Results.Add(new Result("ClassificationSolution", "The classification solution.", solution));
+      else Results["ClassificationSolution"].Value = solution;
+    }
+
+    public static INcaClassificationSolution CreateClassificationSolution(IClassificationProblemData data, int k, int dimensions, int neighborSamples, double regularization, int iterations, INCAInitializer initializer) {
       var clonedProblem = (IClassificationProblemData)data.Clone();
-      var model = Train(clonedProblem, k, dimensions, neighborSamples, iterations, initializer);
+      var model = Train(clonedProblem, k, dimensions, neighborSamples, regularization, iterations, initializer);
       return model.CreateClassificationSolution(clonedProblem);
     }
 
-    public static INcaModel Train(IClassificationProblemData problemData, int k, int dimensions, int neighborSamples, int iterations, INCAInitializer initializer) {
-      return Train(problemData, k, dimensions, neighborSamples, iterations, initializer.Initialize(problemData, dimensions), null, CancellationToken.None);
-    }
-
-    public static INcaModel Train(IClassificationProblemData problemData, int k, int neighborSamples, int iterations, double[,] initalMatrix) {
+    public static INcaModel Train(IClassificationProblemData problemData, int k, int dimensions, int neighborSamples, double regularization, int iterations, INCAInitializer initializer) {
+      return Train(problemData, k, dimensions, neighborSamples, regularization, iterations, initializer.Initialize(problemData, dimensions), null, CancellationToken.None);
+    }
+
+    public static INcaModel Train(IClassificationProblemData problemData, int k, int neighborSamples, double regularization, int iterations, double[,] initalMatrix) {
       var matrix = new double[initalMatrix.Length];
       for (int i = 0; i < initalMatrix.GetLength(0); i++)
         for (int j = 0; j < initalMatrix.GetLength(1); j++)
           matrix[i * initalMatrix.GetLength(1) + j] = initalMatrix[i, j];
-      return Train(problemData, k, initalMatrix.GetLength(1), neighborSamples, iterations, matrix, null, CancellationToken.None);
-    }
-
-    private static INcaModel Train(IClassificationProblemData data, int k, int dimensions, int neighborSamples, int iterations, double[] matrix, Reporter reporter, CancellationToken cancellation) {
+      return Train(problemData, k, initalMatrix.GetLength(1), neighborSamples, regularization, iterations, matrix, null, CancellationToken.None);
+    }
+
+    private static INcaModel Train(IClassificationProblemData data, int k, int dimensions, int neighborSamples, double regularization, int iterations, double[] matrix, Reporter reporter, CancellationToken cancellation) {
       var scaling = new Scaling(data.Dataset, data.AllowedInputVariables, data.TrainingIndices);
       var scaledData = AlglibUtil.PrepareAndScaleInputMatrix(data.Dataset, data.AllowedInputVariables, data.TrainingIndices, scaling);
       var classes = data.Dataset.GetDoubleValues(data.TargetVariable, data.TrainingIndices).ToArray();
       var attributes = scaledData.GetLength(1);
-
-      var penalties = new Dictionary<double, Dictionary<double, double>>();
-      foreach (var c in data.ClassValues) {
-        penalties[c] = new Dictionary<double, double>();
-        foreach (var r in data.ClassValues)
-          penalties[c][r] = data.GetClassificationPenalty(c, r);
-      }
 
       alglib.mincgstate state;
@@ -153 +168 @@
       alglib.mincgsetcond(state, 0, 0, 0, iterations);
       alglib.mincgsetxrep(state, true);
+      //alglib.mincgsetgradientcheck(state, 0.01);
       int neighborSampleSize = neighborSamples;
-      Optimize(state, scaledData, classes, penalties, dimensions, neighborSampleSize, cancellation, reporter);
+      Optimize(state, scaledData, classes, dimensions, neighborSampleSize, regularization, cancellation, reporter);
       alglib.mincgresults(state, out matrix, out rep);
+      if (rep.terminationtype == -7) throw new InvalidOperationException("Gradient verification failed.");
 
       var transformationMatrix = new double[attributes, dimensions];
@@ -166 +183 @@
     }
 
-    private static void Optimize(alglib.mincgstate state, double[,] data, double[] classes, Dictionary<double, Dictionary<double, double>> penalties, int dimensions, int neighborSampleSize, CancellationToken cancellation, Reporter reporter) {
+    private static void Optimize(alglib.mincgstate state, double[,] data, double[] classes, int dimensions, int neighborSampleSize, double lambda, CancellationToken cancellation, Reporter reporter) {
       while (alglib.mincgiteration(state)) {
         if (cancellation.IsCancellationRequested) break;
         if (state.needfg) {
-          Gradient(state.x, ref state.innerobj.f, state.innerobj.g, data, classes, penalties, dimensions, neighborSampleSize);
+          Gradient(state.x, ref state.innerobj.f, state.innerobj.g, data, classes, dimensions, neighborSampleSize, lambda);
           continue;
         }
         if (state.innerobj.xupdated) {
           if (reporter != null)
-            reporter(state.innerobj.f, state.innerobj.x);
+            reporter(state.innerobj.f, state.innerobj.x, state.innerobj.g);
           continue;
         }
@@ -182 +199 @@
     }
 
-    private static void Gradient(double[] A, ref double func, double[] grad, double[,] data, double[] classes, Dictionary<double, Dictionary<double, double>> penalties, int dimensions, int neighborSampleSize) {
+    private static void Gradient(double[] A, ref double func, double[] grad, double[,] data, double[] classes, int dimensions, int neighborSampleSize, double lambda) {
       var instances = data.GetLength(0);
       var attributes = data.GetLength(1);
@@ -199 +216 @@
           }
           var kVector = new Matrix(GetRow(data, k));
-          transformedDistances[k] = Math.Exp(-iVector.Multiply(AMatrix).Subtract(kVector.Multiply(AMatrix)).SquaredVectorLength());
-        }
-        var sample = transformedDistances.OrderByDescending(x => x.Value).Take(neighborSampleSize).ToArray();
-        var normalization = sample.Sum(x => x.Value);
-        if (normalization > 0) {
-          foreach (var s in sample) {
-            if (s.Value <= 0) break;
-            alglib.sparseset(probabilities, i, s.Key, s.Value / normalization);
-          }
+          transformedDistances[k] = Math.Exp(-iVector.Multiply(AMatrix).Subtract(kVector.Multiply(AMatrix)).SumOfSquares());
+        }
+        var normalization = transformedDistances.Sum(x => x.Value);
+        if (normalization <= 0) continue;
+        foreach (var s in transformedDistances.Where(x => x.Value > 0).OrderByDescending(x => x.Value).Take(neighborSampleSize)) {
+          alglib.sparseset(probabilities, i, s.Key, s.Value / normalization);
         }
       }
@@ -215 +229 @@
       double val;
       var pi = new double[instances];
-      func = 0;
       while (alglib.sparseenumerate(probabilities, ref t0, ref t1, out r, out c, out val)) {
-        double vp = val * penalties[classes[r]][classes[c]];
-        pi[r] += vp;
-        func += vp;
-      }
-
-      t0 = 0; t1 = 0;
+        if (classes[r].IsAlmost(classes[c])) {
+          pi[r] += val;
+        }
+      }
+
       var innerSum = new double[attributes, attributes];
       while (alglib.sparseenumerate(probabilities, ref t0, ref t1, out r, out c, out val)) {
         var vector = new Matrix(GetRow(data, r)).Subtract(new Matrix(GetRow(data, c)));
         vector.OuterProduct(vector).Multiply(val * pi[r]).AddTo(innerSum);
-        vector.OuterProduct(vector).Multiply(-val * penalties[classes[r]][classes[c]]).AddTo(innerSum);
-      }
+
+        if (classes[r].IsAlmost(classes[c])) {
+          vector.OuterProduct(vector).Multiply(-val).AddTo(innerSum);
+        }
+      }
+
+      func = -pi.Sum() + lambda * AMatrix.SumOfSquares();
 
       r = 0;
-      var newGrad = AMatrix.Multiply(2.0).Transpose().Multiply(new Matrix(innerSum)).Transpose();
+      var newGrad = AMatrix.Multiply(-2.0).Transpose().Multiply(new Matrix(innerSum)).Transpose();
       foreach (var g in newGrad) {
-        grad[r++] = g;
-      }
-    }
-
-    private void ReportQuality(double func, double[] coefficients) {
+        grad[r] = g + lambda * 2 * A[r];
+        r++;
+      }
+    }
+
+    private void ReportQuality(double func, double[] coefficients, double[] gradients) {
       var instances = Problem.ProblemData.TrainingIndices.Count();
       DataTable qualities;
       if (!Results.ContainsKey("Optimization")) {
         qualities = new DataTable("Optimization");
-        qualities.Rows.Add(new DataRow("Penalty", string.Empty));
+        qualities.Rows.Add(new DataRow("Quality", string.Empty));
         Results.Add(new Result("Optimization", qualities));
       } else qualities = (DataTable)Results["Optimization"].Value;
-      qualities.Rows["Penalty"].Values.Add(func / instances);
-
-      if (!Results.ContainsKey("Penalty")) {
-        Results.Add(new Result("Penalty", new DoubleValue(func / instances)));
-      } else ((DoubleValue)Results["Penalty"].Value).Value = func / instances;
+      qualities.Rows["Quality"].Values.Add(-func / instances);
+
+      string[] attributNames = Problem.ProblemData.AllowedInputVariables.ToArray();
+      if (gradients != null) {
+        DataTable grads;
+        if (!Results.ContainsKey("Gradients")) {
+          grads = new DataTable("Gradients");
+          for (int i = 0; i < gradients.Length; i++)
+            grads.Rows.Add(new DataRow(attributNames[i / Dimensions] + "-" + (i % Dimensions), string.Empty));
+          Results.Add(new Result("Gradients", grads));
+        } else grads = (DataTable)Results["Gradients"].Value;
+        for (int i = 0; i < gradients.Length; i++)
+          grads.Rows[attributNames[i / Dimensions] + "-" + (i % Dimensions)].Values.Add(gradients[i]);
+      }
+
+      if (!Results.ContainsKey("Quality")) {
+        Results.Add(new Result("Quality", new DoubleValue(-func / instances)));
+      } else ((DoubleValue)Results["Quality"].Value).Value = -func / instances;
+
+      var attributes = attributNames.Length;
+      var transformationMatrix = new double[attributes, Dimensions];
+      var counter = 0;
+      for (var i = 0; i < attributes; i++)
+        for (var j = 0; j < Dimensions; j++)
+          transformationMatrix[i, j] = coefficients[counter++];
+
+      var scaling = new Scaling(Problem.ProblemData.Dataset, attributNames, Problem.ProblemData.TrainingIndices);
+      var model = new NcaModel(K, transformationMatrix, Problem.ProblemData.Dataset, Problem.ProblemData.TrainingIndices, Problem.ProblemData.TargetVariable, attributNames, scaling, Problem.ProblemData.ClassValues.ToArray());
+
+      IClassificationSolution solution = model.CreateClassificationSolution(Problem.ProblemData);
+      if (!Results.ContainsKey("ClassificationSolution")) {
+        Results.Add(new Result("ClassificationSolution", solution));
+      } else {
+        Results["ClassificationSolution"].Value = solution;
+      }
     }
 
@@ -256 +304 @@
         yield return data[row, i];
     }
+
   }
 }