%TSNEM tSNE mapping
%
% W = TSNEM(A,K,N,P,MAX)
% W = A*TSNEM([],K,N,P,MAX)
% W = A*TSNEM(K,N,P,MAX)
% D = B*W
%
% INPUT
% A Dataset or matrix of doubles, used for training the mapping
% B Dataset, same dimensionality as A, to be mapped
% K Target dimension of mapping (default 2)
% N Initial dimension (default 30)
% P Perplexity (default 30)
% MAX Maximum number of iterations, default 1000
%
% OUTPUT
% W Trained mapping
% D 2D dataset
%
% DESCRIPTION
% This is PRTools inteface to the t-SNE Simple Matlab routine for high
% dimensional data visualisation. The output is a non-linear projection of
% the original vector space to a K-dimensional target space. The procedure
% starts with a preprocessing to N dimensions by PCA. The perplexity
% determines the number of neighbors taken into account, see references.
%
% The test dataset B is mapped on the target space by a linear mapping
% between the dissimilarity representation of A and the target space. See
% also multi-dimensional scaling by MDS or SAMMONM.
%
% EXAMPLE
% prdatasets; % make sure prdatasets is in the path
% a = satellite; % 36D dataset, 6 classes, 6435 objects
% [x,y] = gendat(a,0.5); % split in train and test set
% w = x*tsnem; % compute mapping
% figure; scattern(x*w); % show train set mapped to 2D: looks overtrained
% figure; scattern((x+randn(size(x))*1e-5)*w): % some noise helps
% figure; scattern(y*w); % show test set mapped to 2D
%
% LINK
% <a href="http://homepage.tudelft.nl/19j49/t-SNE.html">t-sne website</a>
%
% REFERENCES
% 1. L.J.P. van der Maaten and G.E. Hinton. Visualizing High-Dimensional
% Data Using t-SNE. J. of ML Research, 2579-2605, 2008.
% 2. L.J.P. van der Maaten. Learning a Parametric Embedding by Preserving
% Local Structure. Proc. 12th Int. Conf. on AI and Stats. (AI-STATS),
% JMLR W&CP 5:384-391, 2009.
% 3. L.J.P. van der Maaten. Barnes-Hut-SNE. Proc. Int. Conf. on Learning
% Representations.
% 4. E. Pekalska, D. de Ridder, R.P.W. Duin, and M.A. Kraaijveld, A new
% method of generalizing Sammon mapping with application to algorithm
% speed-up, ASCI99, Proc. 5th Annual ASCI Conf., 1999, 221-228. [<a href="http://rduin.nl/papers/asci_99_mds.pdf">pdf</a>]
%
% SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
% DATASETS, MAPPINGS, PCAM, MDS, SAMMONM, SCATTERD, SCATTERN
% Copyright: L. van der Maaten, R.P.W. Duin, r.p.w.duin@37steps.com
function out = tsnem(varargin)
argin = shiftargin(varargin,'scalar');
argin = setdefaults(argin,[],2,30,30,1000,0.01);
if mapping_task(argin,'definition')
% define untrained mapping
out = define_mapping(argin,'untrained');
elseif mapping_task(argin,'training')
% training the mapping
% retrieve input parameters
[a,n,indim,pers,max_iter,pow] = deal(argin{:});
aa = +a; % dataset construct is not needed
indim = min(indim,size(aa,2)); % we cannot go over input dim
% ready for calling tsne
y = tsne(aa,[],n,indim,pers,max_iter);
% refer to distances with a low power to get some smoothness
% compute linear transform between dissimilarity space and target
v = prpinv(sqrt(distm(aa)).^pow)*y;
% return output as a trained mapping
out = trained_mapping(a,{v,a,pow},2);
elseif mapping_task(argin,'trained execution')
% execution the mapping on new data
% retrieve inputs as given by PRTools
[b,w] = deal(argin{1:2}); % test data and mapping
[v,a,pow] = getdata(w); % transform, train data and power
% map and return
if isdataset(b)
% if input is dataset, return dataset
out = setdata(b,(sqrt(distm(b,a)).^pow)*v);
else
% otherwise return doubles
out = (sqrt(distm(b,a)).^pow)*v;
end
else
error('Illegal call');
end
return
function ydata = tsne(X, labels, no_dims, initial_dims, perplexity,max_iter)
%TSNE Performs symmetric t-SNE on dataset X
%
% mappedX = tsne(X, labels, no_dims, initial_dims, perplexity)
% mappedX = tsne(X, labels, initial_solution, perplexity)
%
% The function performs symmetric t-SNE on the NxD dataset X to reduce its
% dimensionality to no_dims dimensions (default = 2). The data is
% preprocessed using PCA, reducing the dimensionality to initial_dims
% dimensions (default = 30). Alternatively, an initial solution obtained
% from an other dimensionality reduction technique may be specified in
% initial_solution. The perplexity of the Gaussian kernel that is employed
% can be specified through perplexity (default = 30). The labels of the
% data are not used by t-SNE itself, however, they are used to color
% intermediate plots. Please provide an empty labels matrix [] if you
% don't want to plot results during the optimization.
% The low-dimensional data representation is returned in mappedX.
%
%
% (C) Laurens van der Maaten, 2010
% University of California, San Diego
if ~exist('labels', 'var')
labels = [];
end
if ~exist('no_dims', 'var') || isempty(no_dims)
no_dims = 2;
end
if ~exist('initial_dims', 'var') || isempty(initial_dims)
initial_dims = min(50, size(X, 2));
end
if ~exist('perplexity', 'var') || isempty(perplexity)
perplexity = 30;
end
% First check whether we already have an initial solution
if numel(no_dims) > 1
initial_solution = true;
ydata = no_dims;
no_dims = size(ydata, 2);
perplexity = initial_dims;
else
initial_solution = false;
end
% Normalize input data
X = X - min(X(:));
X = X / max(X(:));
X = bsxfun(@minus, X, mean(X, 1));
% Perform preprocessing using PCA
prwaitbaronce('Preprocessing %id data using PCA',min(size(X)));
if ~initial_solution
%disp('Preprocessing data using PCA...');
if size(X, 2) < size(X, 1)
C = X' * X;
else
C = (1 / size(X, 1)) * (X * X');
end
[M, lambda] = eig(C);
[lambda, ind] = sort(diag(lambda), 'descend');
M = M(:,ind(1:initial_dims));
lambda = lambda(1:initial_dims);
if ~(size(X, 2) < size(X, 1))
M = bsxfun(@times, X' * M, (1 ./ sqrt(size(X, 1) .* lambda))');
end
X = bsxfun(@minus, X, mean(X, 1)) * M;
clear M lambda ind
end
% Compute pairwise distance matrix
sum_X = sum(X .^ 2, 2);
D = bsxfun(@plus, sum_X, bsxfun(@plus, sum_X', -2 * (X * X')));
% Compute joint probabilities
P = d2p(D, perplexity, 1e-5); % compute affinities using fixed perplexity
clear D
prwaitbar(0)
% Run t-SNE
if initial_solution
ydata = tsne_p(P, labels, ydata,max_iter);
else
ydata = tsne_p(P, labels, no_dims,max_iter);
end
return
function ydata = tsne_p(P, labels, no_dims,max_iter)
%TSNE_P Performs symmetric t-SNE on affinity matrix P
%
% mappedX = tsne_p(P, labels, no_dims)
%
% The function performs symmetric t-SNE on pairwise similarity matrix P
% to create a low-dimensional map of no_dims dimensions (default = 2).
% The matrix P is assumed to be symmetric, sum up to 1, and have zeros
% on the diagonal.
% The labels of the data are not used by t-SNE itself, however, they
% are used to color intermediate plots. Please provide an empty labels
% matrix [] if you don't want to plot results during the optimization.
% The low-dimensional data representation is returned in mappedX.
%
%
% (C) Laurens van der Maaten, 2010
% University of California, San Diego
if ~exist('labels', 'var')
labels = [];
end
if ~exist('no_dims', 'var') || isempty(no_dims)
no_dims = 2;
end
% First check whether we already have an initial solution
if numel(no_dims) > 1
initial_solution = true;
ydata = no_dims;
no_dims = size(ydata, 2);
else
initial_solution = false;
end
% Initialize some variables
n = size(P, 1); % number of instances
momentum = 0.5; % initial momentum
final_momentum = 0.8; % value to which momentum is changed
mom_switch_iter = 250; % iteration at which momentum is changed
stop_lying_iter = 100; % iteration at which lying about P-values is stopped
%max_iter = 1000; % maximum number of iterations
epsilon = 500; % initial learning rate
min_gain = .01; % minimum gain for delta-bar-delta
% Make sure P-vals are set properly
P(1:n + 1:end) = 0; % set diagonal to zero
P = 0.5 * (P + P'); % symmetrize P-values
P = max(P ./ sum(P(:)), realmin); % make sure P-values sum to one
const = sum(P(:) .* log(P(:))); % constant in KL divergence
if ~initial_solution
P = P * 4; % lie about the P-vals to find better local minima
end
% Initialize the solution
if ~initial_solution
ydata = .0001 * randn(n, no_dims);
end
y_incs = zeros(size(ydata));
gains = ones(size(ydata));
% Run the iterations
tt = sprintf('Processing %i iterations: ',max_iter);
prwaitbar(max_iter,tt)
for iter=1:max_iter
% Compute joint probability that point i and j are neighbors
sum_ydata = sum(ydata .^ 2, 2);
num = 1 ./ (1 + bsxfun(@plus, sum_ydata, bsxfun(@plus, sum_ydata', -2 * (ydata * ydata')))); % Student-t distribution
num(1:n+1:end) = 0; % set diagonal to zero
Q = max(num ./ sum(num(:)), realmin); % normalize to get probabilities
% Compute the gradients (faster implementation)
L = (P - Q) .* num;
y_grads = 4 * (diag(sum(L, 1)) - L) * ydata;
% Update the solution
gains = (gains + .2) .* (sign(y_grads) ~= sign(y_incs)) ... % note that the y_grads are actually -y_grads
+ (gains * .8) .* (sign(y_grads) == sign(y_incs));
gains(gains < min_gain) = min_gain;
y_incs = momentum * y_incs - epsilon * (gains .* y_grads);
ydata = ydata + y_incs;
ydata = bsxfun(@minus, ydata, mean(ydata, 1));
% Update the momentum if necessary
if iter == mom_switch_iter
momentum = final_momentum;
end
if iter == stop_lying_iter && ~initial_solution
P = P ./ 4;
end
% Print out progress
if ~rem(iter, 10)
cost = const - sum(P(:) .* log(Q(:)));
prwaitbar(max_iter,iter,[tt num2str(iter) ', error: ' num2str(cost)]);
% cost = const - sum(P(:) .* log(Q(:)));
% disp(['Iteration ' num2str(iter) ': error is ' num2str(cost)]);
end
% Display scatter plot (maximally first three dimensions)
% if ~rem(iter, 10) && ~isempty(labels)
% if no_dims == 1
% scatter(ydata, ydata, 9, labels, 'filled');
% elseif no_dims == 2
% scatter(ydata(:,1), ydata(:,2), 9, labels, 'filled');
% else
% scatter3(ydata(:,1), ydata(:,2), ydata(:,3), 40, labels, 'filled');
% end
% axis tight
% axis off
% drawnow
% end
end
prwaitbar(0);
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [P, beta] = d2p(D, u, tol)
%D2P Identifies appropriate sigma's to get kk NNs up to some tolerance
%
% [P, beta] = d2p(D, kk, tol)
%
% Identifies the required precision (= 1 / variance^2) to obtain a Gaussian
% kernel with a certain uncertainty for every datapoint. The desired
% uncertainty can be specified through the perplexity u (default = 15). The
% desired perplexity is obtained up to some tolerance that can be specified
% by tol (default = 1e-4).
% The function returns the final Gaussian kernel in P, as well as the
% employed precisions per instance in beta.
%
%
% (C) Laurens van der Maaten, 2008
% Maastricht University
if ~exist('u', 'var') || isempty(u)
u = 15;
end
if ~exist('tol', 'var') || isempty(tol)
tol = 1e-4;
end
% Initialize some variables
n = size(D, 1); % number of instances
P = zeros(n, n); % empty probability matrix
beta = ones(n, 1); % empty precision vector
logU = log(u); % log of perplexity (= entropy)
% Run over all datapoints
for i=1:n
% if ~rem(i, 500)
% disp(['Computed P-values ' num2str(i) ' of ' num2str(n) ' datapoints...']);
% end
% Set minimum and maximum values for precision
betamin = -Inf;
betamax = Inf;
% Compute the Gaussian kernel and entropy for the current precision
[H, thisP] = Hbeta(D(i, [1:i - 1, i + 1:end]), beta(i));
% Evaluate whether the perplexity is within tolerance
Hdiff = H - logU;
tries = 0;
while abs(Hdiff) > tol && tries < 50
% If not, increase or decrease precision
if Hdiff > 0
betamin = beta(i);
if isinf(betamax)
beta(i) = beta(i) * 2;
else
beta(i) = (beta(i) + betamax) / 2;
end
else
betamax = beta(i);
if isinf(betamin)
beta(i) = beta(i) / 2;
else
beta(i) = (beta(i) + betamin) / 2;
end
end
% Recompute the values
[H, thisP] = Hbeta(D(i, [1:i - 1, i + 1:end]), beta(i));
Hdiff = H - logU;
tries = tries + 1;
end
% Set the final row of P
P(i, [1:i - 1, i + 1:end]) = thisP;
end
% disp(['Mean value of sigma: ' num2str(mean(sqrt(1 ./ beta)))]);
% disp(['Minimum value of sigma: ' num2str(min(sqrt(1 ./ beta)))]);
% disp(['Maximum value of sigma: ' num2str(max(sqrt(1 ./ beta)))]);
return
% Function that computes the Gaussian kernel values given a vector of
% squared Euclidean distances, and the precision of the Gaussian kernel.
% The function also computes the perplexity of the distribution.
function [H, P] = Hbeta(D, beta)
P = exp(-D * beta);
sumP = sum(P);
H = log(sumP) + beta * sum(D .* P) / sumP;
% why not: H = exp(-sum(P(P > 1e-5) .* log(P(P > 1e-5)))); ???
P = P / sumP;
return