Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs @ 14923

Last change on this file since 14923 was 14862, checked in by gkronber, 8 years ago

#2700 copied tSNE from branch to trunk

File size: 28.4 KB
RevLine 
[14414]1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 *
5 * This file is part of HeuristicLab.
6 *
7 * HeuristicLab is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * HeuristicLab is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
19 */
20
21//Code is based on an implementation from Laurens van der Maaten
22
23/*
24*
25* Copyright (c) 2014, Laurens van der Maaten (Delft University of Technology)
26* All rights reserved.
27*
28* Redistribution and use in source and binary forms, with or without
29* modification, are permitted provided that the following conditions are met:
30* 1. Redistributions of source code must retain the above copyright
31*    notice, this list of conditions and the following disclaimer.
32* 2. Redistributions in binary form must reproduce the above copyright
33*    notice, this list of conditions and the following disclaimer in the
34*    documentation and/or other materials provided with the distribution.
35* 3. All advertising materials mentioning features or use of this software
36*    must display the following acknowledgement:
37*    This product includes software developed by the Delft University of Technology.
38* 4. Neither the name of the Delft University of Technology nor the names of
39*    its contributors may be used to endorse or promote products derived from
40*    this software without specific prior written permission.
41*
42* THIS SOFTWARE IS PROVIDED BY LAURENS VAN DER MAATEN ''AS IS'' AND ANY EXPRESS
43* OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
44* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
45* EVENT SHALL LAURENS VAN DER MAATEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
46* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
47* PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
48* BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
49* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
50* IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
51* OF SUCH DAMAGE.
52*
53*/
54#endregion
55
56using System;
57using System.Collections.Generic;
58using System.Linq;
[14785]59using HeuristicLab.Collections;
[14414]60using HeuristicLab.Common;
61using HeuristicLab.Core;
[14807]62using HeuristicLab.Optimization;
[14414]63using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
64using HeuristicLab.Random;
65
66namespace HeuristicLab.Algorithms.DataAnalysis {
67  [StorableClass]
[14807]68  public class TSNEStatic<T> {
[14414]69
[14788]70    [StorableClass]
71    public sealed class TSNEState : DeepCloneable {
72      // initialized once
[14806]73      [Storable]
[14788]74      public IDistance<T> distance;
[14806]75      [Storable]
[14788]76      public IRandom random;
[14806]77      [Storable]
[14788]78      public double perplexity;
[14806]79      [Storable]
[14788]80      public bool exact;
[14806]81      [Storable]
[14788]82      public int noDatapoints;
[14806]83      [Storable]
[14788]84      public double finalMomentum;
[14806]85      [Storable]
[14788]86      public int momSwitchIter;
[14806]87      [Storable]
[14788]88      public int stopLyingIter;
[14806]89      [Storable]
[14788]90      public double theta;
[14806]91      [Storable]
[14788]92      public double eta;
[14806]93      [Storable]
[14788]94      public int newDimensions;
[14414]95
[14788]96      // for approximate version: sparse representation of similarity/distance matrix
[14806]97      [Storable]
[14788]98      public double[] valP; // similarity/distance
[14806]99      [Storable]
[14788]100      public int[] rowP; // row index
[14806]101      [Storable]
[14788]102      public int[] colP; // col index
[14414]103
[14788]104      // for exact version: dense representation of distance/similarity matrix
[14806]105      [Storable]
[14788]106      public double[,] p;
[14512]107
[14788]108      // mapped data
[14806]109      [Storable]
[14788]110      public double[,] newData;
[14414]111
[14806]112      [Storable]
[14788]113      public int iter;
[14806]114      [Storable]
[14788]115      public double currentMomentum;
[14414]116
[14788]117      // helper variables (updated in each iteration)
[14806]118      [Storable]
[14788]119      public double[,] gains;
[14806]120      [Storable]
[14788]121      public double[,] uY;
[14806]122      [Storable]
[14788]123      public double[,] dY;
[14512]124
[14788]125      private TSNEState(TSNEState original, Cloner cloner) : base(original, cloner) {
[14806]126        this.distance = cloner.Clone(original.distance);
127        this.random = cloner.Clone(original.random);
128        this.perplexity = original.perplexity;
129        this.exact = original.exact;
130        this.noDatapoints = original.noDatapoints;
131        this.finalMomentum = original.finalMomentum;
132        this.momSwitchIter = original.momSwitchIter;
133        this.stopLyingIter = original.stopLyingIter;
134        this.theta = original.theta;
135        this.eta = original.eta;
136        this.newDimensions = original.newDimensions;
137        if(original.valP != null) {
138          this.valP = new double[original.valP.Length];
139          Array.Copy(original.valP, this.valP, this.valP.Length);
140        }
141        if(original.rowP != null) {
142          this.rowP = new int[original.rowP.Length];
143          Array.Copy(original.rowP, this.rowP, this.rowP.Length);
144        }
145        if(original.colP != null) {
146          this.colP = new int[original.colP.Length];
147          Array.Copy(original.colP, this.colP, this.colP.Length);
148        }
149        if(original.p != null) {
150          this.p = new double[original.p.GetLength(0), original.p.GetLength(1)];
151          Array.Copy(original.p, this.p, this.p.Length);
152        }
153        this.newData = new double[original.newData.GetLength(0), original.newData.GetLength(1)];
154        Array.Copy(original.newData, this.newData, this.newData.Length);
155        this.iter = original.iter;
156        this.currentMomentum = original.currentMomentum;
157        this.gains = new double[original.gains.GetLength(0), original.gains.GetLength(1)];
158        Array.Copy(original.gains, this.gains, this.gains.Length);
159        this.uY = new double[original.uY.GetLength(0), original.uY.GetLength(1)];
160        Array.Copy(original.uY, this.uY, this.uY.Length);
161        this.dY = new double[original.dY.GetLength(0), original.dY.GetLength(1)];
162        Array.Copy(original.dY, this.dY, this.dY.Length);
[14788]163      }
[14806]164
[14788]165      public override IDeepCloneable Clone(Cloner cloner) {
166        return new TSNEState(this, cloner);
167      }
[14414]168
[14807]169      [StorableConstructor]
[14837]170      public TSNEState(bool deserializing) { }
[14788]171      public TSNEState(T[] data, IDistance<T> distance, IRandom random, int newDimensions, double perplexity, double theta, int stopLyingIter, int momSwitchIter, double momentum, double finalMomentum, double eta) {
172        this.distance = distance;
173        this.random = random;
174        this.newDimensions = newDimensions;
175        this.perplexity = perplexity;
176        this.theta = theta;
177        this.stopLyingIter = stopLyingIter;
178        this.momSwitchIter = momSwitchIter;
179        this.currentMomentum = momentum;
180        this.finalMomentum = finalMomentum;
181        this.eta = eta;
[14414]182
183
[14788]184        // initialize
185        noDatapoints = data.Length;
[14806]186        if(noDatapoints - 1 < 3 * perplexity)
187          throw new ArgumentException("Perplexity too large for the number of data points!");
[14788]188
189        exact = Math.Abs(theta) < double.Epsilon;
190        newData = new double[noDatapoints, newDimensions];
191        dY = new double[noDatapoints, newDimensions];
192        uY = new double[noDatapoints, newDimensions];
193        gains = new double[noDatapoints, newDimensions];
[14806]194        for(var i = 0; i < noDatapoints; i++)
195          for(var j = 0; j < newDimensions; j++)
[14788]196            gains[i, j] = 1.0;
197
198        p = null;
199        rowP = null;
200        colP = null;
201        valP = null;
202
203        //Calculate Similarities
[14858]204        if (exact) p = CalculateExactSimilarites(data, distance, perplexity);
[14788]205        else CalculateApproximateSimilarities(data, distance, perplexity, out rowP, out colP, out valP);
206
[14837]207        // Lie about the P-values (factor is 4 in the MATLAB implementation)
[14806]208        if(exact) for(var i = 0; i < noDatapoints; i++) for(var j = 0; j < noDatapoints; j++) p[i, j] *= 12.0;
209        else for(var i = 0; i < rowP[noDatapoints]; i++) valP[i] *= 12.0;
[14788]210
211        // Initialize solution (randomly)
212        var rand = new NormalDistributedRandom(random, 0, 1);
[14806]213        for(var i = 0; i < noDatapoints; i++)
214          for(var j = 0; j < newDimensions; j++)
[14837]215            newData[i, j] = rand.NextDouble() * .0001;
[14414]216      }
217
[14788]218      public double EvaluateError() {
[14806]219        return exact ?
220          EvaluateErrorExact(p, newData, noDatapoints, newDimensions) :
221          EvaluateErrorApproximate(rowP, colP, valP, newData, theta);
[14788]222      }
[14512]223
[14788]224      private static void CalculateApproximateSimilarities(T[] data, IDistance<T> distance, double perplexity, out int[] rowP, out int[] colP, out double[] valP) {
225        // Compute asymmetric pairwise input similarities
[14837]226        ComputeGaussianPerplexity(data, distance, out rowP, out colP, out valP, perplexity, (int)(3 * perplexity));
[14788]227        // Symmetrize input similarities
228        int[] sRowP, symColP;
229        double[] sValP;
230        SymmetrizeMatrix(rowP, colP, valP, out sRowP, out symColP, out sValP);
231        rowP = sRowP;
232        colP = symColP;
233        valP = sValP;
234        var sumP = .0;
[14806]235        for(var i = 0; i < rowP[data.Length]; i++) sumP += valP[i];
236        for(var i = 0; i < rowP[data.Length]; i++) valP[i] /= sumP;
[14788]237      }
[14806]238
[14788]239      private static double[,] CalculateExactSimilarites(T[] data, IDistance<T> distance, double perplexity) {
240        // Compute similarities
241        var p = new double[data.Length, data.Length];
242        ComputeGaussianPerplexity(data, distance, p, perplexity);
243        // Symmetrize input similarities
[14806]244        for(var n = 0; n < data.Length; n++) {
245          for(var m = n + 1; m < data.Length; m++) {
[14788]246            p[n, m] += p[m, n];
247            p[m, n] = p[n, m];
248          }
249        }
250        var sumP = .0;
[14806]251        for(var i = 0; i < data.Length; i++) for(var j = 0; j < data.Length; j++) sumP += p[i, j];
252        for(var i = 0; i < data.Length; i++) for(var j = 0; j < data.Length; j++) p[i, j] /= sumP;
[14788]253        return p;
254      }
[14742]255
[14788]256      private static void ComputeGaussianPerplexity(IReadOnlyList<T> x, IDistance<T> distance, out int[] rowP, out int[] colP, out double[] valP, double perplexity, int k) {
[14806]257        if(perplexity > k) throw new ArgumentException("Perplexity should be lower than k!");
[14512]258
[14788]259        int n = x.Count;
260        // Allocate the memory we need
261        rowP = new int[n + 1];
262        colP = new int[n * k];
263        valP = new double[n * k];
264        var curP = new double[n - 1];
265        rowP[0] = 0;
[14806]266        for(var i = 0; i < n; i++) rowP[i + 1] = rowP[i] + k;
[14512]267
[14788]268        var objX = new List<IndexedItem<T>>();
[14806]269        for(var i = 0; i < n; i++) objX.Add(new IndexedItem<T>(i, x[i]));
[14512]270
[14788]271        // Build ball tree on data set
[14837]272        var tree = new VantagePointTree<IndexedItem<T>>(new IndexedItemDistance<T>(distance), objX);
[14742]273
[14788]274        // Loop over all points to find nearest neighbors
[14806]275        for(var i = 0; i < n; i++) {
[14788]276          IList<IndexedItem<T>> indices;
277          IList<double> distances;
[14742]278
[14788]279          // Find nearest neighbors
280          tree.Search(objX[i], k + 1, out indices, out distances);
[14512]281
[14788]282          // Initialize some variables for binary search
283          var found = false;
284          var beta = 1.0;
285          var minBeta = double.MinValue;
286          var maxBeta = double.MaxValue;
[14855]287          const double tol = 1e-5; 
[14512]288
[14788]289          // Iterate until we found a good perplexity
290          var iter = 0; double sumP = 0;
[14855]291          while(!found && iter < 200) { 
[14512]292
[14788]293            // Compute Gaussian kernel row
[14855]294            for(var m = 0; m < k; m++) curP[m] = Math.Exp(-beta * distances[m + 1]);
[14512]295
[14788]296            // Compute entropy of current row
297            sumP = double.Epsilon;
[14806]298            for(var m = 0; m < k; m++) sumP += curP[m];
[14788]299            var h = .0;
[14855]300            for(var m = 0; m < k; m++) h += beta * (distances[m + 1] * curP[m]);
[14788]301            h = h / sumP + Math.Log(sumP);
302
303            // Evaluate whether the entropy is within the tolerance level
304            var hdiff = h - Math.Log(perplexity);
[14806]305            if(hdiff < tol && -hdiff < tol) {
[14788]306              found = true;
307            } else {
[14806]308              if(hdiff > 0) {
[14788]309                minBeta = beta;
[14806]310                if(maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
[14788]311                  beta *= 2.0;
312                else
313                  beta = (beta + maxBeta) / 2.0;
314              } else {
315                maxBeta = beta;
[14806]316                if(minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
[14788]317                  beta /= 2.0;
318                else
319                  beta = (beta + minBeta) / 2.0;
320              }
321            }
322
323            // Update iteration counter
324            iter++;
325          }
326
327          // Row-normalize current row of P and store in matrix
[14806]328          for(var m = 0; m < k; m++) curP[m] /= sumP;
329          for(var m = 0; m < k; m++) {
[14788]330            colP[rowP[i] + m] = indices[m + 1].Index;
331            valP[rowP[i] + m] = curP[m];
332          }
[14512]333        }
334      }
[14788]335      private static void ComputeGaussianPerplexity(T[] x, IDistance<T> distance, double[,] p, double perplexity) {
336        // Compute the distance matrix
337        var dd = ComputeDistances(x, distance);
338
339        int n = x.Length;
340        // Compute the Gaussian kernel row by row
[14806]341        for(var i = 0; i < n; i++) {
[14788]342          // Initialize some variables
343          var found = false;
344          var beta = 1.0;
[14837]345          var minBeta = double.MinValue;
[14788]346          var maxBeta = double.MaxValue;
347          const double tol = 1e-5;
348          double sumP = 0;
349
350          // Iterate until we found a good perplexity
351          var iter = 0;
[14837]352          while(!found && iter < 200) {      // 200 iterations as in tSNE implementation by van der Maarten
[14788]353
354            // Compute Gaussian kernel row
[14806]355            for(var m = 0; m < n; m++) p[i, m] = Math.Exp(-beta * dd[i][m]);
[14788]356            p[i, i] = double.Epsilon;
357
358            // Compute entropy of current row
359            sumP = double.Epsilon;
[14806]360            for(var m = 0; m < n; m++) sumP += p[i, m];
[14788]361            var h = 0.0;
[14806]362            for(var m = 0; m < n; m++) h += beta * (dd[i][m] * p[i, m]);
[14788]363            h = h / sumP + Math.Log(sumP);
364
365            // Evaluate whether the entropy is within the tolerance level
366            var hdiff = h - Math.Log(perplexity);
[14806]367            if(hdiff < tol && -hdiff < tol) {
[14788]368              found = true;
369            } else {
[14806]370              if(hdiff > 0) {
[14788]371                minBeta = beta;
[14806]372                if(maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
[14788]373                  beta *= 2.0;
374                else
375                  beta = (beta + maxBeta) / 2.0;
376              } else {
377                maxBeta = beta;
[14806]378                if(minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
[14788]379                  beta /= 2.0;
380                else
381                  beta = (beta + minBeta) / 2.0;
382              }
383            }
384
385            // Update iteration counter
386            iter++;
387          }
388
389          // Row normalize P
[14806]390          for(var m = 0; m < n; m++) p[i, m] /= sumP;
[14512]391        }
392      }
393
[14788]394      private static double[][] ComputeDistances(T[] x, IDistance<T> distance) {
[14806]395        var res = new double[x.Length][];
396        for(int r = 0; r < x.Length; r++) {
397          var rowV = new double[x.Length];
398          // all distances must be symmetric
399          for(int c = 0; c < r; c++) {
400            rowV[c] = res[c][r];
401          }
402          rowV[r] = 0.0; // distance to self is zero for all distances
403          for(int c = r + 1; c < x.Length; c++) {
404            rowV[c] = distance.Get(x[r], x[c]);
405          }
406          res[r] = rowV;
407        }
408        return res;
409        // return x.Select(m => x.Select(n => distance.Get(m, n)).ToArray()).ToArray();
[14788]410      }
[14414]411
[14788]412      private static double EvaluateErrorExact(double[,] p, double[,] y, int n, int d) {
413        // Compute the squared Euclidean distance matrix
414        var dd = new double[n, n];
415        var q = new double[n, n];
[14837]416        ComputeSquaredEuclideanDistance(y, n, d, dd);
[14414]417
[14788]418        // Compute Q-matrix and normalization sum
419        var sumQ = double.Epsilon;
[14806]420        for(var n1 = 0; n1 < n; n1++) {
421          for(var m = 0; m < n; m++) {
422            if(n1 != m) {
[14788]423              q[n1, m] = 1 / (1 + dd[n1, m]);
424              sumQ += q[n1, m];
425            } else q[n1, m] = double.Epsilon;
426          }
427        }
[14806]428        for(var i = 0; i < n; i++) for(var j = 0; j < n; j++) q[i, j] /= sumQ;
[14414]429
[14788]430        // Sum t-SNE error
431        var c = .0;
[14806]432        for(var i = 0; i < n; i++)
433          for(var j = 0; j < n; j++) {
[14788]434            c += p[i, j] * Math.Log((p[i, j] + float.Epsilon) / (q[i, j] + float.Epsilon));
435          }
436        return c;
437      }
[14806]438
439      // The mapping of the approximate tSNE looks good but the error curve never changes.
[14788]440      private static double EvaluateErrorApproximate(IReadOnlyList<int> rowP, IReadOnlyList<int> colP, IReadOnlyList<double> valP, double[,] y, double theta) {
441        // Get estimate of normalization term
442        var n = y.GetLength(0);
443        var d = y.GetLength(1);
444        var tree = new SpacePartitioningTree(y);
445        var buff = new double[d];
446        double sumQ = 0.0;
[14806]447        for(var i = 0; i < n; i++) tree.ComputeNonEdgeForces(i, theta, buff, ref sumQ);
[14414]448
[14788]449        // Loop over all edges to compute t-SNE error
450        var c = .0;
[14806]451        for(var k = 0; k < n; k++) {
452          for(var i = rowP[k]; i < rowP[k + 1]; i++) {
[14788]453            var q = .0;
[14806]454            for(var j = 0; j < d; j++) buff[j] = y[k, j];
455            for(var j = 0; j < d; j++) buff[j] -= y[colP[i], j];
[14837]456            for(var j = 0; j < d; j++) q += buff[j] * buff[j];
457            q = (1.0 / (1.0 + q)) / sumQ;
[14788]458            c += valP[i] * Math.Log((valP[i] + float.Epsilon) / (q + float.Epsilon));
459          }
460        }
461        return c;
462      }
463      private static void SymmetrizeMatrix(IReadOnlyList<int> rowP, IReadOnlyList<int> colP, IReadOnlyList<double> valP, out int[] symRowP, out int[] symColP, out double[] symValP) {
[14414]464
[14788]465        // Count number of elements and row counts of symmetric matrix
466        var n = rowP.Count - 1;
467        var rowCounts = new int[n];
[14806]468        for(var j = 0; j < n; j++) {
469          for(var i = rowP[j]; i < rowP[j + 1]; i++) {
[14414]470
[14788]471            // Check whether element (col_P[i], n) is present
472            var present = false;
[14806]473            for(var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
474              if(colP[m] == j) present = true;
[14788]475            }
[14806]476            if(present) rowCounts[j]++;
[14788]477            else {
478              rowCounts[j]++;
479              rowCounts[colP[i]]++;
480            }
481          }
482        }
483        var noElem = 0;
[14806]484        for(var i = 0; i < n; i++) noElem += rowCounts[i];
[14414]485
[14788]486        // Allocate memory for symmetrized matrix
487        symRowP = new int[n + 1];
488        symColP = new int[noElem];
489        symValP = new double[noElem];
[14414]490
[14788]491        // Construct new row indices for symmetric matrix
492        symRowP[0] = 0;
[14806]493        for(var i = 0; i < n; i++) symRowP[i + 1] = symRowP[i] + rowCounts[i];
[14788]494
495        // Fill the result matrix
496        var offset = new int[n];
[14806]497        for(var j = 0; j < n; j++) {
498          for(var i = rowP[j]; i < rowP[j + 1]; i++) {                                  // considering element(n, colP[i])
[14788]499
500            // Check whether element (col_P[i], n) is present
501            var present = false;
[14806]502            for(var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
503              if(colP[m] != j) continue;
[14788]504              present = true;
[14806]505              if(j > colP[i]) continue; // make sure we do not add elements twice
[14788]506              symColP[symRowP[j] + offset[j]] = colP[i];
507              symColP[symRowP[colP[i]] + offset[colP[i]]] = j;
508              symValP[symRowP[j] + offset[j]] = valP[i] + valP[m];
509              symValP[symRowP[colP[i]] + offset[colP[i]]] = valP[i] + valP[m];
[14414]510            }
[14788]511
512            // If (colP[i], n) is not present, there is no addition involved
[14806]513            if(!present) {
[14788]514              symColP[symRowP[j] + offset[j]] = colP[i];
515              symColP[symRowP[colP[i]] + offset[colP[i]]] = j;
516              symValP[symRowP[j] + offset[j]] = valP[i];
517              symValP[symRowP[colP[i]] + offset[colP[i]]] = valP[i];
518            }
519
520            // Update offsets
[14806]521            if(present && (j > colP[i])) continue;
[14788]522            offset[j]++;
[14806]523            if(colP[i] != j) offset[colP[i]]++;
[14414]524          }
525        }
526
[14806]527        for(var i = 0; i < noElem; i++) symValP[i] /= 2.0;
[14414]528      }
[14807]529    }
[14788]530
[14807]531    /// <summary>
532    /// Simple interface to tSNE
533    /// </summary>
534    /// <param name="data"></param>
535    /// <param name="distance">The distance function used to differentiate similar from non-similar points, e.g. Euclidean distance.</param>
536    /// <param name="random">Random number generator</param>
537    /// <param name="newDimensions">Dimensionality of projected space (usually 2 for easy visual analysis).</param>
538    /// <param name="perplexity">Perplexity parameter of tSNE. Comparable to k in a k-nearest neighbour algorithm. Recommended value is floor(number of points /3) or lower</param>
539    /// <param name="iterations">Maximum number of iterations for gradient descent.</param>
540    /// <param name="theta">Value describing how much appoximated gradients my differ from exact gradients. Set to 0 for exact calculation and in [0,1] otherwise. CAUTION: exact calculation of forces requires building a non-sparse N*N matrix where N is the number of data points. This may exceed memory limitations.</param>
541    /// <param name="stopLyingIter">Number of iterations after which p is no longer approximated.</param>
542    /// <param name="momSwitchIter">Number of iterations after which the momentum in the gradient descent is switched.</param>
543    /// <param name="momentum">The initial momentum in the gradient descent.</param>
544    /// <param name="finalMomentum">The final momentum in gradient descent (after momentum switch).</param>
545    /// <param name="eta">Gradient descent learning rate.</param>
546    /// <returns></returns>
547    public static double[,] Run(T[] data, IDistance<T> distance, IRandom random,
548      int newDimensions = 2, double perplexity = 25, int iterations = 1000,
549      double theta = 0,
550      int stopLyingIter = 250, int momSwitchIter = 250, double momentum = .5,
551      double finalMomentum = .8, double eta = 200.0
552      ) {
553      var state = CreateState(data, distance, random, newDimensions, perplexity,
554        theta, stopLyingIter, momSwitchIter, momentum, finalMomentum, eta);
555
556      for(int i = 0; i < iterations - 1; i++) {
557        Iterate(state);
558      }
559      return Iterate(state);
[14414]560    }
[14785]561
[14807]562    public static TSNEState CreateState(T[] data, IDistance<T> distance, IRandom random,
563      int newDimensions = 2, double perplexity = 25, double theta = 0,
564      int stopLyingIter = 250, int momSwitchIter = 250, double momentum = .5,
565      double finalMomentum = .8, double eta = 200.0
[14788]566      ) {
567      return new TSNEState(data, distance, random, newDimensions, perplexity, theta, stopLyingIter, momSwitchIter, momentum, finalMomentum, eta);
568    }
[14414]569
570
[14788]571    public static double[,] Iterate(TSNEState state) {
[14806]572      if(state.exact)
[14788]573        ComputeExactGradient(state.p, state.newData, state.noDatapoints, state.newDimensions, state.dY);
574      else
575        ComputeApproximateGradient(state.rowP, state.colP, state.valP, state.newData, state.noDatapoints, state.newDimensions, state.dY, state.theta);
[14414]576
[14788]577      // Update gains
[14806]578      for(var i = 0; i < state.noDatapoints; i++) {
579        for(var j = 0; j < state.newDimensions; j++) {
[14788]580          state.gains[i, j] = Math.Sign(state.dY[i, j]) != Math.Sign(state.uY[i, j])
[14837]581            ? state.gains[i, j] + .2  // +0.2 nd *0.8 are used in two separate implementations of tSNE -> seems to be correct
582            : state.gains[i, j] * .8;
[14414]583
[14837]584          if(state.gains[i, j] < .01) state.gains[i, j] = .01;
[14414]585        }
[14788]586      }
[14414]587
[14788]588
589      // Perform gradient update (with momentum and gains)
[14806]590      for(var i = 0; i < state.noDatapoints; i++)
591        for(var j = 0; j < state.newDimensions; j++)
[14788]592          state.uY[i, j] = state.currentMomentum * state.uY[i, j] - state.eta * state.gains[i, j] * state.dY[i, j];
593
[14806]594      for(var i = 0; i < state.noDatapoints; i++)
595        for(var j = 0; j < state.newDimensions; j++)
[14788]596          state.newData[i, j] = state.newData[i, j] + state.uY[i, j];
597
598      // Make solution zero-mean
599      ZeroMean(state.newData);
[14807]600
[14788]601      // Stop lying about the P-values after a while, and switch momentum
[14806]602      if(state.iter == state.stopLyingIter) {
603        if(state.exact)
[14807]604          for(var i = 0; i < state.noDatapoints; i++)
605            for(var j = 0; j < state.noDatapoints; j++)
[14837]606              state.p[i, j] /= 12.0;
[14788]607        else
[14807]608          for(var i = 0; i < state.rowP[state.noDatapoints]; i++)
[14837]609            state.valP[i] /= 12.0;
[14414]610      }
[14788]611
[14806]612      if(state.iter == state.momSwitchIter)
[14788]613        state.currentMomentum = state.finalMomentum;
614
615      state.iter++;
616      return state.newData;
[14414]617    }
[14785]618
[14788]619
[14858]620
[14788]621    private static void ComputeApproximateGradient(int[] rowP, int[] colP, double[] valP, double[,] y, int n, int d, double[,] dC, double theta) {
622      var tree = new SpacePartitioningTree(y);
623      double sumQ = 0.0;
624      var posF = new double[n, d];
625      var negF = new double[n, d];
626      tree.ComputeEdgeForces(rowP, colP, valP, n, posF);
627      var row = new double[d];
[14806]628      for(var n1 = 0; n1 < n; n1++) {
[14858]629        Array.Clear(row,0,row.Length);
[14788]630        tree.ComputeNonEdgeForces(n1, theta, row, ref sumQ);
[14856]631        Buffer.BlockCopy(row, 0, negF, (sizeof(double) * n1 * d), d*sizeof(double));
[14788]632      }
633
634      // Compute final t-SNE gradient
[14856]635      for (var i = 0; i < n; i++)
[14806]636        for(var j = 0; j < d; j++) {
[14788]637          dC[i, j] = posF[i, j] - negF[i, j] / sumQ;
638        }
[14414]639    }
[14785]640
[14414]641    private static void ComputeExactGradient(double[,] p, double[,] y, int n, int d, double[,] dC) {
642
643      // Make sure the current gradient contains zeros
[14806]644      for(var i = 0; i < n; i++) for(var j = 0; j < d; j++) dC[i, j] = 0.0;
[14414]645
646      // Compute the squared Euclidean distance matrix
647      var dd = new double[n, n];
[14837]648      ComputeSquaredEuclideanDistance(y, n, d, dd);
[14414]649
650      // Compute Q-matrix and normalization sum
651      var q = new double[n, n];
652      var sumQ = .0;
[14806]653      for(var n1 = 0; n1 < n; n1++) {
654        for(var m = 0; m < n; m++) {
655          if(n1 == m) continue;
[14414]656          q[n1, m] = 1 / (1 + dd[n1, m]);
657          sumQ += q[n1, m];
658        }
659      }
660
661      // Perform the computation of the gradient
[14806]662      for(var n1 = 0; n1 < n; n1++) {
663        for(var m = 0; m < n; m++) {
664          if(n1 == m) continue;
[14414]665          var mult = (p[n1, m] - q[n1, m] / sumQ) * q[n1, m];
[14806]666          for(var d1 = 0; d1 < d; d1++) {
[14414]667            dC[n1, d1] += (y[n1, d1] - y[m, d1]) * mult;
668          }
669        }
670      }
671    }
[14788]672
[14414]673    private static void ComputeSquaredEuclideanDistance(double[,] x, int n, int d, double[,] dd) {
674      var dataSums = new double[n];
[14806]675      for(var i = 0; i < n; i++) {
676        for(var j = 0; j < d; j++) {
[14414]677          dataSums[i] += x[i, j] * x[i, j];
678        }
679      }
[14806]680      for(var i = 0; i < n; i++) {
681        for(var m = 0; m < n; m++) {
[14414]682          dd[i, m] = dataSums[i] + dataSums[m];
683        }
684      }
[14806]685      for(var i = 0; i < n; i++) {
[14414]686        dd[i, i] = 0.0;
[14806]687        for(var m = i + 1; m < n; m++) {
[14414]688          dd[i, m] = 0.0;
[14806]689          for(var j = 0; j < d; j++) {
[14414]690            dd[i, m] += (x[i, j] - x[m, j]) * (x[i, j] - x[m, j]);
691          }
692          dd[m, i] = dd[i, m];
693        }
694      }
695    }
696
697    private static void ZeroMean(double[,] x) {
698      // Compute data mean
699      var n = x.GetLength(0);
700      var d = x.GetLength(1);
701      var mean = new double[d];
[14806]702      for(var i = 0; i < n; i++) {
703        for(var j = 0; j < d; j++) {
[14414]704          mean[j] += x[i, j];
705        }
706      }
[14806]707      for(var i = 0; i < d; i++) {
[14414]708        mean[i] /= n;
709      }
710      // Subtract data mean
[14806]711      for(var i = 0; i < n; i++) {
712        for(var j = 0; j < d; j++) {
[14414]713          x[i, j] -= mean[j];
714        }
715      }
716    }
717  }
718}
Note: See TracBrowser for help on using the repository browser.