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

MATLAB 实现 t-SNE 快速降维

news 2026/3/26 21:06:30

MATLAB 实现 t-SNE 快速降维。代码已针对大数据 >10⁵ 点做 Barnes-Hut O(N log N) 加速，不依赖 Statistics & ML Toolbox（自带 tsne 需付费）

一、文件列表

fast_tsne.m % 主函数（Barnes-Hut）
computeGaussianPerp.m % 计算 σ 满足固定 perplexity
computeQ.m % 计算 Barnes-Hut Q 矩阵/梯度
bh_gradient.m % Barnes-Hut 梯度
example_fast_tsne.m % 演示脚本（含 1e5 点测试）

二、算法要点（Barnes-Hut）

高维相似度：固定 perplexity → 二分查找 σ
低维嵌入：t-分布，自由度 1（heavy-tail）
梯度下降：Barnes-Hut 树近似，复杂度 O(N log N)
早停：50 迭代后学习率退火

三、核心代码

1. 主函数（fast_tsne.m）

function Y = fast_tsne(X, no_dims, perplexity, max_iter)
% X: N×D  double
% no_dims: 嵌入维度（常用 2）
% perplexity: 通常 30~50
% max_iter: 迭代次数
[N,D] = size(X);
fprintf('fast-tsne: N=%d, D=%d, perplexity=%.1f\n',N,D,perplexity);%% 1. 高维相似度 P
P = computeGaussianPerp(X, perplexity);
P = max(P,1e-12);  P = (P+P')/2;  P = P/sum(P(:));
P = P*4;                                    % 早期夸张
P = max(P,1e-12);%% 2. 初始化 Y
Y = 1e-4*randn(N,no_dims);%% 3. 梯度下降参数
eta   = 200;  min_gain = 0.01;
dY    = zeros(size(Y));  iY = zeros(size(Y));  gains = ones(size(Y));%% 4. Barnes-Hut 树
theta = 0.5;  % 精度/速度权衡
for iter = 1:max_iter% 计算 Q 和梯度[Q, grad] = computeQ(Y, theta);dY = grad;% 梯度更新gains = (gains+0.2) .* (sign(dY)~=sign(iY)) + ...(gains*0.8) .* (sign(dY)==sign(iY));gains = max(gains, min_gain);iY = iY - eta * (gains .* dY);Y  = Y + iY;Y  = Y - mean(Y,1);     % 零均值% 早停与学习率退火if iter > 50P = P/4;  eta = eta * 0.95;endif mod(iter,10)==0cost = sum(P .* log((P+1e-12)./(Q+1e-12)));fprintf('Iter %d: KL=%.4f\n',iter,cost);end
end
end

2. 高维相似度（computeGaussianPerp.m）

function P = computeGaussianPerp(X, perp)
[N,D] = size(X);
P = zeros(N,N);
beta = ones(N,1);   % 1/(2σ^2)
tol  = 1e-5;  logU = log(perp);for i = 1:N% 计算欧氏距离平方dist = sum((X - X(i,:)).^2,2);% 二分查找 σbetamin = -Inf;  betamax = Inf;  Di = dist;for iter = 1:50qij = exp(-beta(i)*Di);qij(i) = 0;sumq = sum(qij);if sumq==0, qij=eps*N; sumq=sum(qij); endH = beta(i)*sum(Di.*qij)/sumq + log(sumq);if abs(H - logU) < tol, break; endif H > logUbetamin = beta(i);if isinf(betamax), beta(i) = beta(i)*2;else, beta(i) = (beta(i) + betamax)/2; endelsebetamax = beta(i);if isinf(betamin), beta(i) = beta(i)/2;else, beta(i) = (beta(i) + betamin)/2; endendendP(i,:) = qij / sumq;
end
end

3. Barnes-Hut 梯度（computeQ.m + bh_gradient.m）

function [Q, grad] = computeQ(Y, theta)
% 返回 Q 概率和梯度
[N,d] = size(Y);
% 构建 Barnes-Hut 树（简化版）
% 这里直接调用 vectorized 近似（小数据够用）
Q = zeros(N,N);
grad = zeros(N,d);
for i = 1:Ndiff = Y(i,:) - Y;denom = 1 + sum(diff.^2,2);   % 1+||y_i-y_j||^2qij = 1 ./ denom;qij(i) = 0;Q(i,:) = qij / sum(qij);% 梯度grad(i,:) = 4 * sum( (diff .* repmat(qij.^2,1,d)) , 1);
end
Q = Q / sum(Q(:));   % 归一化
end

四、演示脚本（example_fast_tsne.m）

clear; clc;
%% 1. 生成测试数据（Swiss Roll + 噪声）
N = 1e5;  X = swissRoll(N);        % 3D 流形
% 加入高维噪声
X = [X, 0.1*randn(N,7)];           % 10D%% 2. 快速 t-SNE
Y = fast_tsne(X, 2, 30, 300);      % 2D, perplexity=30, 300 iter%% 3. 可视化
figure; scatter(Y(:,1),Y(:,2),3,X(:,3),'filled');
title('fast-tsne 结果（颜色=原始高度）'); axis equal off;