UFLDL教程答案:Exercise:PCA_in_2D&PCA_and_Whitening

Posted 2022-12-15 slim1017

教程地址：http://deeplearning.stanford.edu/wiki/index.php/UFLDL%E6%95%99%E7%A8%8B

练习地址1：http://deeplearning.stanford.edu/wiki/index.php/Exercise:PCA_in_2D

练习地址2： http://deeplearning.stanford.edu/wiki/index.php/Exercise:PCA_and_Whitening

1.习惯性啰嗦几句

（1）通过svd来避免求高维度数据的协方差矩阵

       PCA的推导有挺多版本的，ng的CS229中PCA那章讲过，还介绍了通过svd来避免求高维度数据的协方差矩阵(高维度协方差矩阵计算量很大)，所以这次PCA_in_2D练习中我写了两个版本：

       一个是pca_2d.m是按教程写的标准答案（但是我觉得教程绕了几个弯子,也可能是另有深意，我没理解到）。

       一个是another_pca_2d.m，这是我按照我自己的理解，使用svd绕开求协方差矩阵这一步来计算pca：

       建议先看看矩阵分析课本上关于奇异值分解(svd)的推导过程。

       首先说说[u,s,v]=svd(x),s对角线上为x的奇异值，也是矩阵xx'特征值的平方根；u为矩阵xx'的特征向量。我们用Sigma代表数据x的协方差矩阵。

       若X为n*m矩阵，n维特征，m个样本。m表示样本数，(xx'）/m其实就是Sigma(实际上应该是Sigma的样本估计，样本多就可认为就是Sigma)，也就是说可以直接通过s求得Sigma的特征值（把s对角线元素平方再除以m即可），这样计算，过程中并没有算Sigma具体是多少。

       此外，PCA得出的最佳投影方向为Sigma的特征向量方向，xx'与(xx'）/m特征向量一样，所以u就是旋转矩阵。

       而在教程中，有代码：

sigma = x * x' / size(x, 2)
[U,S,V] = svd(sigma);

       U为sigma*sigma’的特征向量矩阵，S对角线为sigma的奇异值(也等于sigma*sigma‘的特征值开平方根)。

       注意：先计算了协方差矩阵sigma，再用svd计算sigma的奇异值(S的对角线元素)，sigma奇异值就等于sigma的特征值。理由：根据矩阵分析中svd的推导过程，会先计算sigma*sigma‘的特征值，再开平方根，就得到了S(的对角线元素)；而由于sigma是实对称矩阵，sigma=sigma’，其实求得是sigma^2的特征值，再开根，其实就是sigma的特征值。此外，根据svd推导过程，U为sigma*sigma‘（也就是sigma^2）的特征向量组成的矩阵，那U同时也是sigma的特征向量矩阵，所以u和U中各向量其实是一样的（符号可能相反）。（矩阵A与A^2的特征向量都相同，特征值后者是前者的平方）。

       综上：教程中使用svd是为了求sigma的特征值和特征向量(用eig函数也可以，只是svd数值计算上更稳定)，而实际上，根据CS229-PCA那章，svd可以直接绕过求sigma这步，而直接求得sigma的特征值和特征向量，具体实现看后面的another_pca_2d.m。

 （2）PCA只有在x各行均值为0的前提下才成立

     关于数据均值规整化：对x行求和再平均是使x各列均值为0；对x的列求和再平均是使x各行均值为0。

(a) 教程中有这样两句代码，作用是把数据均值归0。若x为n*m矩阵，每列为n维特征，每行表示m个样本。教程用的是教程中有一小节“对图像数据应用PCA算法”中提到的自然图像常用的均值规整化方法，可是，在pca_2d的练习中，数据并不是图像，我认为不应该使用这种方法。在pca_2d中，我认为应该对x各行的均值归0，最后的运行结果和教程标准答案完全一样，证明了我的想法。

     其实PCA的推导过程中是要求x各行均值为0的，理由见(b)，而对于某自然图像数据x，其各行的均值天生就接近0的，所以就算不把x各行均值归0也可以使用PCA。其实pca_2d中的数据虽然不是图像数据，其各行均值也恰好接近0，所以就算不把各行的均值归0也可以（勉强）用PCA，但是理论上是应该各行均值归0的。

教程的均值归0（教程其实是针对的图像数据，而不适用于pca_2d中的数据）：

avg = mean(x, 1); % 对行求和再平均，使各列均值归0。
x = x - repmat(avg, size(x, 1), 1);

我的均值归0（PCA只有在x各行均值为0的前提下才成立）：

mean_x=mean(x,2); % 对列求和再平均，使各行均值归0。
x=x-repmat(mean_x,1,size(x,2));

 

(b) 此外，我认为应该一开始就把x各行均值归0化，因为协方差矩阵sigma可以用sigma = x * x' / size(x, 2)来估计其实是在均值为0的条件下才成立。

理由：X为n维随机变量向量，表示n个特征，x的m列中每一列都是X的一个样本，Σ=E（X-E[X]）(X-E[X])' ，只有均值做了归0处理，使x各行均值为0，即E[X]=0时，才有Σ=EX*X'= x* x' / size(x, 2)，第二个等号其实是对期望的样本估计，在size(x,2)=m=样本数 很大时才成立 。(X是随机变量构成的向量，x是n*m的矩阵)

 

2.Exercise:PCA_in_2D

(1)完全按照教程步骤的实现过程：pca_2d.m

这个版本所有结果图与教程答案完全一致

Step 1a: Finding the PCA basis

u = zeros(size(x, 1)); % You need to compute this %zeros(2)就生成2阶方阵
[u,s,v] = svd(x);

效果图：左为教程标准图，右为我的结果图。

Step 1b: Check xRot

xRot = zeros(size(x)); % You need to compute this
xRot = u'*x;

效果图：左为教程标准图，右为我的结果图。

 

Step 2: Dimension reduce and replot

k = 1; % Use k = 1 and project the data onto the first eigenbasis
xHat = zeros(size(x)); % You need to compute this
xHat=u*[xRot(1:k,:);zeros(size(x,1)-k,size(x,2))];

效果图：左为教程标准图，右为我的结果图。

 

Step 3: PCA Whitening

xPCAWhite = zeros(size(x)); % You need to compute this
mean_x=mean(x,2);
x=x-repmat(mean_x,1,size(x,2));
%------------------------------------------
% avg = mean(x, 1);     % 分别为每个图像块计算像素强度的均值
% x = x - repmat(avg, size(x, 1), 1);
%这是教程上介绍的对图像PCA的常用处理方法，但是对本例数据，应该求所有列的和再除以样本数求均值：
% mean_x=mean(x,2);
% x=x-repmat(mean_x,1,size(x,2));
%这样求均值，绘出来的xPCAWhite，xZCAWhite散点图和教程上的结果图一模一样，证明我的想法正确，
%另外，AndrewNGCS229PCA那章中也提到了这种均值规整化，所以不同数据的规整化方法并不一样
%--------------------------------------------
sigma = x * x' / size(x, 2)
[U,S,V] = svd(sigma);
xTilde = U(:,1:k)' * x; % reduced dimension representation of the data,                       
xPCAWhite = diag(1./sqrt(diag(S)+epsilon)) * U' * x;

效果图：左为教程标准图，右为我的结果图。

 

Step 4: ZCA Whitening

xZCAWhite = zeros(size(x)); % You need to compute this
xZCAWhite = U * diag(1./sqrt(diag(S) + epsilon)) * U' * x;

效果图：左为教程标准图，右为我的结果图。

(2)自己的实现过程：another_pca_2d.m(这部分具体分析看1.习惯性啰嗦几句)

clear,close all; %简化版计算
%%================================================================
%% Step 0: Load data
x = load('pcaData.txt','-ascii');
%---------------
mean_x=mean(x,2);
x=x-repmat(mean_x,1,size(x,2));
%---------------一开始就进行了对x各行均值归0化，因为只要有关于协方差的地方，其实都是需要均值为0这个条件的：
%Σ=E（X-E[X]）(X-E[X])' 在PCA的推导中（详见ng CS229 PCA那章）只有均值做了归0处理，即E[X]=0
%才能用sigma = x * x' / size(x, 2)估计协方差矩阵
figure(1);
subplot(231);scatter(x(1, :), x(2, :));
title('Raw data');
%%================================================================
%% Step 1a: Implement PCA to obtain U 
u = zeros(size(x, 1)); 
[u,s,v] = svd(x);
hold on
plot([0 u(1,1)], [0 u(2,1)]);
plot([0 u(1,2)], [0 u(2,2)]);
scatter(x(1, :), x(2, :));
hold off
%%================================================================
%% Step 1b: Compute xRot, the projection on to the eigenbasis
xRot = zeros(size(x)); 
xRot = u'*x;
%figure(2);
subplot(232);scatter(xRot(1, :), xRot(2, :));
title('xRot');
%%================================================================
%% Step 2: Reduce the number of dimensions from 2 to 1. 
k = 1;
xHat = zeros(size(x)); 
xHat=u*[xRot(1:k,:);zeros(size(x,1)-k,size(x,2))];
%figure(3);
subplot(233);scatter(xHat(1, :), xHat(2, :));
title('xHat');
%%================================================================
%上面几步和原版一样
%================================================================
%% Step 3: PCA Whitening
%PCA就是求x的协方差矩阵sigma的特征向量和特征值
%可以直接用前面的s来计算而不用计算sigma(CS229中PCA那章提出svd正是这个目的)，s里的值是x的奇异值，同时也是x的协方差矩阵的特征值乘以m（样本数）再开根号
%U和u也是一样的，具体参考前文 1.习惯性啰嗦几句 ，这种计算方法省去了很多步骤，和原结果基本一样。
epsilon = 1e-5;
xPCAWhite = zeros(size(x)); 
lambda=(diag(s).^2)/size(x,2);%协方差矩阵的特征值，用s平方再除以m(样本数)即可求得
xPCAWhite = diag(1./sqrt(lambda+epsilon)) * u' * x;
%figure(4);
subplot(234);scatter(xPCAWhite(1, :), xPCAWhite(2, :));
title('xPCAWhite');
%%================================================================
%% Step 4: ZCA Whitening
xZCAWhite = zeros(size(x)); 
xZCAWhite = u * diag(1./sqrt(lambda+epsilon)) * u' * x;
%figure(5);
subplot(235);scatter(xZCAWhite(1, :), xZCAWhite(2, :));
title('xZCAWhite');

效果图如下：


结果分析：

 关于这个版本我的各种改动，请看前文1.习惯性啰嗦几句

（1）可以看到前3幅图和标准答案有很小的差别，这是因为我在一开始就对x各行进行了均值归0化，而教程应该是到xPCAWhite这步才均值归0的（参见我的第一种实现），而这次练习用的数据x（2*45）的各行均值接近0但不等于0，所以均值归0和不归0有一些差别，这就是前3幅图有一点差别的原因。

（2）后俩幅图和教程答案完全一样了，证明我的算法应该是正确的，大家有不同意见请和我讨论。

 

3.PCA_and_Whitening

Step 0b: Zero mean the data

（1）这次的数据各行均值又接近0，所以不用使各行均值归0也可以使用PCA

（2）数据为自然图像，所以需要根据教程中“对图像数据应用PCA算法”那小节，对数据各列均值归0化。

avg=mean(x,1);
x=x-repmat(avg,size(x,1),1);

Step 1a: Implement PCA

（1）按照教程的步骤：

xRot = zeros(size(x)); % You need to compute this
sigma = x * x' / size(x, 2);
[U,S,V]=svd(sigma);
xRot=U'*x;

（2）使用的1.(1)中提到的用svd绕开求sigma的方法：

xRot = zeros(size(x)); % You need to compute this
[u,s,v]=svd(x);
xRot=u'*x;

Step 1b: Check covariance

covar = zeros(size(x, 1)); % You need to compute this
covar = xRot * xRot' / size(xRot, 2);

可以看到两种方法求得的协方差矩阵相同：除对角线元素外，其余均非常小，基本为0。（为了与教程保持一致性，并且pca_2d练习中也给出了两种方法的完整实现，所以后面步骤我都按照教程的方法来实现了）

 

注意：（1）多次运行，结果不一样，原因是样本是由sampleIMAGESRAW()函数随机生成的，x每次运行都不一样。

（2）凡是用这种 xRot * xRot' / size(xRot, 2) 方式来计算协方差的，都是在假设数据x各行均值为0的情况下才成立（具体见1.(2).(b)）。练习中数据集的各行均值恰好接近0，所以可以勉强这样算，大家可以用matlab的cov(xRot')来算算，会发现运气不好的时候（因为数据集是随机生成的），不少对角线以外的元素值其实挺大的，并不非常接近0。      

 

Step 2: Find number of components to retain

k = 0; % Set k accordingly
diag_sum=trace(S);
sum=0;
rate=0.99;
while((sum/diag_sum)<rate)
    k=k+1;
    sum=sum+S(k,k);
end

可以修改rate为0.9。

Step 3: PCA with dimension reduction

xHat = zeros(size(x));  % You need to compute this
xTilde = U(:,1:k)' * x;
xHat(1:k,:)=xTilde;
xHat=U*xHat;

结果图如下，左为原始数据；中间rate=99%，k=117；右边rate=90%，k=43；

 

Step 4a: Implement PCA with whitening and regularization

epsilon = 0.1;
xPCAWhite = zeros(size(x));
xPCAWhite = diag(1./sqrt(diag(S) + epsilon)) * U' * x;

可以修改epsilon为0，即不进行正则化。

Step 4b: Check covariance

covar = zeros(size(xPCAWhite, 1));
covar = xPCAWhite * xPCAWhite' / size(xPCAWhite, 2);

使用了epsilon正则化时，协方差矩阵对角线元素不是1，而是会比1小，而且沿着对角线向右下角方向，元素会越来越小到（0.1左右）。

Step 5: ZCA whitening

xZCAWhite=U * diag(1./sqrt(diag(S) + epsilon)) * U' * x;

结果图，左为ZCA白化结果，右为原图：