Stochastic Neighbor Embedding

1. Knowl 1 — Probabilistic Objective for Stochastic Neighbor Embedding

   model/method

   Given a set of $n$ high-dimensional data vectors $\mathbf{x}_1, \dots, \mathbf{x}_n \in \mathbb{R}^D$ (or asymmetric pairwise dissimilarities $d_{ij}^2$), Stochastic Neighbor Embedding (SNE) models the neighborhood identity around each point probabilistically. The asymmetric probability $p_{j|i}$ that data point $\mathbf{x}_i$ picks data point $\mathbf{x}_j$ ($j \neq i$) as its neighbor is defined as:

   $p_{j|i} = \frac{\exp(-d_{ij}^2)}{\sum_{k \neq i} \exp(-d_{ik}^2)}, \quad p_{i|i} = 0$

   where the squared dissimilarity is typically computed as the scaled squared Euclidean distance $d_{ij}^2 = \frac{\|\mathbf{x}_i - \mathbf{x}_j\|^2}{2\sigma_i^2}$, with $\sigma_i^2$ being the variance of a Gaussian centered at $\mathbf{x}_i$.

   In the low-dimensional embedding space $\mathbb{R}^d$ ($d \ll D$), the induced probability $q_{j|i}$ that low-dimensional coordinate $\mathbf{y}_i$ picks $\mathbf{y}_j$ as its neighbor is defined using a Gaussian with fixed variance $\sigma^2 = \frac{1}{2}$:

   $q_{j|i} = \frac{\exp(-\|\mathbf{y}_i - \mathbf{y}_j\|^2)}{\sum_{k \neq i} \exp(-\|\mathbf{y}_i - \mathbf{y}_k\|^2)}, \quad q_{i|i} = 0$

   The optimal low-dimensional coordinates $\{\mathbf{y}_i\}_{i=1}^n$ are obtained by minimizing the sum of Kullback-Leibler divergences across all data points:

   $C = \sum_{i=1}^n \text{KL}(P_i \parallel Q_i) = \sum_{i=1}^n \sum_{j \neq i} p_{j|i} \log \frac{p_{j|i}}{q_{j|i}}$

   This objective penalizes modeling small high-dimensional distances with large low-dimensional distances ($p_{j|i}$ large, $q_{j|i}$ small), while placing a relatively smaller penalty on modeling large high-dimensional distances with small low-dimensional distances ($p_{j|i}$ small, $q_{j|i}$ large).

2. Knowl 2 — Gradient of the SNE Objective Function

   equation

   For the Stochastic Neighbor Embedding cost function $C = \sum_{i=1}^n \sum_{j \neq i} p_{j|i} \log \frac{p_{j|i}}{q_{j|i}}$, the analytic gradient with respect to the low-dimensional coordinate vector $\mathbf{y}_i \in \mathbb{R}^d$ is:

   $\frac{\partial C}{\partial \mathbf{y}_i} = 2 \sum_{j \neq i} (\mathbf{y}_i - \mathbf{y}_j)(p_{j|i} - q_{j|i} + p_{i|j} - q_{i|j})$

   where $p_{j|i}$ and $q_{j|i}$ are the high-dimensional and low-dimensional conditional neighborhood probabilities, respectively:

   $q_{j|i} = \frac{\exp(-\|\mathbf{y}_i - \mathbf{y}_j\|^2)}{\sum_{k \neq i} \exp(-\|\mathbf{y}_i - \mathbf{y}_k\|^2)}$

   The gradient represents a sum of forces pulling $\mathbf{y}_i$ toward $\mathbf{y}_j$ when $j$ is under-represented in the low-dimensional neighborhood ($p_{j|i} + p_{i|j} > q_{j|i} + q_{i|j}$) and pushing $\mathbf{y}_i$ away from $\mathbf{y}_j$ when $j$ is over-represented ($p_{j|i} + p_{i|j} < q_{j|i} + q_{i|j}$).

3. Knowl 3 — Gaussian Bandwidth Calibration via Perplexity

   model/method

   To handle variations in local density across the high-dimensional dataset, the variance $\sigma_i^2$ of the Gaussian centered at each data point $\mathbf{x}_i$ is calibrated individually. The variance is set such that the Shannon entropy of the resulting distribution over neighbors $P_i$ matches a target perplexity $k$:

   $H(P_i) = -\sum_{j \neq i} p_{j|i} \log_2 p_{j|i} = \log_2 k$

   where the conditional probabilities are given by:

   $p_{j|i} = \frac{\exp\left(-\frac{\|\mathbf{x}_i - \mathbf{x}_j\|^2}{2\sigma_i^2}\right)}{\sum_{l \neq i} \exp\left(-\frac{\|\mathbf{x}_i - \mathbf{x}_l\|^2}{2\sigma_i^2}\right)}$

   The hyperparameter $k$ is chosen by the user and represents the effective number of local neighbors. Because the entropy $H(P_i)$ is a strictly monotonic function of $\sigma_i$, the required value of $\sigma_i$ for each point is found efficiently via binary search.

4. Knowl 4 — Mixture Stochastic Neighbor Embedding for Multi-Valued Representations

   model/method

   To represent ambiguous objects (such as homonyms or manifold intersections) that map naturally to multiple distinct low-dimensional locations, SNE extends each high-dimensional object $i$ to possess $B$ low-dimensional image locations $\{\mathbf{y}_{i_b}\}_{b=1}^B$ with mixing proportions $\{\pi_{i_b}\}_{b=1}^B$ satisfying $\sum_{b=1}^B \pi_{i_b} = 1$ and $\pi_{i_b} \ge 0$.

   The low-dimensional neighborhood conditional probability $q_{j|i}$ is defined as a mixture across all image versions of object $i$ and all image versions of object $j$:

   $q_{j|i} = \sum_{b=1}^B \pi_{i_b} \sum_{c=1}^B \frac{\pi_{j_c} \exp(-\|\mathbf{y}_{i_b} - \mathbf{y}_{j_c}\|^2)}{\sum_{k \neq i} \sum_{d=1}^B \pi_{k_d} \exp(-\|\mathbf{y}_{i_b} - \mathbf{y}_{k_d}\|^2)}$

   Defining $r_{i_b j_c}$ as the probability that version $b$ of object $i$ picks version $c$ of object $j$ as its neighbor:

   $r_{i_b j_c} = \frac{\pi_{j_c} \exp(-\|\mathbf{y}_{i_b} - \mathbf{y}_{j_c}\|^2)}{\sum_{k \neq i} \sum_{d=1}^B \pi_{k_d} \exp(-\|\mathbf{y}_{i_b} - \mathbf{y}_{k_d}\|^2)}$

   the mixture probability can be written concisely as $q_{j|i} = \sum_{b=1}^B \pi_{i_b} \sum_{c=1}^B r_{i_b j_c}$.

5. Knowl 5 — Parameter Derivatives and Softmax Optimization for Mixture SNE

   equation

   In the mixture formulation of Stochastic Neighbor Embedding, the derivative of the cost function $C = \sum_i \sum_{j \neq i} p_{j|i} \log \frac{p_{j|i}}{q_{j|i}}$ with respect to the mixing proportion $\pi_{m_g}$ of version $g$ of object $m$ is given by:

   $\frac{\partial C}{\partial \pi_{m_g}} = -\sum_i \sum_{j \neq i} \frac{p_{j|i}}{q_{j|i}} \frac{\partial q_{j|i}}{\partial \pi_{m_g}}$

   where the derivative of the low-dimensional neighborhood distribution $q_{j|i}$ with respect to $\pi_{m_g}$ is:

   $\frac{\partial q_{j|i}}{\partial \pi_{m_g}} = \delta_{mi} \sum_c r_{m_g j_c} + \sum_b \frac{\pi_{i_b}}{\pi_{m_g}} r_{i_b m_g} \left[\delta_{mj} - \sum_c r_{i_b j_c}\right]$

   where $\delta_{ab} = 1$ if $a = b$ and $0$ otherwise, and $r_{i_b j_c} = \frac{\pi_{j_c} \exp(-\|\mathbf{y}_{i_b} - \mathbf{y}_{j_c}\|^2)}{\sum_{k \neq i} \sum_d \pi_{k_d} \exp(-\|\mathbf{y}_{i_b} - \mathbf{y}_{k_d}\|^2)}$.

   To ensure positivity and summation to 1 without constrained optimization, mixing proportions are parameterized via softmax weights $w_{i_b} \in \mathbb{R}$:

   $\pi_{i_b} = \frac{\exp(w_{i_b})}{\sum_{b'} \exp(w_{i_{b'}})}$

6. Knowl 6 — Stochastic Neighbor Embedding Optimization Procedure

   algorithm

   SNE computes low-dimensional coordinates $\mathbf{y}_i$ using batch gradient descent combined with decaying additive Gaussian jitter to escape poor local minima.

   Input: High-dimensional points {x_i}_{i=1}^n subset R^D, perplexity k, learning rate \eta, initial jitter std \sigma_jitter, jitter iterations T_jitter, total iterations T

   Output: Low-dimensional points {y_i}_{i=1}^n subset R^d

   for each point i = 1 to n:

       Find \sigma_i via binary search such that -\sum_{j != i} p_{j|i} log_2 p_{j|i} = log_2 k

       for each point j != i:

           p_{j|i} = exp(-||x_i - x_j||^2 / (2 \sigma_i^2)) / \sum_{l != i} exp(-||x_i - x_l||^2 / (2 \sigma_i^2))

       p_{i|i} = 0

   Initialize each y_i randomly in R^d close to the origin

   for iteration t = 1 to T:

       for each pair (i, j) with j != i:

           q_{j|i} = exp(-||y_i - y_j||^2) / \sum_{l != i} exp(-||y_i - y_l||^2)

       for each point i = 1 to n:

           g_i = 2 \sum_{j != i} (y_i - y_j) (p_{j|i} - q_{j|i} + p_{i|j} - q_{i|j})

           y_i = y_i - \eta * g_i

           if t <= T_jitter:

               Sample \epsilon_i ~ Normal(0, \sigma_jitter^2 I_d)

               y_i = y_i + \epsilon_i

   return {y_i}_{i=1}^n

   In standard USPS digit experiments ($n = 3000, d = 2$), the hyperparameters were $\eta = 0.2$, $\sigma_{\text{jitter}} = 0.3$ for $T_{\text{jitter}} = 3500$ iterations, followed by $\sigma_{\text{jitter}} = 0$ for an additional 500 iterations ($T = 4000$).

7. Knowl 7 — Large-Variance Limit of SNE as Metric Mismatch Minimization

   theoretical result

   When the high-dimensional Gaussian variance $\sigma_i^2 \to \infty$, the distribution over high-dimensional neighbors $P_i$ approaches a uniform distribution. In this high-variance regime, minimizing the SNE Kullback-Leibler cost function is mathematically equivalent to minimizing the squared distance mismatch between high- and low-dimensional representations, centered by their antigeometric means:

   $\text{Mismatch} = \sum_{i=1}^n \sum_{j=1}^n \left[(d_{ij}^2 - g_i^2) - (\tilde{d}_{ij}^2 - \tilde{g}_i^2)\right]^2$

   where $n$ is the total number of objects, and the scaled squared distances and antigeometric means are defined as:

   $d_{ij}^2 = \frac{\|\mathbf{x}_i - \mathbf{x}_j\|^2}{\sigma^2}, \quad g_i^2 = -\log \sum_{k \neq i} \frac{\exp(-d_{ik}^2)}{n-1}$

   $\tilde{d}_{ij}^2 = \frac{\|\mathbf{y}_i - \mathbf{y}_j\|^2}{\sigma^2}, \quad \tilde{g}_i^2 = -\log \sum_{k \neq i} \frac{\exp(-\tilde{d}_{ik}^2)}{n-1}$

   This equivalence links SNE in its large-variance limit to stress objective functions used in nonmetric multidimensional scaling (MDS).

8. Knowl 8 — SNE Dimension Reduction on USPS Handwritten Digit Bitmaps

   empirical result

   SNE was evaluated on 3000 grayscale $16 \times 16$ pixel digit bitmaps ($D = 256$) from the USPS database, comprising 600 examples from each of five classes (0, 1, 2, 3, 4), embedded into $d = 2$ dimensions with perplexity $k = 15$.

   Key empirical findings:

   • SNE separated the five digit classes into cleanly segregated clusters without access to class labels, while smoothly ordering continuous within-class geometric variations (orientation, skew, stroke thickness) across the 2D plane.
   • The SNE cost function converged to a value of $6719\text{ nats}$, compared to a theoretical baseline cost of $3000 \log_e(2999/15) \approx 15894\text{ nats}$ under a uniform low-dimensional neighbor distribution.
   • In comparison, Principal Component Analysis (PCA) projection onto the first two principal components failed to separate classes cleanly and jumbled boundaries between similar classes because PCA optimizes global variance preservation rather than local neighborhood preservation.

9. Knowl 9 — Mixture SNE Separation of Hybrid Ambiguous Objects and Topology Unfolding

   empirical result

   A simplified mixture SNE model was evaluated where each object had exactly two low-dimensional components ($B=2$) with equal mixing proportions $\pi_{i_1} = \pi_{i_2} = 0.5$, constrained by an attractive spring force pulling the two components together with a force increasing linearly up to a separation threshold of 0.05 (with a force of 0.025 nats per unit length at threshold) and remaining constant beyond it.

   Two experiments demonstrated its properties:

   1. Hybrid Digit Disambiguation: Tested on 300 USPS digits (100 each of classes 2, 3, 4) plus 100 synthetic hybrid digits generated by randomly sampling pixels from two different parent classes (perplexity $k = 10$, gradient multiplier 0.7, constant jitter 0.05). After optimization, $66\%$ of the hybrid digits had their two mixture components placed in widely separated locations corresponding precisely to the two parent digit clusters (and never in the third class cluster), compared to only $19\%$ of non-hybrid digits having separated components.
   2. Circle-to-Line Unfolding: Mixture SNE mapped a 2D circle onto a 1D line without altering neighbor relationships by splitting the point at the cut into two mixture components located at opposite ends of the line.

10. Knowl 10 — Author-Topic Mapping on NIPS Conference Text Corpus

    empirical result

    SNE was applied to embed 676 authors who published more than one paper in NIPS volumes 0--12 into $d = 2$ dimensions based on log aggregate author word counts.

    Experimental setup and results:

    • Representation: Author representations were generated from a vocabulary of 13,649 non-stop words appearing in $\ge 6$ papers, with co-authored counts allocated fractionally. Distances $d_{ij}$ were calculated as the Euclidean norm of differences between log aggregate word-count vectors.
    • Hyperparameter: Perplexity was set to $k = 25$ neighbors.
    • Outcome: SNE mapped authors into distinct spatial clusters that corresponded cleanly to broad research subfields, including generative models, support vector machines, neuroscience, reinforcement learning, and VLSI.

11. Knowl 11 — SNE Computational Complexity and Proposed Acceleration Strategies

    limitation

    The standard SNE gradient evaluation requires $O(n^2)$ operations per iteration due to all-pairs normalization in $q_{j|i}$ and force updates, requiring several hours of computation for $n = 3000$ points under gradient descent with noise annealing.

    The paper proposes three strategies to accelerate and improve optimization:

    1. Sparsification: Since high-dimensional probabilities $p_{j|i}$ are static, values below a threshold can be zeroed out and renormalized. Point pairs where both $p_{j|i} = 0$ and $p_{i|j} = 0$ can be ignored in gradient calculations when both $q_{j|i}$ and $q_{i|j}$ are also small.
    2. Spatial Data Structures: Identifying pairs with negligible $q_{j|i}$ can be performed in $O(\log n)$ time using geometric spatial trees (e.g., K-D trees, ball-trees, AD-trees) or fast multipole algorithms from $N$-body physics.
    3. Perplexity Annealing and Extra Dimensions: Gradually annealing the perplexity $k$ from large to small values, or optimizing in extra dimensions with an increasing penalty on non-zero coordinates, helps circumvent low-dimensional topological barriers without getting trapped in local optima.

Coverage note — The brief conceptual discussion interpreting SNE as a special case of Linear Relational Embedding (LRE) with an identity relation matrix was omitted as it is an interpretive analogy rather than a concrete contributed model or result.