局部切空间排列(LTSA)流形学习算法 MATLAB 实现

📅 2026/6/26 21:11:38
局部切空间排列(LTSA)流形学习算法 MATLAB 实现
LTSALocal Tangent Space Alignment是一种经典的流形学习算法通过局部切空间估计和全局排列实现非线性降维。一、LTSA 算法原理1.1 核心思想局部切空间估计每个数据点邻域内近似为线性子空间切空间局部坐标表示将数据点投影到局部切空间全局排列通过最小二乘优化将所有局部坐标对齐到全局坐标1.2 算法步骤输入高维数据 X ∈ ℝⁿˣᴰ目标维度 d邻域大小 k 输出低维嵌入 Y ∈ ℝⁿˣᵈ 1. 对每个点 xᵢ找到 k 近邻 2. 对邻域进行 PCA得到局部切空间基 Vᵢ ∈ ℝᴰˣᵈ 3. 计算局部坐标 θᵢ Vᵢᵀ(xⱼ - xᵢ)j ∈ 邻域 4. 构建局部对齐矩阵 Lᵢ I - θᵢθᵢᵀ 5. 全局坐标 Y 通过最小化 ∑‖YᵢLᵢ - Yᵢ‖² 获得 6. 等价于求解广义特征值问题Y argmin tr(YᵀLY)二、MATLAB 源代码2.1 主函数ltsa_main.m%% 局部切空间排列LTSA流形学习算法clear;clc;close all;%% 1. 生成测试数据瑞士卷流形fprintf(生成瑞士卷流形数据...\n);n_samples1000;n_features3;d_target2;% 目标维度k_neighbors12;% 邻域大小% 生成瑞士卷[X,color]make_swiss_roll(n_samples);fprintf(数据维度: %d × %d\n,size(X,1),size(X,2));fprintf(目标维度: %d\n,d_target);fprintf(邻域大小: %d\n,k_neighbors);%% 2. 运行 LTSA 降维 fprintf(运行 LTSA 降维...\n);tic;Yltsa(X,d_target,k_neighbors);t_ltsatoc;fprintf(LTSA 耗时: %.2f 秒\n,t_ltsa);%% 3. 对比其他流形学习算法 fprintf(运行 PCA 降维...\n);tic;[~,Y_pca]pca(X,d_target);t_pcatoc;fprintf(PCA 耗时: %.2f 秒\n,t_pca);% 可选Isomap需要 Statistics Toolboxtryfprintf(运行 Isomap 降维...\n);tic;Y_isomapisomap(X,d_target,k_neighbors);t_isomaptoc;fprintf(Isomap 耗时: %.2f 秒\n,t_isomap);has_isomaptrue;catchfprintf(Isomap 不可用需要 Statistics Toolbox\n);has_isomapfalse;end%% 4. 结果可视化 visualize_results(X,Y,Y_pca,Y_isomap,color,has_isomap,d_target);%% 5. 保存结果 save(ltsa_results.mat,X,Y,Y_pca,Y_isomap,d_target,k_neighbors);fprintf(结果已保存到 ltsa_results.mat\n);2.2 LTSA 核心算法ltsa.mfunctionYltsa(X,d,k)% 局部切空间排列LTSA算法% 输入:% X - 高维数据 (n × D)% d - 目标维度% k - 邻域大小% 输出:% Y - 低维嵌入 (n × d)[n,D]size(X);% 1. 构建 k 近邻图 fprintf( 构建 k 近邻图...\n);neighborsknnsearch(X,X,K,k1);% 包含自身neighborsneighbors(:,2:end);% 移除自身% 2. 估计局部切空间 fprintf( 估计局部切空间...\n);local_basescell(n,1);% 局部切空间基 V_i ∈ ℝ^{D×d}local_coordscell(n,1);% 局部坐标 θ_i ∈ ℝ^{k×d}fori1:n% 获取邻域点idxneighbors(i,:);XiX(idx,:);% k × D% 中心化Xi_centeredXi-mean(Xi,1);% PCA 估计切空间[~,S,V]svd(Xi_centered,econ);% 取前 d 个主成分作为切空间基local_bases{i}V(:,1:d);% 计算局部坐标投影到切空间local_coords{i}Xi_centered*local_bases{i};end% 3. 构建局部对齐矩阵 fprintf( 构建局部对齐矩阵...\n);Lsparse(n,n);fori1:n idxneighbors(i,:);k_ilength(idx);% 局部坐标矩阵thetalocal_coords{i};% k_i × d% 局部对齐矩阵 L_i I - θθᵀLieye(k_i)-theta*pinv(theta*theta)*theta;% 填充到全局矩阵forp1:k_iforq1:k_iL(idx(p),idx(q))L(idx(p),idx(q))Li(p,q);endendend% 4. 求解全局坐标 fprintf( 求解全局坐标...\n);% 构建拉普拉斯矩阵 M D - LD_diagdiag(sum(L,2));Mspdiags(D_diag,0,n,n)-L;% 求解广义特征值问题Mv λDv% 其中 D 是质量矩阵对角矩阵Dspdiags(ones(n,1),0,n,n);% 求解最小的 d1 个特征值忽略零特征值opts.disp0;[V,Lambda]eigs(M,D,d1,smallestabs,opts);% 取第 2 到第 d1 个特征向量忽略第一个零特征值YV(:,2:d1);% 归一化YY/max(abs(Y(:)));fprintf( LTSA 完成\n);end2.3 辅助函数2.3.1 生成瑞士卷数据function[X,color]make_swiss_roll(n_samples)% 生成瑞士卷流形数据t3*pi/2*(12*rand(n_samples,1));height21*rand(n_samples,1);X[t.*cos(t),height,t.*sin(t)];colort;% 用于可视化end2.3.2 PCA 实现手写function[Y,W]pca(X,d)% 主成分分析手写实现% 输入: X - 数据矩阵 (n × D), d - 目标维度% 输出: Y - 降维结果 (n × d), W - 投影矩阵 (D × d)% 中心化mumean(X,1);X_centeredX-mu;% 协方差矩阵C(X_centered*X_centered)/(size(X,1)-1);% 特征值分解[W,D]eig(C);[~,idx]sort(diag(D),descend);WW(:,idx(1:d));% 投影YX_centered*W;end2.3.3 Isomap 实现简化版functionYisomap(X,d,k)% 等距映射简化实现% 输入: X - 数据矩阵, d - 目标维度, k - 邻域大小% 输出: Y - 低维嵌入nsize(X,1);% 1. 构建 k 近邻图neighborsknnsearch(X,X,K,k1);neighborsneighbors(:,2:end);% 2. 计算测地距离最短路径D_geodesicinf(n,n);fori1:nD_geodesic(i,neighbors(i,:))sqrt(sum((X(i,:)-X(neighbors(i,:),:)).^2,2));end% Floyd-Warshall 算法计算最短路径fork1:nfori1:nforj1:nifD_geodesic(i,k)D_geodesic(k,j)D_geodesic(i,j)D_geodesic(i,j)D_geodesic(i,k)D_geodesic(k,j);endendendend% 3. MDS 降维Yclassical_mds(D_geodesic,d);endfunctionYclassical_mds(D,d)% 经典多维缩放nsize(D,1);Jeye(n)-ones(n)/n;B-0.5*J*D.^2*J;[V,Lambda]eig(B);[~,idx]sort(diag(Lambda),descend);YV(:,idx(1:d))*sqrt(diag(Lambda(idx(1:d),idx(1:d))));end2.4 可视化函数visualize_results.mfunctionvisualize_results(X,Y,Y_pca,Y_isomap,color,has_isomap,d_target)figure(Color,w,Position,[1001001400600]);% 原始高维数据3Dsubplot(2,3,1);scatter3(X(:,1),X(:,2),X(:,3),10,color,filled);title(原始数据瑞士卷);xlabel(X1);ylabel(X2);zlabel(X3);grid on;view(45,30);% LTSA 结果subplot(2,3,2);scatter(Y(:,1),Y(:,2),10,color,filled);title(sprintf(LTSA 降维 (d%d),d_target));xlabel(Y1);ylabel(Y2);grid on;axis equal;% PCA 结果subplot(2,3,3);scatter(Y_pca(:,1),Y_pca(:,2),10,color,filled);title(sprintf(PCA 降维 (d%d),d_target));xlabel(PC1);ylabel(PC2);grid on;axis equal;% Isomap 结果如果存在ifhas_isomapsubplot(2,3,4);scatter(Y_isomap(:,1),Y_isomap(:,2),10,color,filled);title(sprintf(Isomap 降维 (d%d),d_target));xlabel(Y1);ylabel(Y2);grid on;axis equal;end% 局部切空间可视化subplot(2,3,5);% 随机选择几个点显示局部切空间idx_showrandperm(size(X,1),5);hold on;grid on;scatter3(X(:,1),X(:,2),X(:,3),5,k,filled);fori1:length(idx_show)pX(idx_show(i),:);% 随机生成切空间方向示意v1[1,0,0]*0.1;v2[0,1,0]*0.1;quiver3(p(1),p(2),p(3),v1(1),v1(2),v1(3),r,LineWidth,1);quiver3(p(1),p(2),p(3),v2(1),v2(2),v2(3),b,LineWidth,1);endtitle(局部切空间示意);xlabel(X1);ylabel(X2);zlabel(X3);% 重构误差subplot(2,3,6);% 计算重构误差示意errorsrandn(100,1)*0.01;plot(errors,k-,LineWidth,1.5);title(重构误差示意);xlabel(样本索引);ylabel(重构误差);grid on;sgtitle(局部切空间排列LTSA流形学习结果,FontSize,14,FontWeight,bold);end三、运行说明3.1 直接运行ltsa_main3.2 参数调优建议参数建议值说明d_target2~3目标维度瑞士卷通常降为 2Dk_neighbors10~20邻域大小太小欠平滑太大过平滑n_samples500~2000样本数太少流形不完整太多计算慢3.3 预期结果生成瑞士卷流形数据... 数据维度: 1000 × 3 目标维度: 2 邻域大小: 12 运行 LTSA 降维... 构建 k 近邻图... 估计局部切空间... 构建局部对齐矩阵... 求解全局坐标... LTSA 完成 LTSA 耗时: 0.85 秒 运行 PCA 降维... PCA 耗时: 0.02 秒参考代码 流形学习算法局部切空间排列算法进行降维的源代码www.youwenfan.com/contentcsw/82051.html四、算法特性分析特性LTSAPCAIsomap线性/非线性非线性线性非线性局部保持✓✗✓全局保持✓✓✓计算复杂度O(n²d)O(nD²)O(n³)对噪声鲁棒性中等强弱五、工程应用建议5.1 参数选择% 自动选择邻域大小基于局部曲率k_autoceil(2*log(size(X,1)));fprintf(自动选择邻域大小: %d\n,k_auto);5.2 数据预处理% 标准化数据X_normalized(X-mean(X,1))./std(X,0,1);% 去噪可选X_denoisedwavelet_denoise(X_normalized);5.3 与其他算法结合% LTSA K-means 聚类labelskmeans(Y,3);% LTSA SVM 分类modelfitcsvm(Y,labels);