当前位置：首页 > article >正文

Python实现t-分布随机邻域嵌入（t-SNE）降维算法

article 2024/9/22 4:14:53

- - Python实现t-分布随机邻域嵌入（t-SNE）降维算法的博客
  - - 引言
    - t-SNE算法原理
    - t-SNE的优势与局限
    - Python实现t-SNE算法
    - - 1. 创建t-SNE类
      - 2. 示例场景：MNIST手写数字数据集
      - 3. 结果分析
    - 结论
    - 运行结果

Python实现t-分布随机邻域嵌入（t-SNE）降维算法的博客

引言

在数据科学和机器学习中，降维技术是用于降低数据维度并且保留重要特征的关键方法。当我们处理高维数据时，常规的线性降维方法如PCA（主成分分析）可能不足以捕捉数据的非线性结构。为此，t-分布随机邻域嵌入（t-SNE）算法作为一种强大的非线性降维工具，被广泛用于高维数据的可视化和聚类分析。本文将详细介绍t-SNE算法的原理，并使用Python进行实现。我们将通过一个具体的场景来展示如何使用t-SNE算法实现降维。

t-SNE算法原理

t-SNE（t-Distributed Stochastic Neighbor Embedding）是一种非线性降维算法，旨在将高维数据映射到低维空间（通常是二维或三维），以便进行数据可视化。t-SNE通过保留高维空间中数据点的局部邻域结构，使得降维后的数据点在低维空间中保持相似的局部关系。t-SNE算法的核心思想可以分为以下几步：

高维空间中的相似度计算：
t-SNE首先在高维空间中计算数据点之间的相似度。具体而言，给定两个数据点 $x_i$ 和 $x_j$ ，其在高维空间中的相似度由条件概率 $p_{j|i}$ 表示。这个概率反映了在高维空间中选择点 $x_j$ 作为点 $x_i$ 邻居的概率。t-SNE假设高维数据分布遵循高斯分布，并通过高斯核函数计算条件概率：
$p_{j|i} = \frac{\exp(-\|x_i - x_j\|^2 / 2\sigma_i^2)}{\sum_{k \neq i} \exp(-\|x_i - x_k\|^2 / 2\sigma_i^2)}$
其中， $\sigma_i$ 是根据点 $x_i$ 的局部密度自适应调整的参数。
低维空间中的相似度计算：
在低维空间中，t-SNE使用t-分布而非高斯分布来计算相似度。t-分布在低维空间中具有较长的尾部，这使得t-SNE能够更好地处理数据点间的较大距离。给定两个低维数据点 $y_i$ 和 $y_j$ ，其相似度由条件概率 $q_{j|i}$ 表示：
$q_{ij} = \frac{(1 + \|y_i - y_j\|^2)^{-1}}{\sum_{k \neq l} (1 + \|y_k - y_l\|^2)^{-1}}$
t-分布通过其较长的尾部，更容易将远距离的数据点推远，从而使得局部结构更加清晰。
最小化KL散度：
t-SNE通过最小化高维空间中的条件概率分布 $P$ 和低维空间中的条件概率分布 $Q$ 之间的Kullback-Leibler散度（KL散度）来优化低维嵌入：
$\| Q) = \sum_i \sum_j p_{j|i} \log \frac{p_{j|i}}{q_{j|i}}$
通过最小化KL散度，t-SNE确保在低维空间中，具有相似关系的数据点保持紧密，而不相似的数据点分离开。
梯度下降优化：
为了最小化KL散度，t-SNE使用梯度下降法进行优化。算法通过逐步调整低维空间中的数据点位置，来减少KL散度的值，从而得到低维空间中的嵌入。

t-SNE的优势与局限

t-SNE的最大优势在于它能够很好地保留数据的局部结构，非常适合高维数据的可视化。它常用于揭示数据中的聚类或流形结构，尤其适用于处理复杂的非线性数据。

然而，t-SNE也存在一些局限性。首先，它的计算复杂度较高，尤其是在处理大规模数据时。其次，t-SNE对参数（如学习率和邻居数）的选择较为敏感，参数设置不当可能导致结果不理想。此外，t-SNE在降维时可能会丢失全局结构信息，这意味着它更适合用于局部结构的探索，而非全局结构的分析。

Python实现t-SNE算法

接下来，我们将使用Python实现t-SNE算法，并将其封装到一个面向对象的类中。

1. 创建t-SNE类

import numpy as np
from scipy.spatial.distance import pdist, squareform

class TSNE:
    def __init__(self, n_components=2, perplexity=30.0, learning_rate=200.0, n_iter=1000):
        self.n_components = n_components
        self.perplexity = perplexity
        self.learning_rate = learning_rate
        self.n_iter = n_iter
        self.Y = None  # Embeddings in low-dimensional space

    def _h_beta(self, D, beta):
        P = np.exp(-D * beta)
        sumP = np.sum(P)
        if sumP == 0:
            sumP = 1e-10  # Avoid division by zero
        H = np.log(sumP) + beta * np.sum(D * P) / sumP
        P = P / sumP
        return H, P

    def _binary_search(self, D, target, tol=1e-5, max_iter=50):
        beta_min = -np.inf
        beta_max = np.inf
        beta = 1.0
        H, P = self._h_beta(D, beta)
        Hdiff = H - np.log(target)
        iter_count = 0
        
        while np.abs(Hdiff) > tol and iter_count < max_iter:
            if Hdiff > 0:
                beta_min = beta
                if beta_max == np.inf:
                    beta *= 2
                else:
                    beta = (beta + beta_max) / 2
            else:
                beta_max = beta
                if beta_min == -np.inf:
                    beta /= 2
                else:
                    beta = (beta + beta_min) / 2
                    
            H, P = self._h_beta(D, beta)
            Hdiff = H - np.log(target)
            iter_count += 1
        
        return P

    def fit(self, X):
        n_samples = X.shape[0]
        target_perplexity = np.log(self.perplexity)
        P = np.zeros((n_samples, n_samples), dtype=np.float64)
        
        for i in range(n_samples):
            Di = np.square(X[i] - X).sum(axis=1).astype(np.float64)
            Di[i] = np.inf  # Exclude self-distance
            P[i] = self._binary_search(Di, target_perplexity)
        
        P = (P + P.T) / (2 * n_samples)
        P = np.maximum(P, 1e-12)  # Prevent log(0)
        
        self.Y = np.random.randn(n_samples, self.n_components)
        
        for iter in range(self.n_iter):
            D_Y = squareform(pdist(self.Y, 'sqeuclidean'))
            Q = np.exp(-D_Y / 2)
            Q[np.diag_indices_from(Q)] = 0  # Zero out diagonal (self-similarity)
            Q = np.maximum(Q, 1e-12)  # Prevent log(0)
            Q /= Q.sum()
            
            PQ = P - Q
            
            for i in range(n_samples):
                gradient = np.sum((PQ[:, i][:, np.newaxis] * (self.Y[i] - self.Y)), axis=0)
                self.Y[i] += self.learning_rate * gradient
            
            if iter % 100 == 0:
                cost = np.sum(P * np.log(P / Q))
                print(f"Iteration {iter}: cost = {cost}")

    def transform(self):
        """
        返回降维后的数据
        :return: 降维后的数据
        """
        return self.Y

2. 示例场景：MNIST手写数字数据集

为了演示如何使用t-SNE算法，我们将使用经典的MNIST手写数字数据集。该数据集包含10个类别的手写数字，每个样本是一个28x28的灰度图像。我们将使用t-SNE将784维的图像数据降维到2维空间，并进行可视化。

from sklearn.datasets import fetch_openml
import matplotlib.pyplot as plt

# 加载MNIST数据集
mnist = fetch_openml('mnist_784', version=1)
X, y = mnist["data"], mnist["target"]

# 随机抽取1000个样本进行降维
np.random.seed(42)
indices = np.random.choice(X.shape[0], 1000, replace=False)
X_subset, y_subset = X.iloc[indices], y.iloc[indices]

# 使用TSNE类进行降维
tsne = TSNE(n_components=2, perplexity=30.0, learning_rate=200.0, n_iter=1000)
tsne.fit(X_subset.values)
Y = tsne.transform()

# 可视化
plt.scatter(Y[:, 0], Y[:, 1], c=y_subset.astype(int), cmap='tab10', s=50, alpha=0.8)
plt.colorbar()
plt.title("t-SNE visualization of MNIST")
plt.show()

3. 结果分析

通过上述代码，我们可以将高维的MNIST数据集降维到2维，并对不同的数字类别进行可视化。我们可以观察到，t-SNE算法能够有效地将不同类别的数据点分离开来，形成不同的簇。每个簇代表一个手写数字类别，从而使我们能够直观地识别数据的聚类结构。

结论

t-SNE是一种强大的非线性降维技术，特别适用于高维数据的可视化和聚类分析。在本篇博客中，我们详细介绍了t-SNE算法的原理，并使用Python实现了该算法。通过对MNIST手写数字数据集的降维和可视化，我们展示了t-SNE在揭示数据结构方面的优势。然而，t-SNE在大规模数据集上的计算开销较高，对参数设置敏感，因此在实际应用中需要仔细调优和考虑其局限性。希望本文对您理解t-SNE算法有所帮助，并能在实际项目中应用这一技术。