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 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';
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降维后的二维平面图