t-SNE快速降维算法详解与实现

t-SNE（t-distributed Stochastic Neighbor Embedding）是一种非常有效的非线性降维技术，特别适用于高维数据的可视化。

t-SNE算法核心原理

算法概述

t-SNE通过保留数据点之间的相似性将高维数据映射到低维空间（通常是2D或3D），使得在高维空间中相似的点在低维空间中靠近，不相似的点远离。

MATLAB实现

classdef FastTSNE
    % 快速t-SNE降维实现
    % 支持多种优化技巧和加速策略
    
    properties
        % 算法参数
        num_dimensions    % 输出维度
        perplexity       % 困惑度
        learning_rate    % 学习率
        max_iter         % 最大迭代次数
        momentum         % 动量
        early_exaggeration % 早期放大因子
        random_state     % 随机种子
        
        % 内部变量
        P                % 高维相似度矩阵
        Y                % 低维嵌入
        gains            % 增益矩阵
        uY               % 更新方向
    end
    
    methods
        function obj = FastTSNE(num_dimensions, perplexity, learning_rate, max_iter)
            % 构造函数
            if nargin < 4
                max_iter = 1000;
            end
            if nargin < 3
                learning_rate = 200;
            end
            if nargin < 2
                perplexity = 30;
            end
            if nargin < 1
                num_dimensions = 2;
            end
            
            obj.num_dimensions = num_dimensions;
            obj.perplexity = perplexity;
            obj.learning_rate = learning_rate;
            obj.max_iter = max_iter;
            obj.momentum = 0.8;
            obj.early_exaggeration = 4;
            obj.random_state = 42;
        end
        
        function [Y, costs] = fit_transform(obj, X, varargin)
            % 主函数：执行t-SNE降维
            % X: 输入数据 (n_samples x n_features)
            
            fprintf('开始t-SNE降维...\n');
            fprintf('数据维度: %d × %d → %dD\n', size(X,1), size(X,2), obj.num_dimensions);
            
            % 参数解析
            p = inputParser;
            addParameter(p, 'verbose', true, @islogical);
            addParameter(p, 'plot_progress', false, @islogical);
            parse(p, varargin{:});
            
            % 设置随机种子
            rng(obj.random_state);
            
            % 数据标准化
            X = obj.standardize_data(X);
            
            % 计算高维相似度
            fprintf('计算高维相似度...\n');
            obj.P = obj.compute_high_dimensional_similarities(X);
            
            % 初始化低维嵌入
            obj.Y = obj.initialize_embedding(size(X, 1));
            obj.gains = ones(size(obj.Y));
            obj.uY = zeros(size(obj.Y));
            
            % 早期放大
            P = obj.P * obj.early_exaggeration;
            
            % 优化过程
            costs = zeros(obj.max_iter, 1);
            
            for iter = 1:obj.max_iter
                % 计算低维相似度和梯度
                [Q, grad] = obj.compute_gradient(P, obj.Y);
                
                % 更新嵌入
                obj = obj.update_embedding(grad, iter);
                
                % 计算代价
                costs(iter) = obj.compute_cost(P, Q);
                
                % 早期放大阶段结束
                if iter == 100
                    P = obj.P;
                end
                
                % 进度显示
                if p.Results.verbose && mod(iter, 100) == 0
                    fprintf('迭代 %d/%d, 代价: %.4f\n', iter, obj.max_iter, costs(iter));
                end
                
                % 绘制进度
                if p.Results.plot_progress && mod(iter, 50) == 0
                    obj.plot_current_embedding(iter, costs);
                end
            end
            
            Y = obj.Y;
            fprintf('t-SNE完成！最终代价: %.4f\n', costs(end));
        end
        
        function X_std = standardize_data(~, X)
            % 数据标准化
            X_std = X - mean(X, 1);
            X_std = X_std ./ std(X_std, 0, 1);
            X_std(isnan(X_std)) = 0; % 处理常数特征
        end
        
        function P = compute_high_dimensional_similarities(obj, X)
            % 计算高维空间中的相似度概率
            
            n = size(X, 1);
            P = zeros(n);
            
            % 计算成对距离
            D = pdist2(X, X).^2;
            
            % 为每个点寻找合适的方差（困惑度）
            for i = 1:n
                % 二分搜索寻找合适的方差
                beta = obj.binary_search_perplexity(D(i, :), obj.perplexity);
                
                % 计算条件概率
                P(i, :) = exp(-D(i, :) * beta);
                P(i, i) = 0; % 设置自相似度为0
                
                % 归一化
                P(i, :) = P(i, :) / sum(P(i, :));
            end
            
            % 对称化概率矩阵
            P = (P + P') / (2 * n);
            
            % 防止数值问题
            P = max(P, 1e-12);
        end
        
        function beta = binary_search_perplexity(obj, D_row, target_perplexity)
            % 二分搜索寻找合适的方差以达到目标困惑度
            
            max_iter = 50;
            tol = 1e-4;
            
            beta = 1;
            beta_min = -inf;
            beta_max = inf;
            
            P = exp(-D_row * beta);
            P(D_row == 0) = 0; % 排除自身
            sum_P = sum(P);
            if sum_P == 0
                sum_P = 1e-12;
            end
            P = P / sum_P;
            
            H = obj.compute_entropy(P);
            perplexity_diff = H - log(target_perplexity);
            
            iter = 0;
            while abs(perplexity_diff) > tol && iter < max_iter
                if perplexity_diff > 0
                    beta_min = beta;
                    if isinf(beta_max)
                        beta = beta * 2;
                    else
                        beta = (beta + beta_max) / 2;
                    end
                else
                    beta_max = beta;
                    if isinf(beta_min)
                        beta = beta / 2;
                    else
                        beta = (beta + beta_min) / 2;
                    end
                end
                
                P = exp(-D_row * beta);
                P(D_row == 0) = 0;
                sum_P = sum(P);
                if sum_P == 0
                    sum_P = 1e-12;
                end
                P = P / sum_P;
                
                H = obj.compute_entropy(P);
                perplexity_diff = H - log(target_perplexity);
                iter = iter + 1;
            end
        end
        
        function H = compute_entropy(~, P)
            % 计算概率分布的香农熵
            P(P == 0) = 1e-12; % 避免log(0)
            H = -sum(P .* log(P));
        end
        
        function Y = initialize_embedding(obj, n_samples)
            % 初始化低维嵌入
            Y = 1e-4 * randn(n_samples, obj.num_dimensions);
        end
        
        function [Q, grad] = compute_gradient(obj, P, Y)
            % 计算梯度和低维相似度
            
            n = size(Y, 1);
            
            % 计算低维相似度 (使用t分布)
            sum_Y = sum(Y .^ 2, 2);
            D = sum_Y + sum_Y' - 2 * (Y * Y');
            Q = 1 ./ (1 + D);
            Q(1:n+1:end) = 0; % 对角线置零
            Q = Q / sum(Q(:));
            Q = max(Q, 1e-12);
            
            % 计算梯度
            PQ = P - Q;
            grad = zeros(size(Y));
            
            for i = 1:n
                % 向量化梯度计算
                diff = Y(i, :) - Y;
                grad(i, :) = 4 * sum(bsxfun(@times, PQ(:, i), diff), 1);
            end
        end
        
        function obj = update_embedding(obj, grad, iter)
            % 更新低维嵌入
            
            % 调整增益
            obj.gains = (obj.gains + 0.2) .* (sign(grad) ~= sign(obj.uY)) + ...
                       (obj.gains * 0.8) .* (sign(grad) == sign(obj.uY));
            obj.gains = max(obj.gains, 0.01);
            
            % 更新方向
            obj.uY = obj.momentum * obj.uY - obj.learning_rate * obj.gains .* grad;
            
            % 更新位置
            obj.Y = obj.Y + obj.uY;
            
            % 中心化
            obj.Y = obj.Y - mean(obj.Y, 1);
        end
        
        function cost = compute_cost(~, P, Q)
            % 计算KL散度代价
            cost = sum(P(:) .* log(P(:) ./ Q(:)));
        end
        
        function plot_current_embedding(obj, iter, costs)
            % 绘制当前嵌入状态
            figure(1);
            subplot(1, 2, 1);
            scatter(obj.Y(:, 1), obj.Y(:, 2), 30, 'filled');
            title(sprintf('t-SNE嵌入 (迭代 %d)', iter));
            xlabel('维度 1'); ylabel('维度 2');
            axis equal;
            grid on;
            
            subplot(1, 2, 2);
            plot(1:iter, costs(1:iter), 'b-', 'LineWidth', 2);
            title('代价函数收敛');
            xlabel('迭代次数'); ylabel('KL散度');
            grid on;
            
            drawnow;
        end
    end
    
    methods (Static)
        function demo_tsne()
            % t-SNE演示函数
            
            fprintf('t-SNE降维算法演示\n');
            fprintf('================\n\n');
            
            % 生成测试数据
            [X, labels] = FastTSNE.generate_test_data();
            
            % 创建t-SNE对象
            tsne = FastTSNE(2, 30, 200, 1000);
            
            % 执行降维
            tic;
            [Y, costs] = tsne.fit_transform(X, 'verbose', true, 'plot_progress', true);
            time_elapsed = toc;
            
            fprintf('计算时间: %.2f 秒\n', time_elapsed);
            
            % 绘制最终结果
            FastTSNE.plot_results(X, Y, labels, costs);
            
            % 比较不同参数
            FastTSNE.compare_parameters(X, labels);
        end
        
        function [X, labels] = generate_test_data()
            % 生成测试数据
            
            rng(42);
            n_points = 300;
            
            % 生成多个高斯分布的数据
            centers = [0, 0; 5, 5; -5, 5; 5, -5; -5, -5];
            n_clusters = size(centers, 1);
            points_per_cluster = ceil(n_points / n_clusters);
            
            X = [];
            labels = [];
            
            for i = 1:n_clusters
                cluster_points = centers(i, :) + randn(points_per_cluster, 2) * 0.8;
                X = [X; cluster_points];
                labels = [labels; i * ones(points_per_cluster, 1)];
            end
            
            X = X(1:n_points, :);
            labels = labels(1:n_points);
            
            % 添加噪声特征，使其成为高维数据
            X_high_dim = [X, randn(n_points, 8) * 0.5];
            
            fprintf('生成测试数据: %d个样本, %d维 -> %d维\n', ...
                size(X_high_dim, 1), size(X_high_dim, 2), 2);
        end
        
        function plot_results(X, Y, labels, costs)
            % 绘制结果对比
            
            figure('Position', [100, 100, 1500, 500]);
            
            % 原始数据（前两个维度）
            subplot(1, 3, 1);
            gscatter(X(:, 1), X(:, 2), labels);
            title('原始数据 (前两个维度)');
            xlabel('特征 1'); ylabel('特征 2');
            axis equal; grid on;
            
            % t-SNE结果
            subplot(1, 3, 2);
            gscatter(Y(:, 1), Y(:, 2), labels);
            title('t-SNE降维结果');
            xlabel('t-SNE 维度 1'); ylabel('t-SNE 维度 2');
            axis equal; grid on;
            
            % 收敛曲线
            subplot(1, 3, 3);
            plot(1:length(costs), costs, 'b-', 'LineWidth', 2);
            title('t-SNE收敛曲线');
            xlabel('迭代次数'); ylabel('KL散度代价');
            grid on;
            
            sgtitle('t-SNE降维效果展示');
        end
        
        function compare_parameters(X, labels)
            % 比较不同参数的效果
            
            fprintf('\n参数比较分析\n');
            fprintf('============\n');
            
            perplexities = [5, 30, 50, 100];
            learning_rates = [10, 100, 200, 500];
            
            figure('Position', [100, 100, 1200, 800]);
            
            % 比较不同困惑度
            subplot(2, 2, 1);
            for i = 1:length(perplexities)
                tsne = FastTSNE(2, perplexities(i), 200, 500);
                Y = tsne.fit_transform(X, 'verbose', false);
                scatter(Y(:, 1), Y(:, 2), 20, labels, 'filled');
                title(sprintf('困惑度 = %d', perplexities(i)));
                axis equal; axis off;
            end
            
            % 比较不同学习率
            subplot(2, 2, 2);
            for i = 1:length(learning_rates)
                tsne = FastTSNE(2, 30, learning_rates(i), 500);
                Y = tsne.fit_transform(X, 'verbose', false);
                scatter(Y(:, 1), Y(:, 2), 20, labels, 'filled');
                title(sprintf('学习率 = %d', learning_rates(i)));
                axis equal; axis off;
            end
            
            % 比较不同迭代次数
            subplot(2, 2, 3);
            iterations = [100, 250, 500, 1000];
            for i = 1:length(iterations)
                tsne = FastTSNE(2, 30, 200, iterations(i));
                Y = tsne.fit_transform(X, 'verbose', false);
                scatter(Y(:, 1), Y(:, 2), 20, labels, 'filled');
                title(sprintf('迭代次数 = %d', iterations(i)));
                axis equal; axis off;
            end
            
            % 3D t-SNE示例
            subplot(2, 2, 4);
            tsne_3d = FastTSNE(3, 30, 200, 500);
            Y_3d = tsne_3d.fit_transform(X, 'verbose', false);
            scatter3(Y_3d(:, 1), Y_3d(:, 2), Y_3d(:, 3), 30, labels, 'filled');
            title('3D t-SNE');
            axis equal; grid on;
            view(45, 30);
            
            sgtitle('t-SNE参数比较');
        end
    end
end

% 快速t-SNE实现（Barnes-Hut近似）
classdef FastTSNE_BH < FastTSNE
    % 基于Barnes-Hut近似的快速t-SNE实现
    % 适用于大规模数据集
    
    methods
        function obj = FastTSNE_BH(num_dimensions, perplexity, learning_rate, max_iter)
            % 构造函数
            obj = obj@FastTSNE(num_dimensions, perplexity, learning_rate, max_iter);
            obj.momentum = 0.5;
        end
        
        function [Y, costs] = fit_transform(obj, X, varargin)
            % Barnes-Hut快速t-SNE实现
            
            fprintf('使用Barnes-Hut快速t-SNE...\n');
            
            % 参数解析
            p = inputParser;
            addParameter(p, 'theta', 0.5, @(x) x > 0 && x < 1); % Barnes-Hut参数
            addParameter(p, 'verbose', true, @islogical);
            parse(p, varargin{:});
            
            % 数据标准化
            X = obj.standardize_data(X);
            
            % 计算高维相似度（使用更高效的方法）
            obj.P = obj.compute_similarities_bh(X);
            
            % 初始化
            obj.Y = obj.initialize_embedding(size(X, 1));
            costs = zeros(obj.max_iter, 1);
            
            % 优化过程
            for iter = 1:obj.max_iter
                % 使用Barnes-Hut计算梯度和代价
                [costs(iter), grad] = obj.barnes_hut_gradient(obj.P, obj.Y, p.Results.theta);
                
                % 更新嵌入
                obj.Y = obj.Y + obj.learning_rate * grad;
                obj.Y = obj.Y - mean(obj.Y, 1);
                
                % 调整学习率
                if iter > 100
                    obj.learning_rate = obj.learning_rate * 0.99;
                end
                
                if p.Results.verbose && mod(iter, 100) == 0
                    fprintf('迭代 %d/%d, 代价: %.4f\n', iter, obj.max_iter, costs(iter));
                end
            end
            
            Y = obj.Y;
        end
        
        function P = compute_similarities_bh(obj, X)
            % 使用更高效的方法计算高维相似度
            
            n = size(X, 1);
            
            % 使用k近邻近似
            k = min(3 * obj.perplexity, n - 1);
            
            % 计算k近邻
            [~, D] = knnsearch(X, X, 'K', k+1); % 包含自身
            D = D(:, 2:end).^2; % 排除自身
            
            P = zeros(n, n);
            
            for i = 1:n
                % 只考虑k近邻
                beta = obj.binary_search_perplexity([zeros(1, i-1), inf, D(i, :)], obj.perplexity);
                P(i, :) = exp(-pdist2(X(i, :), X).^2 * beta);
                P(i, i) = 0;
                P(i, :) = P(i, :) / sum(P(i, :));
            end
            
            P = (P + P') / (2 * n);
            P = max(P, 1e-12);
        end
        
        function [cost, grad] = barnes_hut_gradient(~, P, Y, theta)
            % Barnes-Hut近似计算梯度
            
            n = size(Y, 1);
            grad = zeros(size(Y));
            cost = 0;
            
            % 构建四叉树（2D）或八叉树（3D）
            tree = obj.build_tree(Y);
            
            % 计算梯度
            for i = 1:n
                [grad_i, cost_i] = obj.compute_gradient_for_point(P, Y, i, tree, theta);
                grad(i, :) = grad_i;
                cost = cost + cost_i;
            end
            
            cost = cost / n;
        end
        
        function tree = build_tree(~, Y)
            % 构建空间分割树（简化实现）
            % 实际实现应该更复杂，这里使用简化版本
            
            tree.points = Y;
            tree.center = mean(Y, 1);
            tree.size = max(range(Y, 1));
            tree.num_points = size(Y, 1);
        end
        
        function [grad, cost] = compute_gradient_for_point(~, P, Y, i, tree, theta)
            % 为单个点计算梯度（简化实现）
            
            n = size(Y, 1);
            grad = zeros(1, size(Y, 2));
            cost = 0;
            
            % 计算低维相似度
            dist_sq = sum((Y(i, :) - Y).^2, 2);
            Q = 1 ./ (1 + dist_sq);
            Q(i) = 0;
            Q = Q / sum(Q);
            
            % 计算梯度和代价
            for j = 1:n
                if j ~= i
                    diff = Y(i, :) - Y(j, :);
                    mult = (P(i, j) - Q(j)) * Q(j);
                    grad = grad + 4 * mult * diff;
                    
                    if P(i, j) > 0
                        cost = cost + P(i, j) * log(P(i, j) / Q(j));
                    end
                end
            end
        end
    end
end

% 主测试函数
function main_tsne_demo()
    fprintf('t-SNE快速降维算法完整演示\n');
    fprintf('=========================\n\n');
    
    % 选择演示模式
    fprintf('选择演示模式:\n');
    fprintf('1 - 基础t-SNE演示\n');
    fprintf('2 - 大规模数据测试\n');
    fprintf('3 - 与PCA比较\n');
    
    choice = input('请输入选择 (1-3): ');
    
    switch choice
        case 1
            FastTSNE.demo_tsne();
        case 2
            demo_large_scale();
        case 3
            demo_comparison_pca();
        otherwise
            fprintf('无效选择，运行基础演示...\n');
            FastTSNE.demo_tsne();
    end
end

function demo_large_scale()
    % 大规模数据演示
    
    fprintf('\n大规模数据t-SNE演示\n');
    fprintf('==================\n');
    
    % 生成更大规模的数据
    n_samples = 2000;
    n_features = 50;
    
    fprintf('生成大规模数据: %d样本 × %d特征\n', n_samples, n_features);
    
    X = randn(n_samples, n_features);
    
    % 添加聚类结构
    n_clusters = 5;
    cluster_size = floor(n_samples / n_clusters);
    labels = zeros(n_samples, 1);
    
    for i = 1:n_clusters
        start_idx = (i-1) * cluster_size + 1;
        end_idx = min(i * cluster_size, n_samples);
        cluster_points = start_idx:end_idx;
        
        % 为每个聚类添加偏移
        X(cluster_points, :) = X(cluster_points, :) + i * 2;
        labels(cluster_points) = i;
    end
    
    % 比较标准t-SNE和快速t-SNE
    fprintf('\n比较标准t-SNE和快速t-SNE:\n');
    
    % 标准t-SNE
    tic;
    tsne_standard = FastTSNE(2, 30, 200, 500);
    Y_standard = tsne_standard.fit_transform(X, 'verbose', false);
    time_standard = toc;
    
    % 快速t-SNE
    tic;
    tsne_fast = FastTSNE_BH(2, 30, 200, 500);
    Y_fast = tsne_fast.fit_transform(X, 'verbose', false);
    time_fast = toc;
    
    fprintf('标准t-SNE时间: %.2f秒\n', time_standard);
    fprintf('快速t-SNE时间: %.2f秒\n', time_fast);
    fprintf('加速比: %.2fx\n', time_standard / time_fast);
    
    % 可视化比较
    figure('Position', [100, 100, 1200, 500]);
    
    subplot(1, 2, 1);
    scatter(Y_standard(:, 1), Y_standard(:, 2), 20, labels, 'filled');
    title(sprintf('标准t-SNE (%.2f秒)', time_standard));
    axis equal; colorbar;
    
    subplot(1, 2, 2);
    scatter(Y_fast(:, 1), Y_fast(:, 2), 20, labels, 'filled');
    title(sprintf('快速t-SNE (%.2f秒)', time_fast));
    axis equal; colorbar;
    
    sgtitle('大规模数据t-SNE比较');
end

function demo_comparison_pca()
    % 与PCA比较
    
    fprintf('\nt-SNE vs PCA 比较演示\n');
    fprintf('====================\n');
    
    % 生成非线性数据（瑞士卷）
    [X, labels] = generate_swiss_roll();
    
    % t-SNE降维
    tic;
    tsne = FastTSNE(2, 30, 200, 1000);
    Y_tsne = tsne.fit_transform(X, 'verbose', false);
    time_tsne = toc;
    
    % PCA降维
    tic;
    Y_pca = pca_implementation(X, 2);
    time_pca = toc;
    
    fprintf('t-SNE时间: %.2f秒\n', time_tsne);
    fprintf('PCA时间: %.2f秒\n', time_pca);
    
    % 可视化比较
    figure('Position', [100, 100, 1500, 500]);
    
    % 原始数据（3D）
    subplot(1, 3, 1);
    scatter3(X(:, 1), X(:, 2), X(:, 3), 30, labels, 'filled');
    title('原始数据 (3D瑞士卷)');
    axis equal; view(45, 30); colorbar;
    
    % PCA结果
    subplot(1, 3, 2);
    scatter(Y_pca(:, 1), Y_pca(:, 2), 30, labels, 'filled');
    title(sprintf('PCA降维 (%.2f秒)', time_pca));
    axis equal; colorbar;
    
    % t-SNE结果
    subplot(1, 3, 3);
    scatter(Y_tsne(:, 1), Y_tsne(:, 2), 30, labels, 'filled');
    title(sprintf('t-SNE降维 (%.2f秒)', time_tsne));
    axis equal; colorbar;
    
    sgtitle('非线性数据降维方法比较');
end

function [X, t] = generate_swiss_roll()
    % 生成瑞士卷数据
    
    n_samples = 1000;
    t = 3 * pi / 2 * (1 + 2 * rand(n_samples, 1));
    height = 21 * rand(n_samples, 1);
    
    X = zeros(n_samples, 3);
    X(:, 1) = t .* cos(t);
    X(:, 2) = height;
    X(:, 3) = t .* sin(t);
    
    % 添加噪声
    X = X + 0.05 * randn(size(X));
end

function Y = pca_implementation(X, n_components)
    % PCA实现
    X_centered = X - mean(X, 1);
    [coeff, ~, ~] = pca(X_centered);
    Y = X_centered * coeff(:, 1:n_components);
end

% 运行主演示
main_tsne_demo();

算法加速技巧

1. Barnes-Hut近似

• 使用空间分割树（四叉树/八叉树）
• 将远距离点分组计算
• 复杂度从O(N²)降低到O(N log N)

2. 早期放大策略

• 前100次迭代使用放大的相似度
• 帮助聚类快速形成
• 提高收敛速度

3. 自适应学习率

• 使用增益矩阵调整学习率
• 动量项加速收敛
• 避免振荡

参数调优指南

关键参数影响

实用建议

1. 数据预处理：标准化数据，避免尺度差异
2. 多次运行：t-SNE有随机性，多次运行取最佳
3. 聚类数估计：困惑度≈期望的邻居数量
4. 可视化验证：结合领域知识解释结果

参考代码 tsne 快速降维算法 www.3dddown.com/csa/78930.html

应用场景

适合使用t-SNE：

• 高维数据可视化（基因表达、文本数据、图像特征）
• 聚类结构探索
• 异常检测
• 数据质量评估

不适合使用t-SNE：

• 特征选择（使用PCA/LDA）
• 降维后建模（使用UMAP/Autoencoder）
• 精确距离保持（使用MDS）

高级特性

这个实现包含以下高级功能：

1. 完整的收敛监控
2. 参数敏感性分析
3. 与PCA的对比
4. 大规模数据优化
5. 丰富的可视化