第22章 无监督学习方法总结 - 实现降维可视化

import numpy as np
import matplotlib.pyplot as plt

def tsne(X, n_components=2, perplexity=30.0, n_iter=1000, learning_rate=200.0):
    n_samples = X.shape[0]
    
    distances = np.zeros((n_samples, n_samples))
    for i in range(n_samples):
        for j in range(i+1, n_samples):
            distances[i, j] = distances[j, i] = np.linalg.norm(X[i] - X[j])
    
    P = np.zeros((n_samples, n_samples))
    for i in range(n_samples):
        beta = 1.0
        for _ in range(50):
            exp_dist = np.exp(-distances[i] * beta)
            exp_dist[i] = 0
            sum_exp = exp_dist.sum()
            if sum_exp == 0:
                break
            H = -np.sum(exp_dist / sum_exp * np.log(exp_dist / sum_exp + 1e-10))
            if H > np.log(perplexity):
                beta *= 2
            else:
                beta /= 2
        P[i] = exp_dist / (2 * sum_exp)
    
    P = (P + P.T) / (2 * n_samples)
    
    Y = np.random.randn(n_samples, n_components) * 1e-4
    velocity = np.zeros_like(Y)
    
    for t in range(n_iter):
        Q = np.zeros((n_samples, n_samples))
        for i in range(n_samples):
            for j in range(i+1, n_samples):
                Q[i, j] = Q[j, i] = 1 / (1 + np.linalg.norm(Y[i] - Y[j])**2)
        
        Q = Q / Q.sum()
        
        grad = np.zeros_like(Y)
        for i in range(n_samples):
            for j in range(n_samples):
                if i != j:
                    grad[i] += 4 * (P[i, j] - Q[i, j]) * (Y[i] - Y[j]) * (1 + np.linalg.norm(Y[i] - Y[j])**2)**(-1)
        
        momentum = 0.5 if t < 250 else 0.8
        velocity = momentum * velocity - learning_rate * grad
        Y = Y + velocity
    
    return Y

def plot_tsne(X, labels=None, title='t-SNE Visualization'):
    Y = tsne(X)
    
    plt.figure(figsize=(10, 8))
    if labels is not None:
        scatter = plt.scatter(Y[:, 0], Y[:, 1], c=labels, cmap='tab10', alpha=0.6)
        plt.colorbar(scatter)
    else:
        plt.scatter(Y[:, 0], Y[:, 1], alpha=0.6)
    
    plt.title(title)
    plt.xlabel('Component 1')
    plt.ylabel('Component 2')
    plt.show()