大一线代的综合作业,一个简单运用pca(主成分分析法)的作业,给的数据是一个高光谱数据,不过老师已经把先它变成二维矩阵了,记录一下。
目录
2.求协方差矩阵G的特征值,特征向量。并且降序特征值,同时特征向量一一对应
1.载入数据,对数据去均值,求矩阵的协方差矩阵
1)调用load函数,在工作区得知数据储存在X里,并且X是一个65536*191的矩阵
load('线代作业\03\HyperspectralImage.mat');
2)去均值
%去均值
Q = X - repmat(mean(X),65536,1);
3)求协方差矩阵G,协方差矩阵 = 1/(n-1)*QT*Q,其中n为列数
G = 1/190*(Q')*Q;%求协方差矩阵
2.求协方差矩阵G的特征值,特征向量。并且降序特征值,同时特征向量一一对应
1)调用eig函数,函数的返回值是特征向量和特征值组成的对角矩阵,其中特征向量为V,特征值组成的对角矩阵为D
[V,D] = eig(G);%V表示G的特征向量,D表示特征值组成的对角矩阵
2)对特征值进行降序处理,首先调用diag数,提取矩阵对角线上的数(即提取特征值),赋予到变量y中
y = diag(D);%提取特征值
然后调用sort函数,对y进行降序,结果储存到变量Y中,索引为index。
[Y,index] = sort(y,'descend');%对特征值降序,其中Y为降序后的列向量,index为索引
3)利用索引,将特征向量与变化后的特征值一一对应,结果赋予到变量V_index
V_index = V(:,index);%将原特征向量根据新顺序的特征值进行排序
3.根据信息含量计算主成分数量k
需要保留99.8%的信息,寻找k考虑使用for循环和if语句
%运用循环语句寻找K值
for k=1:1:191
Y_k = Y(1:k,1);%选择前k个特征值
N = sum(Y_k);%前k个特征值的和
M = sum(Y);%所有特征值的和
message = N/M;%信息含量
if message >= 0.998
disp(k)
break
else
continue
end
end
运行,得到k的值是6
4.显示这6个图像
V_index_6 = V_index(:,1:6);%取前六个主成分
U = Q * V_index_6;%把原始数据投影到新的坐标系上
U_scale = rescale(U);%缩放数据
immg = U_scale(:,1:6);%提取前6个图像
mul = reshape(immg,256,[]);%升维
imshow(mul);%显示图像
运行,得到图像
通过主成分变换的结果可以明显观察出,贡献率前6的图像已经可以基本完全表示原191张图像所呈现的不同特征,达到了数据低损耗压缩的效果。可以观察出,输出结果的6张图像基本线性无关,尤其是第1张和第2张图像,清晰地呈现了几乎完全不同的特征。
附完整代码:
load('线代作业\03\HyperspectralImage.mat');
%去均值
Q = X - repmat(mean(X),65536,1);
G = 1/190*(Q')*Q;%求协方差矩阵
%计算协方差矩阵G的特征值和特征向量
[V,D] = eig(G);%V表示G的特征向量,D表示特征值组成的对角矩阵
%对特征值进行降序处理,并且对应特征向量,然后对特征向量施行正交化
y = diag(D);%提取特征值
[Y,index] = sort(y,'descend');%对特征值降序,其中Y为降序后的列向量,index为索引
V_index = V(:,index);%将原特征向量根据新顺序的特征值进行排序
%运用循环语句寻找K值
for k=1:1:191
Y_k = Y(1:k,1);%选择前k个特征值
N = sum(Y_k);%前k个特征值的和
M = sum(Y);%所有特征值的和
message = N/M;%信息含量
if message >= 0.998
disp(k)
break
else
continue
end
end
V_index_6 = V_index(:,1:6);%取前六个主成分
U = Q * V_index_6;%把原始数据投影到新的坐标系上
U_scale = rescale(U);%缩放数据
immg = U_scale(:,1:6);%提取前6个图像
mul = reshape(immg,256,[]);%升维
imshow(mul);%显示图像