t-Distributed Elastic Embedding

This article retains its original experiments and examples. For current package usage, see the usage guide.

Elastic Embedding (PDF) is a method like SNE, but doesn’t use normalization of weights into probabilities. I have looked at its behavior elsewhere. In this discussion, I would like to consider an obvious extension: EE uses Gaussian output weights, like ASNE and SSNE, which we know don’t work as well in the normalized case compared to using the t-distribution with one degree of freedom, as used in t-SNE. Can EE also be improved by using the t-distribution as an output weight function? I was unable to find any example of anyone trying this, which was slightly surprising to me. If someone has tried this, it’s at least not a very popular or well-known technique.

Theory

Cost function

The generic elastic embedding cost function is:

C_{\mathrm{EE}} = -\sum_{ij}^{N} v_{ij}^{+} \log w_{ij} + \lambda \sum_{ij}^{N} v_{ij}^{-} w_{ij}

where N is the number of observations in the dataset, i and j are indexes to observations in the dataset, v_{ij}^{+} is the input weight (similarity) between i and j associated with the attractive part of the cost function, v_{ij}^{-} are input weights associated with the repulsive part of the cost function, w_{ij} are the output weights and \lambda is a positive scalar that weights the repulsive versus attractive interactions in the cost function.

For more details, see the introduction to the discussion of elastic embedding, but the attractive input weights are the output weights are Gaussian-based, as in ASNE and SSNE.

Here, I propose a simplified cost function for t-distributed Elastic Embedding (t-EE)

C_{\text{t-EE}} = -\frac{1}{N}\sum_{ij}^{N} v_{ij} \log w_{ij} + \frac{\lambda}{N} \sum_{ij}^{N} w_{ij}

The t-EE cost function differs from EE by:

only having one set of input weights, v_{ij}, which are equivalent to the unnormalized versions of the probabilities p_{ij} used in SSNE and t-SNE and determined by the same procedure of perplexity calibration, followed by row-normalization and symmetrization.
w_{ij} are the t-distributed output weights (same as in t-SNE).
Both parts of the cost function are weighted by 1 / N, the reason for which we will see below.

EE as the Weighted I-Divergence

It would be nice if the cost function was zero when the output weights perfectly reproduce the input weights, as with the Kullback-Leibler divergence used by t-SNE. As long as we don’t add any terms involving w_{ij}, we can make this happen without affecting the EE gradient. We would like C_{\text{t-EE}} = 0 when v_{ij} = w_{ij}. Let’s add an extra constant, C_{v}, that only involves v_{ij} (and maybe \lambda), and we can write:

\begin{aligned} C_{\text{t-EE}} &= -\frac{1}{N}\sum_{ij}^{N} v_{ij} \log v_{ij} + \frac{\lambda}{N} \sum_{ij}^{N} v_{ij} + C_{v} \\ &= 0 \end{aligned} and then, knowing that \sum_{ij}^{N} v_{ij} = N, we have:

C_{v} = \frac{1}{N}\sum_{ij}^{N} v_{ij} \log v_{ij} - \lambda

so the cost function can be written as:

C_{\text{t-EE}} = \frac{1}{N}\sum_{ij}^{N} v_{ij} \log v_{ij} -\frac{1}{N}\sum_{ij}^{N} v_{ij} \log w_{ij} + \frac{\lambda}{N} \sum_{ij}^{N} w_{ij} - \lambda The first two terms can be combined to give a Kullback-Leibler-like expression, and we can group the \lambda expressions too (and I have also rewritten N as \sum_{ij}^{N} v_{ij})

C_{\text{t-EE}} = \frac{1}{N} \left[ \sum_{ij}^{N} v_{ij} \log \frac{v_{ij}}{w_{ij}} + \lambda \left(\sum_{ij}^{N} w_{ij} - \sum_{ij}^{N} v_{ij}\right) \right]

This expression is basically the same as the I-divergence (also called the generalized Kullback-Leibler divergence), used in non-negative matrix factorization, but with a weighting term introduced by \lambda. The introduction of \lambda does allow the cost to get negative, but empirically, I’ve noticed that a negative cost indicates that \lambda has been set too small: visually the output shows over-compressed clusters compared to a t-SNE result. So this version of the cost function may be superior for helping to at least set a lower bound on \lambda.

The figures below use that cost. In current smallvis, the comparison is with the scaled output weights \lambda w_{ij}. This adds \lambda - 1 - \log \lambda to the expression above and gives the usual nonnegative I-divergence with the same gradient.

Gradient

Whichever version of the cost function you prefer the gradient of t-EE with respect to the output coordinate of point i, \mathbf{y}_i, is:

\frac{\partial C_{\text{t-EE}}}{\partial \mathbf{y}_i} = \frac{4}{N} \sum_{j}^N \left(v_{ij} - \lambda w_{ij} \right) w_{ij} \left(\mathbf{y}_i - \mathbf{y}_j\right)

If we were to define \lambda at each iteration as N / Z where Z =\sum_{ij}^{N} w_{ij}, and then set p_{ij} = v_{ij} / N and q_{ij} = w_{ij} / Z, we would get back the t-SNE gradient. For elastic embedding, we fix \lambda as constant. The connection between the t-EE and the t-SNE gradient is the reason for making the weighting factor in the cost function \lambda / N.

Datasets

See the Datasets page. We shall be using these datasets not just to benchmark the results of t-EE, but also to think about how t-SNE weight distributions change over time.

Choosing lambda

Like EE, t-EE requires choosing the value of an extra parameter, \lambda, specified as lambda in the smallvis interface. This is certainly an irksome requirement. At least with the way t-EE is formulated we can think about \lambda as an approximation to N/Z in t-SNE, i.e. the ratio of the size of the dataset to the sum of the output weights. Because the output weights range from 0 to 1, we can try to get a handle on the sort of range \lambda will have.

If using the random initialization, at the start of t-SNE optimization, all output distances are small, which means all the weights are large and close to 1. As there are N observations, there are N^2 elements in \mathbf{W} so Z \approx N^2. Technically, it’s actually more like N^2 - N because self-weights are always set to 0, but as N grows the difference becomes less important. So we shall say that at the start of a t-SNE optimization, \lambda = N / Z \approx N / N^2 = 1 / N.

If the t-SNE optimization did a perfect job and the error was reduced to zero, the sum of the input and output weights would be the same and Z = N, which would mean that \lambda = 1. Over the course of t-SNE optimization, then, we would expect \lambda to start very small and increase, approaching 1 as the error got close to 0.

So we can at least expect \lambda to very roughly take values between 0 and 1. How good are these bounds in reality? Well, they’re ok. I have looked at how N/Z evolves in t-SNE and N / Z = 1 / N is a very good approximation even if you are initializing from Laplacian Eigenmaps or scaled PCA. The value quickly rise, but eventually levels off.

And of course, we never see an error of zero with t-SNE, so the final value of N / Z is not 1. For most datasets and perplexities, it’s substantially smaller than 1. The exact value depends on N and the perplexity, with N / Z showing an inverse relationship to both: large N and large perplexity results in a lower final value of N / Z. For low perplexities (e.g. 5), N / Z can actually end up being larger than 1, although I’ve not seen it larger than around 5, although this might be because low perplexity requires more iterations for convergence to occur, and I might just not have been patient.

Evolution of N/Z with iteration

Here’s a little extra detail on how N/Z changes during a t-SNE optimization.

Swiss Roll

A synthetic dataset like the Swiss Roll let’s use look at the effect of dataset size in a controlled way. We’ll look at N = 1000 and N = 3000 and several perplexity values. Some relevant parameters used are:

scale = FALSE, Y_init = "lap", eta = 100, max_iter = 5000, epoch = 100, tol = 0

First, here are some lots for N/Z evolves with iteration and different values of perplexity for sr1k on the left and sr3k on the right. Black line is perplexity = 10, red is perplexity = 20, orange is perplexity = 30, blue is perplexity = 40 and purple is perplexity = 50.

As perplexity increases, the final value of N/Z decreases.

Next, some plots of how N/Z evolves for perplexity 10 and perplexity 50 with different N. Black is for sr1k, red is sr2k, blue is sr3k.

This plot shows that as N increases, the final value of N/Z decreases.

Real Datasets

Here are some plots of how 1/Z evolves with iteration for each dataset. For each dataset, the Laplacian Eigenmap initialization is shown in black, scaled PCA is shown in red, random initialization is shown in blue and random initialization with early exaggeration is in cyan.

Results seem to not depend on initialization very much. The most obvious difference occurs with iris at perplexity = 100, but the final value of N/Z is still the same. In general, a steady increase in N/Z is observed with an elbow in every plot when during the transition of low to high momentum at mom_switch_iter = 250.

As with the Swiss Roll results, higher perplexities lead to lower values of N/Z, and datasets with larger N tend to end up with lower values of N/Z. iris bucks this trend however.

With the standard settings (e.g. perplexity = 40 and max_iter = 1000), here is a summary of the N / Z values starting from scaled PCA results:

dataset	N	N / Z
iris	150	0.057
s1k	1000	0.20
oli	400	0.33
frey	1965	0.21
coil20	1440	0.20
mnist	6000	0.12
fashion	6000	0.088

For typical t-SNE settings, we see that N/Z attains values between 0.057 to 0.33. Ignoring iris for a moment, the smallest dataset, oli, has the largest N/Z = 0.33, while mnist6k and fashion6k have N/Z = 0.12 and N/Z = 0.09, respectively. Intermediate-sized datasets s1k, frey and coil20 with N between 1000 and 2000, have intermediate values of N/Z at around 0.2. Unfortunately, iris messes this all up, having the smallest N/Z value despite being the smallest dataset.

Converged Results

4 May 2020 From the plots above, it’s clear that N/Z values haven’t leveled off by 1000 iterations for larger datasets like mnist6k and fashion6k, and no dataset has converged for the perplexity = 5 results. I ran some very long runs with max_iter = 100000, eta = 100, epoch = 100, exaggeration_factor = 4, Y_init = "spca" with perplexity = 40 and perplexity = 100 to ensure that the final results were converged with the expectation that this would be way more than enough iterations. For perplexity = 5, even 100,000 iterations isn’t enough to fully converge the results: the layout doesn’t really change, it just slowly expands with the cost slowly ticking down every few thousand iterations. I assume that this is caused by the low perplexity failing to capture any longer-range connected structure, so you can always decrease the cost by pushing all the mini clusters away from each other forever, so that non-neighbors get less and less probability assigned to them.

As an alternative, for a low perplexity case I used perplexity = 15 and only max_iter = 10000. This still shows similar issues as perplexity = 5, and doesn’t converge with tol = 1e-7, but runs in a more reasonable time-frame.

dataset	N	N / Z (15)	N / Z (40)	N / Z (100)
iris	150	2.40	0.057	0.016
s1k	1000	2.14	0.20	0.035
oli	400	1.69	0.33	0.075
frey	1965	0.89	0.32	0.094
coil20	1440	1.08	0.31	0.067
mnist6k	6000	0.68	0.60	0.19
fashion6k	6000	0.63	0.29	0.065

Compared to the table above, frey and coil20 were slightly under-converged, but mnist6k and fashion6k change a lot and clearly had a lot more expansion to do. Interestingly, mnist6k and fashion6k have quite similar values of N/Z at low perplexity, but show very different behavior at higher perplexity. While clearly N/Z descreases with increased perplexity for all datasets, the exact behavior seems to be very dataset dependent. For the t-EE plots below we’ll show results for multiple \lambda values, with the results here used to establish a sensible set of trial values.

Settings

I was able to use the standard delta-bar-delta settings with t-EE, without having to reduce eta, the learning rate. Therefore we will use the same settings, I tend to use for t-SNE, based off those used in the original t-SNE paper, except with scaled PCA initialization rather than random initialization and without early exaggeration. As well as looking at perplexity = 40, we’ll also repeat results using perplexity = 5 and perplexity = 100 to look at how t-EE performs across a practical range of perplexity values.

Example command for embedding iris with perplexity = 5:

iris_tee_p5 <- smallvis(
  iris,
  scale = "absmax",
  perplexity = 5,
  Y_init = "spca",
  eta = 100,
  method = list("tee", lambda = 0.05),
  min_cost = -Inf
)

I used min_cost = -Inf here because the original cost could go below zero.

Values for lambda will be varied between 0.001 and 0.5. On the plots below, both lambda and the t-EE error (I-divergence style) will be shown. t-SNE results (with the corresponding KL divergence) with the same perplexity and other settings are also shown for comparison.

Results: Perplexity 5

iris

s1k

oli

frey

coil20

mnist

fashion

Results: Perplexity 40

iris

s1k

oli