De-PCA? NO!
Introduction
今天在使用PCA的时候，突然想到一个问题：根据PCA得到的系数矩阵和降维后的主成分的值可以还原到原维度的值吗？在思考这个问题的时候，对之前学习到的知识有了进一步的认识，很有意思，记录一下。
根据PCA得到的系数矩阵和降维后的主成分的值可以还原到原维度的值吗？从几何的角度看，这种想法是想把一个低维空间的向量转化为高维空间的向量，可以预想到在这个过程中，信息的损失似乎是必然的；从矩阵论的角度看，这个问题等价于——从PCA得到的系数矩阵出发，是否能够得到一个线性变换，将低维度的向量转换为高纬度的向量？如果不能，能不能在某种意义下得到一个最好的线性变换？另外，如果不限制还原操作为线性变换呢？
pca Function in MATLAB
PCA很重要的一个用处就是数据降维，将一个高维的特征向量变换成一个低维的特征向量，以MATLAB提供的hald数据集和pca函数为例进行展示。
变量hald中保存的数据：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
hald =
7.0000 26.0000 6.0000 60.0000 78.5000
1.0000 29.0000 15.0000 52.0000 74.3000
11.0000 56.0000 8.0000 20.0000 104.3000
11.0000 31.0000 8.0000 47.0000 87.6000
7.0000 52.0000 6.0000 33.0000 95.9000
11.0000 55.0000 9.0000 22.0000 109.2000
3.0000 71.0000 17.0000 6.0000 102.7000
1.0000 31.0000 22.0000 44.0000 72.5000
2.0000 54.0000 18.0000 22.0000 93.1000
21.0000 47.0000 4.0000 26.0000 115.9000
1.0000 40.0000 23.0000 34.0000 83.8000
11.0000 66.0000 9.0000 12.0000 113.3000
10.0000 68.0000 8.0000 12.0000 109.4000
一共有13条数据，每条数据有5个特征量。
之后，对其进行PCA：
1
[coeff, score, latent, tsquared, explained, mu] = pca(hald);
pca函数的输出：
-
变量
coeff，表示the principal component coefficients, also known as loadings，对于具有$p$个特征的数据，得到的主成分系数矩阵的形状就是$p\times p$。1 2 3 4 5 6
coeff = -0.0951 -0.4881 0.2635 0.5967 0.5721 -0.5735 0.2673 -0.5939 -0.0807 0.4904 0.0611 0.5376 0.6228 -0.2904 0.4848 0.6179 -0.3368 -0.3239 -0.4550 0.4391 -0.5259 -0.5366 0.2917 -0.5883 -0.0659
-
变量
socre，表示the principal component scores，在principal component space中表示数据hald。socre的行对应数据的数量，列对应各成分。1 2 3 4 5 6 7 8 9 10 11 12 13 14
score = 39.8331 -9.8213 -5.2115 -0.5073 0.3642 36.4980 3.6952 -1.6029 -0.8316 -0.4699 -15.9135 -3.0512 -0.2464 1.9006 -0.9318 23.8891 -9.8666 0.9838 1.4567 -0.2345 -0.9111 -3.1142 -6.8314 -0.5559 0.1112 -16.6198 -6.0836 1.7519 -2.1020 -0.3823 -31.0147 15.2749 -1.5894 0.6155 0.1687 31.7817 11.6528 3.6351 1.6734 0.5103 -6.1745 11.5187 0.8831 -0.5327 -0.5966 -14.3400 -20.7331 6.6827 0.1999 0.3059 14.5601 11.9010 5.4483 -1.4408 0.2723 -31.2632 -1.9752 -0.3458 -0.8512 0.3506 -30.3253 0.6025 -3.5575 0.9757 0.5318
-
变量
latent，表示the principal component variances，它是数据hald的协方差矩阵的的特征值；1 2 3 4 5 6
latent = 695.3761 111.7525 15.3443 1.5555 0.2211
-
变量
tsquared，表示数据hald的Hotelling’s T-squared statistic(霍特林$\mathrm{T}^2$统计量)1 2 3 4 5 6 7 8 9 10 11 12 13 14
tsquared = 5.6803 3.6484 6.7002 3.3676 3.3839 4.4300 4.0080 6.5067 3.0849 7.5016 5.1768 2.4701 4.0413
-
变量
explained，the percentage of the total variance explained by each principal component，1 2 3 4 5 6
explained = 84.3648 13.5581 1.8616 0.1887 0.0268
每个成分的贡献度就是根据这个变量来看的，比如：
1 2 3 4 5 6
ans = 84.3648 97.9229 99.7845 99.9732 100.0000
-
变量
mu，the estimated mean of each variable，其实就是样本均值。注：这里没有设置任何的Name-Value，如果进行了其他设置，可能就不是简单的样本均值了。
1 2
mu = 7.4615 48.1538 11.7692 30.0000 95.4231
From point of view of orthogonal matrix
在最开始思考开头提出的问题时，很朴素的想法就是$13\times5$的数据矩阵$A$乘$5\times2$的主成分系数矩阵$B$，就可以降到2维，得到$13\times2$的降维后的数据矩阵$C$，那么将数据矩阵乘$B^T$，那么就可以将数据还原到5维，但是结果并非如此。为了验证这个想法，定义一个循环将原始数据hald进行PCA后，选取前idx个主成分降维，之后再进行还原为数据矩阵$D$，观察矩阵$D-A$二范数的变化：
1
2
3
4
5
6
7
8
A = hald;
[coeff, score, latent, tsquared, explained, mu] = pca(A);
for idx = 1:width(A)
B = coeff(:, 1:idx);
C = A*B;
D = C*B';
disp(norm(D-A, 2))
end
1
2
3
4
5
342.6863
299.8675
299.8373
145.9264
1.7229e-13
可以看到，随着主成分数量的增加，两者之间的差异是原来越小的，当选取的主成分的个数和原始特征的数量一致时，差异就几乎为零，如果考虑到数值计算的误差，它们就是相等的。
注：如果在PCA之前进行了标准化，那么在PCA之后是不需要反标准化的：
如果进行了反标准化，那么结果就会相差极大：
我们期望得到的结果是：选取不同的主成分个数，都有norm(D-A, 2)值为零，那么就回答了最开始提出的问题，并且还原的步骤也很简单。但是很明显，当选取的成分数小于原始特征数量时，结果是否定的。
其实可以从数学上来解释这样的现象。仔细地考虑上面的这个步骤，其实我们是想有：
\[A\times B\times B^T=A\notag\]即：
\[BB^T=E\label{eq1}\]也就是从几何角度看，先降维后还原的操作等于什么也没有做。
而根据正交矩阵的定义：
正交矩阵(Orthogonal matrix)
如果$AA^T=E$或者$A^TA=E$，则$n$阶实矩阵$A$为正交矩阵。（注意，这里的矩阵$A$一定是方阵，方阵才可能是正交矩阵）
设$B$为PCA的主成分系数矩阵，形状为$m\times n$，$m$为降维前的特征数量，$n$为降维后的特征数量。假设有$BB^T=E$成立，则根据定义，$B$为正交矩阵。当$m>n$时，矩阵$B$不是方阵，因此假设不成立，即此时不可能存在关系$BB^T=E$成立；当$m=n$时，本问题所对应的主成分系数矩阵coeff满足：
1
2
3
4
5
6
7
>> coeff*coeff'
ans =
1.0000 0 0.0000 0.0000 -0.0000
0 1.0000 0.0000 0.0000 -0.0000
0.0000 0.0000 1.0000 -0.0000 0.0000
0.0000 0.0000 -0.0000 1.0000 0.0000
-0.0000 -0.0000 0.0000 0.0000 1.0000
有$BB^T=E$成立，也就解释了上述代码运行的结果。
From point of view of LSM
针对上面这个问题，我们所采用的方法——利用保存的主成分系数矩阵——还原出的矩阵矩阵与原始数据矩阵之差的二范数分别为：
1
2
3
4
5
342.6863
299.8675
299.8373
145.9264
1.7229e-13
这种还原方法得到的结果是最优的吗？大概率不是的。
假设原始数据矩阵为$A_{p\times m}$，降维后的数据矩阵为$B_{p\times n}$(假设不考虑前后维度相同的情况，即$m>n$)，如果不利用保存的主成分系数矩阵，那么我们就是想找到这样一个变换$X$，使得：
\[B_{p\times n}X_{n\times m}=A_{p\times m}\notag\]其中，矩阵$B$和矩阵$A$是已知的。这实际上就是超定方程组的求解问题，我们可以求得其最小二乘解为：
\[\hat{X}_{n\times m}=(B^TB)^{-1}B^TA\notag\]使用MATLAB实现一下：
1
2
3
4
5
6
7
8
9
10
11
12
clc, clear, close all
load hald
A = hald;
[coeff, score, latent, tsquared, explained, mu] = pca(A);
for idx = 1:width(A)
B = A*coeff(:, 1:idx);
X = inv(B'*B)*B'*A;
D = B*X;
disp(norm(D-A, 2))
end
1
2
3
4
5
135.8125
58.0986
57.9865
2.8992
3.1370e-11
可以看到，后面这种不利用保存的主成分系数矩阵，而采用最小二乘解的方法还原的数据矩阵和原数据矩阵之间的差异是更小的！！！
What about nonlinear transformation?
再进一步，既然我们已经放弃了使用主成分系数矩阵的方法，我们可不可以采用非线性变换进行还原呢？
采用非线性变换大概有两大类方法。一类是假设非线性变换的形式，比如说是多项式形式(多项式形式的向量函数)等等，之后根据数据拟合参数；一类是采用机器学习模型，我们所面临的其实是一个多输入(低维特征向量)多输出(还原的高维特征向量)的回归问题。
前一种方式应该属于非线性映射和函数空间的数学范畴，目前我不是很了解其中的知识，但是印象比较深刻的是，在SVM中经常看到说，数据点在低维空间是线性不可分的，可以使用核函数(kernel function)将低维向量映射到高维空间中，就可能做到线性可分。
后一种使用机器学习模型的方法是个比较现代的方法，学术界应该对这个问题已经有了很多的研究。但是可以预想，尽管使用机器学习模型得到的还原效果会比LSM要好一些，但是也无法完全还原到原来的维度。
比LSM效果好的原因在于：尽管我们在训练模型的时候同样会采用最小二乘作为损失函数，但是机器学习模型是在寻找一种非线性映射(尽管没有显式表达式)，其效果会比线性映射好一些，但最终的效果还是依赖于训练集的质量。
Conclusion
综上，本文所想要得出的结论是：原始特征数据经过PCA降维后，我们并不能使用保存的PCA的系数矩阵实现完美的还原，并且最好的(在最小二乘的损失函数下)线性还原是使用超定方程的最小二乘法进行求解。如果不将还原操作限制在线性的情况，则可以采用其他非线性的形式，比如拟合非线性向量函数、采用多输入多输出的机器学习回归模型。但无论如何，正如本文一开始所预期的，最终都不可能实现完美的还原，因为我们在降维的过程中损失了信息。
Reference
[1] pca - MathWorks.