source: branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs @ 14855

Last change on this file since 14855 was 14855, checked in by gkronber, 2 years ago

#2700: made some changes / bug-fixes while reviewing

File size: 28.4 KB
Line 
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;
59using HeuristicLab.Collections;
60using HeuristicLab.Common;
61using HeuristicLab.Core;
62using HeuristicLab.Optimization;
63using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
64using HeuristicLab.Random;
65
66namespace HeuristicLab.Algorithms.DataAnalysis {
67  [StorableClass]
68  public class TSNEStatic<T> {
69
70    [StorableClass]
71    public sealed class TSNEState : DeepCloneable {
72      // initialized once
73      [Storable]
74      public IDistance<T> distance;
75      [Storable]
76      public IRandom random;
77      [Storable]
78      public double perplexity;
79      [Storable]
80      public bool exact;
81      [Storable]
82      public int noDatapoints;
83      [Storable]
84      public double finalMomentum;
85      [Storable]
86      public int momSwitchIter;
87      [Storable]
88      public int stopLyingIter;
89      [Storable]
90      public double theta;
91      [Storable]
92      public double eta;
93      [Storable]
94      public int newDimensions;
95
96      // for approximate version: sparse representation of similarity/distance matrix
97      [Storable]
98      public double[] valP; // similarity/distance
99      [Storable]
100      public int[] rowP; // row index
101      [Storable]
102      public int[] colP; // col index
103
104      // for exact version: dense representation of distance/similarity matrix
105      [Storable]
106      public double[,] p;
107
108      // mapped data
109      [Storable]
110      public double[,] newData;
111
112      [Storable]
113      public int iter;
114      [Storable]
115      public double currentMomentum;
116
117      // helper variables (updated in each iteration)
118      [Storable]
119      public double[,] gains;
120      [Storable]
121      public double[,] uY;
122      [Storable]
123      public double[,] dY;
124
125      private TSNEState(TSNEState original, Cloner cloner) : base(original, cloner) {
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);
163      }
164
165      public override IDeepCloneable Clone(Cloner cloner) {
166        return new TSNEState(this, cloner);
167      }
168
169      [StorableConstructor]
170      public TSNEState(bool deserializing) { }
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;
182
183
184        // initialize
185        noDatapoints = data.Length;
186        if(noDatapoints - 1 < 3 * perplexity)
187          throw new ArgumentException("Perplexity too large for the number of data points!");
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];
194        for(var i = 0; i < noDatapoints; i++)
195          for(var j = 0; j < newDimensions; j++)
196            gains[i, j] = 1.0;
197
198        p = null;
199        rowP = null;
200        colP = null;
201        valP = null;
202
203        //Calculate Similarities
204        if(exact) p = CalculateExactSimilarites(data, distance, perplexity);
205        else CalculateApproximateSimilarities(data, distance, perplexity, out rowP, out colP, out valP);
206
207        // Lie about the P-values (factor is 4 in the MATLAB implementation)
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;
210
211        // Initialize solution (randomly)
212        var rand = new NormalDistributedRandom(random, 0, 1);
213        for(var i = 0; i < noDatapoints; i++)
214          for(var j = 0; j < newDimensions; j++)
215            newData[i, j] = rand.NextDouble() * .0001;
216      }
217
218      public double EvaluateError() {
219        return exact ?
220          EvaluateErrorExact(p, newData, noDatapoints, newDimensions) :
221          EvaluateErrorApproximate(rowP, colP, valP, newData, theta);
222      }
223
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
226        ComputeGaussianPerplexity(data, distance, out rowP, out colP, out valP, perplexity, (int)(3 * perplexity));
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;
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;
237      }
238
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
244        for(var n = 0; n < data.Length; n++) {
245          for(var m = n + 1; m < data.Length; m++) {
246            p[n, m] += p[m, n];
247            p[m, n] = p[n, m];
248          }
249        }
250        var sumP = .0;
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;
253        return p;
254      }
255
256      private static void ComputeGaussianPerplexity(IReadOnlyList<T> x, IDistance<T> distance, out int[] rowP, out int[] colP, out double[] valP, double perplexity, int k) {
257        if(perplexity > k) throw new ArgumentException("Perplexity should be lower than k!");
258
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;
266        for(var i = 0; i < n; i++) rowP[i + 1] = rowP[i] + k;
267
268        var objX = new List<IndexedItem<T>>();
269        for(var i = 0; i < n; i++) objX.Add(new IndexedItem<T>(i, x[i]));
270
271        // Build ball tree on data set
272        var tree = new VantagePointTree<IndexedItem<T>>(new IndexedItemDistance<T>(distance), objX);
273
274        // Loop over all points to find nearest neighbors
275        for(var i = 0; i < n; i++) {
276          IList<IndexedItem<T>> indices;
277          IList<double> distances;
278
279          // Find nearest neighbors
280          tree.Search(objX[i], k + 1, out indices, out distances);
281
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;
287          const double tol = 1e-5; 
288
289          // Iterate until we found a good perplexity
290          var iter = 0; double sumP = 0;
291          while(!found && iter < 200) { 
292
293            // Compute Gaussian kernel row
294            for(var m = 0; m < k; m++) curP[m] = Math.Exp(-beta * distances[m + 1]); 
295
296            // Compute entropy of current row
297            sumP = double.Epsilon;
298            for(var m = 0; m < k; m++) sumP += curP[m];
299            var h = .0;
300            for(var m = 0; m < k; m++) h += beta * (distances[m + 1] * curP[m]); 
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);
305            if(hdiff < tol && -hdiff < tol) {
306              found = true;
307            } else {
308              if(hdiff > 0) {
309                minBeta = beta;
310                if(maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
311                  beta *= 2.0;
312                else
313                  beta = (beta + maxBeta) / 2.0;
314              } else {
315                maxBeta = beta;
316                if(minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
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
328          for(var m = 0; m < k; m++) curP[m] /= sumP;
329          for(var m = 0; m < k; m++) {
330            colP[rowP[i] + m] = indices[m + 1].Index;
331            valP[rowP[i] + m] = curP[m];
332          }
333        }
334      }
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
341        for(var i = 0; i < n; i++) {
342          // Initialize some variables
343          var found = false;
344          var beta = 1.0;
345          var minBeta = double.MinValue;
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;
352          while(!found && iter < 200) {      // 200 iterations as in tSNE implementation by van der Maarten
353
354            // Compute Gaussian kernel row
355            for(var m = 0; m < n; m++) p[i, m] = Math.Exp(-beta * dd[i][m]);
356            p[i, i] = double.Epsilon;
357
358            // Compute entropy of current row
359            sumP = double.Epsilon;
360            for(var m = 0; m < n; m++) sumP += p[i, m];
361            var h = 0.0;
362            for(var m = 0; m < n; m++) h += beta * (dd[i][m] * p[i, m]);
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);
367            if(hdiff < tol && -hdiff < tol) {
368              found = true;
369            } else {
370              if(hdiff > 0) {
371                minBeta = beta;
372                if(maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
373                  beta *= 2.0;
374                else
375                  beta = (beta + maxBeta) / 2.0;
376              } else {
377                maxBeta = beta;
378                if(minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
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
390          for(var m = 0; m < n; m++) p[i, m] /= sumP;
391        }
392      }
393
394      private static double[][] ComputeDistances(T[] x, IDistance<T> distance) {
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();
410      }
411
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];
416        ComputeSquaredEuclideanDistance(y, n, d, dd);
417
418        // Compute Q-matrix and normalization sum
419        var sumQ = double.Epsilon;
420        for(var n1 = 0; n1 < n; n1++) {
421          for(var m = 0; m < n; m++) {
422            if(n1 != m) {
423              q[n1, m] = 1 / (1 + dd[n1, m]);
424              sumQ += q[n1, m];
425            } else q[n1, m] = double.Epsilon;
426          }
427        }
428        for(var i = 0; i < n; i++) for(var j = 0; j < n; j++) q[i, j] /= sumQ;
429
430        // Sum t-SNE error
431        var c = .0;
432        for(var i = 0; i < n; i++)
433          for(var j = 0; j < n; j++) {
434            c += p[i, j] * Math.Log((p[i, j] + float.Epsilon) / (q[i, j] + float.Epsilon));
435          }
436        return c;
437      }
438
439      // TODO: there seems to be a bug in the error approximation.
440      // The mapping of the approximate tSNE looks good but the error curve never changes.
441      private static double EvaluateErrorApproximate(IReadOnlyList<int> rowP, IReadOnlyList<int> colP, IReadOnlyList<double> valP, double[,] y, double theta) {
442        // Get estimate of normalization term
443        var n = y.GetLength(0);
444        var d = y.GetLength(1);
445        var tree = new SpacePartitioningTree(y);
446        var buff = new double[d];
447        double sumQ = 0.0;
448        for(var i = 0; i < n; i++) tree.ComputeNonEdgeForces(i, theta, buff, ref sumQ);
449
450        // Loop over all edges to compute t-SNE error
451        var c = .0;
452        for(var k = 0; k < n; k++) {
453          for(var i = rowP[k]; i < rowP[k + 1]; i++) {
454            var q = .0;
455            for(var j = 0; j < d; j++) buff[j] = y[k, j];
456            for(var j = 0; j < d; j++) buff[j] -= y[colP[i], j];
457            for(var j = 0; j < d; j++) q += buff[j] * buff[j];
458            q = (1.0 / (1.0 + q)) / sumQ;
459            c += valP[i] * Math.Log((valP[i] + float.Epsilon) / (q + float.Epsilon));
460          }
461        }
462        return c;
463      }
464      private static void SymmetrizeMatrix(IReadOnlyList<int> rowP, IReadOnlyList<int> colP, IReadOnlyList<double> valP, out int[] symRowP, out int[] symColP, out double[] symValP) {
465
466        // Count number of elements and row counts of symmetric matrix
467        var n = rowP.Count - 1;
468        var rowCounts = new int[n];
469        for(var j = 0; j < n; j++) {
470          for(var i = rowP[j]; i < rowP[j + 1]; i++) {
471
472            // Check whether element (col_P[i], n) is present
473            var present = false;
474            for(var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
475              if(colP[m] == j) present = true;
476            }
477            if(present) rowCounts[j]++;
478            else {
479              rowCounts[j]++;
480              rowCounts[colP[i]]++;
481            }
482          }
483        }
484        var noElem = 0;
485        for(var i = 0; i < n; i++) noElem += rowCounts[i];
486
487        // Allocate memory for symmetrized matrix
488        symRowP = new int[n + 1];
489        symColP = new int[noElem];
490        symValP = new double[noElem];
491
492        // Construct new row indices for symmetric matrix
493        symRowP[0] = 0;
494        for(var i = 0; i < n; i++) symRowP[i + 1] = symRowP[i] + rowCounts[i];
495
496        // Fill the result matrix
497        var offset = new int[n];
498        for(var j = 0; j < n; j++) {
499          for(var i = rowP[j]; i < rowP[j + 1]; i++) {                                  // considering element(n, colP[i])
500
501            // Check whether element (col_P[i], n) is present
502            var present = false;
503            for(var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
504              if(colP[m] != j) continue;
505              present = true;
506              if(j > colP[i]) continue; // make sure we do not add elements twice
507              symColP[symRowP[j] + offset[j]] = colP[i];
508              symColP[symRowP[colP[i]] + offset[colP[i]]] = j;
509              symValP[symRowP[j] + offset[j]] = valP[i] + valP[m];
510              symValP[symRowP[colP[i]] + offset[colP[i]]] = valP[i] + valP[m];
511            }
512
513            // If (colP[i], n) is not present, there is no addition involved
514            if(!present) {
515              symColP[symRowP[j] + offset[j]] = colP[i];
516              symColP[symRowP[colP[i]] + offset[colP[i]]] = j;
517              symValP[symRowP[j] + offset[j]] = valP[i];
518              symValP[symRowP[colP[i]] + offset[colP[i]]] = valP[i];
519            }
520
521            // Update offsets
522            if(present && (j > colP[i])) continue;
523            offset[j]++;
524            if(colP[i] != j) offset[colP[i]]++;
525          }
526        }
527
528        for(var i = 0; i < noElem; i++) symValP[i] /= 2.0;
529      }
530    }
531
532    /// <summary>
533    /// Simple interface to tSNE
534    /// </summary>
535    /// <param name="data"></param>
536    /// <param name="distance">The distance function used to differentiate similar from non-similar points, e.g. Euclidean distance.</param>
537    /// <param name="random">Random number generator</param>
538    /// <param name="newDimensions">Dimensionality of projected space (usually 2 for easy visual analysis).</param>
539    /// <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>
540    /// <param name="iterations">Maximum number of iterations for gradient descent.</param>
541    /// <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>
542    /// <param name="stopLyingIter">Number of iterations after which p is no longer approximated.</param>
543    /// <param name="momSwitchIter">Number of iterations after which the momentum in the gradient descent is switched.</param>
544    /// <param name="momentum">The initial momentum in the gradient descent.</param>
545    /// <param name="finalMomentum">The final momentum in gradient descent (after momentum switch).</param>
546    /// <param name="eta">Gradient descent learning rate.</param>
547    /// <returns></returns>
548    public static double[,] Run(T[] data, IDistance<T> distance, IRandom random,
549      int newDimensions = 2, double perplexity = 25, int iterations = 1000,
550      double theta = 0,
551      int stopLyingIter = 250, int momSwitchIter = 250, double momentum = .5,
552      double finalMomentum = .8, double eta = 200.0
553      ) {
554      var state = CreateState(data, distance, random, newDimensions, perplexity,
555        theta, stopLyingIter, momSwitchIter, momentum, finalMomentum, eta);
556
557      for(int i = 0; i < iterations - 1; i++) {
558        Iterate(state);
559      }
560      return Iterate(state);
561    }
562
563    public static TSNEState CreateState(T[] data, IDistance<T> distance, IRandom random,
564      int newDimensions = 2, double perplexity = 25, double theta = 0,
565      int stopLyingIter = 250, int momSwitchIter = 250, double momentum = .5,
566      double finalMomentum = .8, double eta = 200.0
567      ) {
568      return new TSNEState(data, distance, random, newDimensions, perplexity, theta, stopLyingIter, momSwitchIter, momentum, finalMomentum, eta);
569    }
570
571
572    public static double[,] Iterate(TSNEState state) {
573      if(state.exact)
574        ComputeExactGradient(state.p, state.newData, state.noDatapoints, state.newDimensions, state.dY);
575      else
576        ComputeApproximateGradient(state.rowP, state.colP, state.valP, state.newData, state.noDatapoints, state.newDimensions, state.dY, state.theta);
577
578      // Update gains
579      for(var i = 0; i < state.noDatapoints; i++) {
580        for(var j = 0; j < state.newDimensions; j++) {
581          state.gains[i, j] = Math.Sign(state.dY[i, j]) != Math.Sign(state.uY[i, j])
582            ? state.gains[i, j] + .2  // +0.2 nd *0.8 are used in two separate implementations of tSNE -> seems to be correct
583            : state.gains[i, j] * .8;
584
585          if(state.gains[i, j] < .01) state.gains[i, j] = .01;
586        }
587      }
588
589
590      // Perform gradient update (with momentum and gains)
591      for(var i = 0; i < state.noDatapoints; i++)
592        for(var j = 0; j < state.newDimensions; j++)
593          state.uY[i, j] = state.currentMomentum * state.uY[i, j] - state.eta * state.gains[i, j] * state.dY[i, j];
594
595      for(var i = 0; i < state.noDatapoints; i++)
596        for(var j = 0; j < state.newDimensions; j++)
597          state.newData[i, j] = state.newData[i, j] + state.uY[i, j];
598
599      // Make solution zero-mean
600      ZeroMean(state.newData);
601
602      // Stop lying about the P-values after a while, and switch momentum
603      if(state.iter == state.stopLyingIter) {
604        if(state.exact)
605          for(var i = 0; i < state.noDatapoints; i++)
606            for(var j = 0; j < state.noDatapoints; j++)
607              state.p[i, j] /= 12.0;
608        else
609          for(var i = 0; i < state.rowP[state.noDatapoints]; i++)
610            state.valP[i] /= 12.0;
611      }
612
613      if(state.iter == state.momSwitchIter)
614        state.currentMomentum = state.finalMomentum;
615
616      state.iter++;
617      return state.newData;
618    }
619
620
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];
628      for(var n1 = 0; n1 < n; n1++) {
629        Buffer.BlockCopy(negF, (sizeof(double) * n1 * d), row, 0, d);
630        tree.ComputeNonEdgeForces(n1, theta, row, ref sumQ);
631      }
632
633      // Compute final t-SNE gradient
634      for(var i = 0; i < n; i++)
635        for(var j = 0; j < d; j++) {
636          dC[i, j] = posF[i, j] - negF[i, j] / sumQ;
637        }
638    }
639
640    private static void ComputeExactGradient(double[,] p, double[,] y, int n, int d, double[,] dC) {
641
642      // Make sure the current gradient contains zeros
643      for(var i = 0; i < n; i++) for(var j = 0; j < d; j++) dC[i, j] = 0.0;
644
645      // Compute the squared Euclidean distance matrix
646      var dd = new double[n, n];
647      ComputeSquaredEuclideanDistance(y, n, d, dd);
648
649      // Compute Q-matrix and normalization sum
650      var q = new double[n, n];
651      var sumQ = .0;
652      for(var n1 = 0; n1 < n; n1++) {
653        for(var m = 0; m < n; m++) {
654          if(n1 == m) continue;
655          q[n1, m] = 1 / (1 + dd[n1, m]);
656          sumQ += q[n1, m];
657        }
658      }
659
660      // Perform the computation of the gradient
661      for(var n1 = 0; n1 < n; n1++) {
662        for(var m = 0; m < n; m++) {
663          if(n1 == m) continue;
664          var mult = (p[n1, m] - q[n1, m] / sumQ) * q[n1, m];
665          for(var d1 = 0; d1 < d; d1++) {
666            dC[n1, d1] += (y[n1, d1] - y[m, d1]) * mult;
667          }
668        }
669      }
670    }
671
672    private static void ComputeSquaredEuclideanDistance(double[,] x, int n, int d, double[,] dd) {
673      var dataSums = new double[n];
674      for(var i = 0; i < n; i++) {
675        for(var j = 0; j < d; j++) {
676          dataSums[i] += x[i, j] * x[i, j];
677        }
678      }
679      for(var i = 0; i < n; i++) {
680        for(var m = 0; m < n; m++) {
681          dd[i, m] = dataSums[i] + dataSums[m];
682        }
683      }
684      for(var i = 0; i < n; i++) {
685        dd[i, i] = 0.0;
686        for(var m = i + 1; m < n; m++) {
687          dd[i, m] = 0.0;
688          for(var j = 0; j < d; j++) {
689            dd[i, m] += (x[i, j] - x[m, j]) * (x[i, j] - x[m, j]);
690          }
691          dd[m, i] = dd[i, m];
692        }
693      }
694    }
695
696    private static void ZeroMean(double[,] x) {
697      // Compute data mean
698      var n = x.GetLength(0);
699      var d = x.GetLength(1);
700      var mean = new double[d];
701      for(var i = 0; i < n; i++) {
702        for(var j = 0; j < d; j++) {
703          mean[j] += x[i, j];
704        }
705      }
706      for(var i = 0; i < d; i++) {
707        mean[i] /= n;
708      }
709      // Subtract data mean
710      for(var i = 0; i < n; i++) {
711        for(var j = 0; j < d; j++) {
712          x[i, j] -= mean[j];
713        }
714      }
715    }
716  }
717}
Note: See TracBrowser for help on using the repository browser.