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% 更新方向endmethodsfunctionobj=FastTSNE(num_dimensions,perplexity,learning_rate,max_iter)% 构造函数ifnargin<4max_iter=1000;endifnargin<3learning_rate=200;endifnargin<2perplexity=30;endifnargin<1num_dimensions=2;endobj.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;endfunction[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);foriter=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);% 早期放大阶段结束ifiter==100P=obj.P;end% 进度显示ifp.Results.verbose&&mod(iter,100)==0fprintf('迭代 %d/%d, 代价: %.4f\n',iter,obj.max_iter,costs(iter));end% 绘制进度ifp.Results.plot_progress&&mod(iter,50)==0obj.plot_current_embedding(iter,costs);endendY=obj.Y;fprintf('t-SNE完成！最终代价: %.4f\n',costs(end));endfunctionX_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;% 处理常数特征endfunctionP=compute_high_dimensional_similarities(obj,X)% 计算高维空间中的相似度概率n=size(X,1);P=zeros(n);% 计算成对距离D=pdist2(X,X).^2;% 为每个点寻找合适的方差（困惑度）fori=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);endfunctionbeta=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);ifsum_P==0sum_P=1e-12;endP=P/sum_P;H=obj.compute_entropy(P);perplexity_diff=H-log(target_perplexity);iter=0;whileabs(perplexity_diff)>tol&&iter<max_iterifperplexity_diff>0beta_min=beta;ifisinf(beta_max)beta=beta*2;elsebeta=(beta+beta_max)/2;endelsebeta_max=beta;ifisinf(beta_min)beta=beta/2;elsebeta=(beta+beta_min)/2;endendP=exp(-D_row*beta);P(D_row==0)=0;sum_P=sum(P);ifsum_P==0sum_P=1e-12;endP=P/sum_P;H=obj.compute_entropy(P);perplexity_diff=H-log(target_perplexity);iter=iter+1;endendfunctionH=compute_entropy(~,P)% 计算概率分布的香农熵P(P==0)=1e-12;% 避免log(0)H=-sum(P.*log(P));endfunctionY=initialize_embedding(obj,n_samples)% 初始化低维嵌入Y=1e-4*randn(n_samples,obj.num_dimensions);endfunction[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));fori=1:n% 向量化梯度计算diff=Y(i,:)-Y;grad(i,:)=4*sum(bsxfun(@times,PQ(:,i),diff),1);endendfunctionobj=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);endfunctioncost=compute_cost(~,P,Q)% 计算KL散度代价cost=sum(P(:).*log(P(:)./Q(:)));endfunctionplot_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;endendmethods(Static)functiondemo_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);endfunction[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=[];fori=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)];endX=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);endfunctionplot_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降维效果展示');endfunctioncompare_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);fori=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);fori=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];fori=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参数比较');endendend% 快速t-SNE实现（Barnes-Hut近似）classdef FastTSNE_BH<FastTSNE% 基于Barnes-Hut近似的快速t-SNE实现% 适用于大规模数据集methodsfunctionobj=FastTSNE_BH(num_dimensions,perplexity,learning_rate,max_iter)% 构造函数obj=obj@FastTSNE(num_dimensions,perplexity,learning_rate,max_iter);obj.momentum=0.5;endfunction[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);% 优化过程foriter=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);% 调整学习率ifiter>100obj.learning_rate=obj.learning_rate*0.99;endifp.Results.verbose&&mod(iter,100)==0fprintf('迭代 %d/%d, 代价: %.4f\n',iter,obj.max_iter,costs(iter));endendY=obj.Y;endfunctionP=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);fori=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,:));endP=(P+P')/(2*n);P=max(P,1e-12);endfunction[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);% 计算梯度fori=1:n[grad_i,cost_i]=obj.compute_gradient_for_point(P,Y,i,tree,theta);grad(i,:)=grad_i;cost=cost+cost_i;endcost=cost/n;endfunctiontree=build_tree(~,Y)% 构建空间分割树（简化实现）% 实际实现应该更复杂，这里使用简化版本tree.points=Y;tree.center=mean(Y,1);tree.size=max(range(Y,1));tree.num_points=size(Y,1);endfunction[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);% 计算梯度和代价forj=1:nifj~=idiff=Y(i,:)-Y(j,:);mult=(P(i,j)-Q(j))*Q(j);grad=grad+4*mult*diff;ifP(i,j)>0cost=cost+P(i,j)*log(P(i,j)/Q(j));endendendendendend% 主测试函数functionmain_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): ');switchchoicecase1FastTSNE.demo_tsne();case2demo_large_scale();case3demo_comparison_pca();otherwisefprintf('无效选择，运行基础演示...\n');FastTSNE.demo_tsne();endendfunctiondemo_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);fori=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-SNEfprintf('\n比较标准t-SNE和快速t-SNE:\n');% 标准t-SNEtic;tsne_standard=FastTSNE(2,30,200,500);Y_standard=tsne_standard.fit_transform(X,'verbose',false);time_standard=toc;% 快速t-SNEtic;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比较');endfunctiondemo_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('非线性数据降维方法比较');endfunction[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));endfunctionY=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.自适应学习率

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

参数调优指南

关键参数影响

参数	推荐范围	影响
困惑度	5-50	控制局部/全局结构平衡
学习率	10-1000	影响收敛速度和稳定性
迭代次数	500-1000	保证充分收敛
动量	0.5-0.8	加速收敛，防止振荡