<jupyter_start><jupyter_text>t-SNE structure in linear data by <NAME>Recently a tweet https://twitter.com/mikelove/status/738021869341839360 shook data scientists awake, showing that one of the most used dimensionality reduction techniques stochastic neighberhood embedding is 'making up' structure in a linear embedding. Their plots show incredible detailed structure in a high dimensional dataset containing only a one dimensional linear embedded input. Here we explore the same problem using python's implementations of tSNE and show that GPLVM and Bayesian GPLVM do not suffer from this problem.<jupyter_code>n = 200
m = 40
np.random.seed(1)
x = np.random.uniform(-1, 1, n)
c = np.digitize(x, np.linspace(-1,1,12))-1
cols = np.asarray(sns.color_palette('spectral_r',12))[c]<jupyter_output><empty_output><jupyter_text>Linear embeddingThe embedded linear structure. It is just a point cloud on a line betwenn -1 and 1.<jupyter_code>plt.scatter(x, np.zeros(n), c=cols, cmap='spectral_r', lw=0)<jupyter_output><empty_output><jupyter_text>Embed into higher dimensionsWe embedd this linear point cloud in a 40 dimensional space using a 40 dimensional set of 40 orthogonal bases.<jupyter_code>r = np.random.normal(0,1,[m,1])
M = np.eye(m)-2*(r.dot(r.T)/r.T.dot(r))
X = np.c_[x, np.zeros((n,m-1))].dot(M)
plt.scatter(*X[:,1:3].T, c=x, cmap='spectral', lw=0)<jupyter_output><empty_output><jupyter_text>PCAThe algorithm PCA finds the linear point cloud and is able to reproduce the original.<jupyter_code>from sklearn.decomposition import PCA
Xpca = PCA(2).fit_transform(X)
plt.scatter(*Xpca.T, c=x, cmap='spectral', lw=0)<jupyter_output><empty_output><jupyter_text>Barnes_hut implementationThis implementation of tSNE is taken from https://github.com/danielfrg/tsne.Note the 'artificial' structure in the embeddings. It is not as detailed ('bad') as in the original R implementation of tSNE, but still it introduces some structure and non-linearities.For other implementations in R have a look at https://gist.github.com/mikelove/74bbf5c41010ae1dc94281cface90d32, where others looked at different implementations of tSNE in R.<jupyter_code>fig, axes = plt.subplots(2,3,tight_layout=True,figsize=(15,10))
axit = axes.flat
for lr in range(6):
    ax = next(axit)
    Xtsne = tsne.bh_sne(X.copy())
    ax.scatter(*Xtsne.T, c=cols, cmap='spectral', lw=0)
    ax.set_title('restart: ${}$'.format(lr))<jupyter_output>/home/maxz/anaconda/lib/python2.7/site-packages/matplotlib/figure.py:1744: UserWarning:This figure includes Axes that are not compatible with tight_layout, so its results might be incorrect.<jupyter_text>scikit-learn emplementationCareful with the scikit-learn implementation, the parameters and convergence seems to be unstable.A too high learning rate can lead to a completely random picture. Note that a learning rate of 1000. is the standard in the library.It seems that the scikit-learn implementation suffers from high dependence on the parameters chosen for learning. Beware of this, when using this implementation of tSNE.<jupyter_code>from sklearn.manifold import TSNE
fig, axes = plt.subplots(2,3,tight_layout=True,figsize=(15,10))
axit = axes.flat

for perplexity in [30,60]:
    for lr in [200, 500, 1000]:
        ax = next(axit)
        Xtsne = TSNE(perplexity=perplexity, learning_rate=lr, init='pca').fit_transform(X.copy())
        ax.scatter(*Xtsne.T, c=cols, cmap='spectral', lw=0)
        ax.set_title('perp=${}$, learn-rate=${}$'.format(perplexity, lr))<jupyter_output><empty_output><jupyter_text>GPLVMWe show the performance of GPLVM on the dataset provided. The top right bar plot shows the sensitivity to the dimensions of the embedding. When the sensitivity is high, the dimension is used.<jupyter_code>from mpl_toolkits.axes_grid.inset_locator import inset_axes

fig, axes = plt.subplots(2,3,tight_layout=True,figsize=(15,10))
axit = axes.flat
for lr in range(6):
    ax = next(axit)
    m = GPy.models.GPLVM(X.copy(), 2)
    m.optimize(messages=1, gtol=0, clear_after_finish=True)
    msi = m.get_most_significant_input_dimensions()[:2]
    
    is_ = m.kern.input_sensitivity().copy()
    is_ /= is_.max()
    
    XBGPLVM = m.X[:,msi] * is_[np.array(msi)]
    #m.kern.plot_ARD(ax=ax)
    ax.scatter(*XBGPLVM.T, c=cols, cmap='spectral', lw=0)
    ax.set_title('restart: ${}$'.format(lr))
    ax.set_xlabel('dimension ${}$'.format(msi[0]))
    ax.set_ylabel('dimension ${}$'.format(msi[1]))
    
    a = inset_axes(ax,
                    width="30%", # width = 30% of parent_bbox
                    height='20%', # height : 1 inch
                    loc=1)
    sns.barplot(np.array(msi), is_[np.array(msi)], label='input-sens', ax=a)
    a.set_title('sensitivity')
    a.set_xlabel('dimension')<jupyter_output>/home/maxz/anaconda/lib/python2.7/site-packages/matplotlib/__init__.py:892: UserWarning:axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.<jupyter_text>Bayesian GPLVMLast, we want to show the ability of Bayesian GPLVM to not only find the linear imbedding, but also determine the number of necessary dimensions to embedd the data in. We start at a dimensionality of 5 and show that Bayesian GPLVM is capable of identifying one dimension as the significant one and switches off the other ones. See for that the significance of the dimensions in the top right of each subplot.<jupyter_code>from mpl_toolkits.axes_grid.inset_locator import inset_axes


fig, axes = plt.subplots(2,3,tight_layout=True,figsize=(15,10))
axit = axes.flat
for lr in range(6):
    ax = next(axit)
    m = GPy.models.BayesianGPLVM(X, 5, num_inducing=25)
    m.optimize(messages=1, gtol=0, clear_after_finish=True)
    msi = m.get_most_significant_input_dimensions()[:2]
    
    is_ = m.kern.input_sensitivity()
    is_ /= is_.max()
    
    XBGPLVM = m.X.mean[:,msi] * is_[np.array(msi)]
    #m.kern.plot_ARD(ax=ax)
    ax.scatter(*XBGPLVM.T, c=cols, cmap='spectral', lw=0)
    ax.set_title('restart: ${}$'.format(lr))
    ax.set_xlabel('dimension ${}$'.format(msi[0]))
    ax.set_ylabel('dimension ${}$'.format(msi[1]))
    
    a = inset_axes(ax,
                    width="30%", # width = 30% of parent_bbox
                    height='20%', # height : 1 inch
                    loc=1)
    sns.barplot(range(m.input_dim), is_, label='input-sens')
    a.set_title('sensitivity')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1218401632909405, 'fraction_numerical': 0.025278693672476055, 'mean_word_length': 4.602462620932278, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 32, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '365586', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2334, 'n_tokens_neox': 2166, 'n_words': 579, 'path': 'compbio/tSNE-over-interpretation.ipynb', 'repo_name': 'olumighty1/notebook', 'size': '6376'}