降维算法

PCA算法

• 算法流程:最小的m-l个特征值被舍弃了，这是降维的结果。舍弃部分信息后使得样本的1信息密度变大；2后面的m-l个特征是噪声，有去噪的效果。
  

对数据集进行PCA后再重构

• 1通过上面算法降维成l个维度l <<m，m为之前的维度。
  
• 2再进行升维重构
  

PCA降维-Matlab代码

例子1-系统pca做降维pca_mat.m

• 改例子有计算过程。

%使用pca函数进行降维
X=[1,2;3,3;3,5;5,4;5,6;6,5;8,7;9,8]; %样本矩阵,每行一个样本向量
[coeff,  score,  latent]=pca(X); %主成分分析
coeff, %主成分分量 (每列为一个变换空间中的基向量)
score, %主成分, score(:,1) 为X的一维表示, score为X在变换空间中的二维表示
latent, %X样本协方差矩阵的特征值

例子2-系统pca对简单数据集降维pca2_mat.m

• 13*4的数据进行降维

clc, clear all, close all
load hald %载入数据集
X=ingredients; %将数据集中的13x4矩阵矩阵赋给X
[coeff,score]=pca(X,'Centered','off'), %调用函数pca, 不对数据中心化

例子3-用fastPCA做降维

• fastPCA.m

%快速PCA函数
function [V,pcaX]=fastPCA(X,k)
%输入: X--样本矩阵,每行为一个样本
%输出: V--主成分向量
%pcaX--降维后的k维样本特征向量组成的矩阵,每一行一个样本,列数k为降维后的样本特征列数
[n,m]=size(X);
meanV=mean(X); %样本均值向量
Z=(X-repmat(meanV,n,1));
CovM=Z*Z';  %计算协方差矩阵的转置
[V,D]=eigs(CovM,k); %计算CovM的前k个特征值和特征向量
V=Z'*V; %得到协方差矩阵(CovM)'的特征向量
%特征向量归一化为单位特征向量
for i=1:k
    V(:,i)=V(:,i)/norm(V(:,i));
end
pcaX=Z*V; %线性变换(投影)降维至k维
save('PCA.mat','V','meanV'); %保存变换矩阵V和变换原点meanV

• pca_iris.m对鸢尾属植物数据集降至2维

%用fastPCA函数对鸢尾属植物数据集降至2维
clc, clear all, close all
load fisheriris
X = meas;  Y = species;
[n,m]=size(X);
[V,pcaX]=fastPCA(X,2);
gscatter(pcaX(:,1),pcaX(:,2),Y,'rbm','*vo');

流行学习

• 由于实际应 用中很多时候数据是非线性的, 因此有必要发展非线性降维算法, 流形学习就是一种非 线性降维技术的典型代表. “流形”是几何学中的一个概念, 它是高维空间中的几何结构, 即空间中的点构成 的集合, 例如二维空间中的曲线、三维空间中的曲面都是流形的例子。

局部线性嵌入

• lle.m

%LLE MATLAB程序
function [Z] = lle(X,k,l)
%X-样本矩阵,m*n 矩阵,m为样本向量的维数,n为样本的个数
%k-近邻数,l-低维空间维数,Z为降维后的数据矩阵,维度为l*n
[m,n] = size(X);
fprintf(1,'LLE运行于 %d 维空间中的 %d 个样本点上\n',m,n);
%第一步: 计算点与点之间的距离并寻找k个近邻点.
fprintf(1,'-->寻找 %d 个近邻点.\n',k);
X2 = sum(X.^2,1);
distance = repmat(X2,n,1)+repmat(X2',1,n)-2*X'*X;
%欧氏距离公式的展开||x_i-x_j||^2=x_i^2+x_j^2-2x_i'x_j
[~,index] = sort(distance); %距离排序(从小到大)
neighborhood = index(2:(1+k),:); %取第2到第k+1个最小距离
%第二步: 计算重构权重矩阵
fprintf(1,'-->求解重构权重矩阵.\n');
if(k>m)  %如果近邻点个数大于样本特征向量维数, 需要正则化
    fprintf(1,'-->近邻点个数大于样本特征向量维数需要正则化.\n');
    tol=1e-3; 
else
    tol=0;
end
W = zeros(k,n); %权重矩阵初始化
for i=1:n
   z = X(:,neighborhood(:,i))-repmat(X(:,i),1,k); %将第i个点平移至坐标原点
   C = z'*z; %局部协方差矩阵
   C = C + eye(k,k)*tol*trace(C);  %k>m时加上正则化项
   W(:,i) = C\ones(k,1);   %解方程组Cw=1
   W(:,i) = W(:,i)/sum(W(:,i)); %按行归一化权重矩阵
end
%第三步: 计算矩阵M=(I-W)'(I-W)
fprintf(1,'-->计算局部线性嵌入.\n');
%M=eye(n,n); %M是一个只有4kn个非零元的稀疏矩阵
M = sparse(1:n,1:n,ones(1,n),n,n,4*k*n);
for i=1:n %计算矩阵M
   w = W(:,i);
   j = neighborhood(:,i);
   M(i,j) = M(i,j) - w';
   M(j,i) = M(j,i) - w;
   M(j,j) = M(j,j) + w*w';
end
%计算局部嵌入
options.disp = 0; options.isreal = 1; %结构变量options赋初值
[Z,~] = eigs(M,l+1,0,options); %计算矩阵M的前l+1个最小特征值对应的特征向量
Z = Z(:,2:l+1)'*sqrt(n); 
fprintf(1,'程序运行结束.\n');

局部线性嵌入-Matlab代码和使用lle_examp.m

%用瑞士卷数据测试LLE算法
clc, close all, clear all;
k = 20 ;%(邻域点个数)
l = 2;%最大嵌入维数
%瑞士卷的生成图
rng (3)
n=1500; 
t=(3*pi/2)*(1+2*rand(1,n));
s=21*rand(1,n);
X=[t.*cos(t); s; t.*sin(t)]; %3*n矩阵
figure (1)
plot3(X(1,:),X(2,:),X(3,:),'.'), %绘制瑞士卷三维曲线图
view([12 12]);  %设置视角
%调用lle函数降维
Z=lle(X,k,l);
figure (2)
plot(Z(1,:),Z(2,:),'.b'), %绘制LLE降维后的二维平面图

MDS降维

• 多维缩放方法(MDS)
  
  
  
• 最后推导
  

MDS降维-Matlab代码

• MDS.m

function Z = MDS(D, d)
%多维缩放方法(MDS)的MATLAB 程序
%功能: 支持n个m维样本降维值d维. 
%输入: D--n*n距离矩阵, d--低维空间维数
%输出: Z--低维空间的d*n 样本矩阵
[n, n1] = size(D);
if n~=n1
    printf('D不是一个方阵');
end
D2 = D .* D;  e = ones(n,1);
H = eye(n) - 1 / n * (e *e');
B =-0.5 * H * D2 * H; 
%B=zeros(n,n);
% for i=1:n
%     for j=1:n
%         B(i,j)=-0.5*(D(i,j)^2 -1/n*D(i,:)*D(i,:)' -1/n*D(:,j)'*D(:,j) +1/n^2*sum(sum(D.^2)));
%     end
% end
[U, Lam] = eig(B);
Lam = diag(Lam);
[~, pos] = sort(Lam,'descend');
index = pos(1:d);
U = U(:, index);
Z = diag(Lam(index))^(0.5) * U';

MDS降维

• Isomap 源于微分几何中测地线的思想, 它希望样本数据在向低 维空间投影之后仍能保持流形上的测地线距离.
• 最后需要指出的是, 近邻图的构建通常有两种方法: 一种是指定近邻点个数 𝑘, 将与 𝑥𝑖 欧氏距离最近的 𝑘 个点作为近邻点. 另一种方法是指定距离阈值 𝜖, 将与 𝑥𝑖 的欧氏距 离小于 𝜖 的点作为近邻点. 显然这两种方法均有不足: 若 𝑘 值过大, 则距离很远的点可 能被误认为是近邻点, 即出现所谓的“短路”现象; 若阈值 𝜖 过小, 则图中有些区域可能 与其他区域不存在连接, 即出现所谓的“断路”现象. 短路与断路都会给后续的最短路 径计算造成误导.
  

MDS降维-Matlab实现和使用

• isomap.m

function  Z = isomap(X, k, d)
%Isomap算法MATLAB程序, 该程序调用了MDS算法
%第一步: 按欧氏距离选取k个近邻点
[n, ~] = size(X); D = zeros(n);
for i =1 : n %构造距离矩阵D
    xx = repmat(X(i, :), n, 1); 
    %复制X的第i行n次叠成一个n*n矩阵, 其每一行元素都是X的第i行元素
    diff = xx - X;
    dist = sum(diff.* diff, 2); %按行求和
    [dd, pos] = sort(dist); %按升序排序
    index = pos(1 : k + 1)'; %前k+1个设为近邻点
    index2 = pos(k + 2 : n); %后n-k-1设为非近邻点
    D(i,index) = sqrt(dd(index)); %距离矩阵第i行近邻点元素值
    D(i, index2) = inf; %距离矩阵第i行非近邻点元素值设为无穷大
end
%第二步:  重新计算最短距离矩阵
for k=1:n
    for i=1:n
        for j=1:n
            if D(i,j)>D(i,k)+D(k,j)
                D(i,j)=D(i,k)+D(k,j);
            end
        end
    end
end
%第三步: 调用MDS算法降维
Z = MDS(D, d);

• isomap_examp.m

%用瑞士卷数据测试Isomap算法
close all, clear all;
k = 20 ;%(邻域点个数)
l = 2;%最大嵌入维数
%瑞士卷的生成图
rng (3)
n=1500; 
t=(3*pi/2)*(1+2*rand(1,n));
s=21*rand(1,n);
X=[t.*cos(t); s; t.*sin(t)]; %3*n矩阵
figure (1)
plot3(X(1,:),X(2,:),X(3,:),'.'), %绘制瑞士卷三维曲线图
view([12 12]);  %设置视角
X = X';
Z=isomap(X,k,l);  
figure (2)
plot(Z(1,:),Z(2,:),'.b'), %绘制Isomap降维后的二维平面图