本文没有从头讲PCA的思想, 建议有一定基础后再阅读.
之后若有时间, 争取补充源码每一步对应的公式(如果有人需要, 可以评论区告诉我, 我会提前安排)
目录
Sklearn的decomposition.PCA类源码注释
PCA逆转换实验
创建pcaInverseDemo.py进行尝试
先试一下调包, 然后对比手刻的代码
sklearn的逆转换
建立简单矩阵, PCA转换, 输出值
| import numpy as np from sklearn import decomposition # 建立简单矩阵 X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]]) # 将含有2个特征的数据经过PCA压缩为1个特征 pca = decomposition.PCA(n_components=1) pca.fit(X) X_pca = pca.transform(X) print("X_pca:\n",X_pca) |
| X_pca: [[ 1.38340578] [ 2.22189802] [ 3.6053038 ] [-1.38340578] [-2.22189802] [-3.6053038 ]] |
再逆转换, 输出值
| X_origin=pca.inverse_transform(X_pca) |
| X_origin: [[-1.15997501 -0.75383654] [-1.86304424 -1.21074232] [-3.02301925 -1.96457886] [ 1.15997501 0.75383654] [ 1.86304424 1.21074232] [ 3.02301925 1.96457886]] |
如果维度降低, 则会损失信息, 如果进行PCA时维度不变, 则逆转换后值与原来相同
手刻PCA以及逆转换
PCA:
| def pca(X, k): # k is the components you want # mean of each feature n_samples, n_features = X.shape mean = np.array([np.mean(X[:, i]) for i in range(n_features)]) # normalization norm_X = X - mean # scatter matrix scatter_matrix = np.dot(np.transpose(norm_X), norm_X) # Calculate the eigenvectors and eigenvalues eig_val, eig_vec = np.linalg.eig(scatter_matrix) eig_pairs = [(np.abs(eig_val[i]), eig_vec[:, i]) for i in range(n_features)] # sort eig_vec based on eig_val from highest to lowest eig_pairs = [((np.abs(eig_val[i])), eig_vec[:, i]) for i in range(len(eig_vec))] eig_pairs.sort(key=lambda x: x[0], reverse=True) # select the top k eig_vec feature = np.array([ele[1] for ele in eig_pairs[:k]]) # get new data new_data = np.dot(norm_X, np.transpose(feature)) return new_data X_pca = pca(X, 1).real#.real取得复数的实部 print("X_pca:\n",X_pca) |
| X_pca: [[-1.38340578] [-2.22189802] [-3.6053038 ] [ 1.38340578] [ 2.22189802] [ 3.6053038 ]] |
和刚刚结果绝对值相同, 正负相反, 但是和原数据的正负还是保持一致的 (原因已找到, sklearn用的是工业界技巧, 统一使最大主成分为正)
逆转换:
| def pca(X, k): # k is the components you want # mean of each feature n_samples, n_features = X.shape mean = np.array([np.mean(X[:, i]) for i in range(n_features)]) # normalization norm_X = X - mean # scatter matrix scatter_matrix = np.dot(np.transpose(norm_X), norm_X) # Calculate the eigenvectors and eigenvalues eig_val, eig_vec = np.linalg.eig(scatter_matrix) eig_pairs = [(np.abs(eig_val[i]), eig_vec[:, i]) for i in range(n_features)] # sort eig_vec based on eig_val from highest to lowest eig_pairs = [((np.abs(eig_val[i])), eig_vec[:, i]) for i in range(len(eig_vec))] eig_pairs.sort(key=lambda x: x[0], reverse=True) # select the top k eig_vec feature = np.array([ele[1] for ele in eig_pairs[:k]]) # get new data new_data = np.dot(norm_X, np.transpose(feature)) return new_data,feature,mean # X_pca = pca(X, 1)#可行 X_pca,X_components,X_mean = pca(X, 1) print("X_pca:\n",X_pca) X_origin=np.dot(X_pca, X_components)+X_mean print("X_origin:\n",X_origin) |
| X_pca: [[-1.38340578] [-2.22189802] [-3.6053038 ] [ 1.38340578] [ 2.22189802] [ 3.6053038 ]] X_origin: [[-1.15997501 -0.75383654] [-1.86304424 -1.21074232] [-3.02301925 -1.96457886] [ 1.15997501 0.75383654] [ 1.86304424 1.21074232] [ 3.02301925 1.96457886]] |
结果与调包相同, 基本没问题
Sklearn用户手册
2.5矩阵分解问题Decomposing signals in components (matrix factorization problems)
主成分分析 （PCA）
标准的PCA和概率解释
Exact PCA and probabilistic interpretation
PCA 用于分解多元数据集，通过一组连续的正交分量, 这些分量解释了方差的最大量。在 scikit-learn 中，PCA 的实现是一个transformer 对象，该对象在其fit方法中学习n个组件，并且可以用于将新数据投影到这些组件上
在应用 SVD 之前，PCA 对于每个特征中心化，但不缩放输入数据。可选参数whiten=True可以将数据投影到奇异空间(singular space)，同时将每个分量(component)缩放到单位方差(unit variance)。如果下游模型对信号的各向同性(isotropy)做出强有力的假设，这通常很有用：例如，对于具有 RBF 核的支持向量机和 K-Means 聚类算法。
| 中心化即变量减去它的均值, 数据中心平移到(0,0) 参考: https://blog.csdn.net/qq_36523839/article/details/82919412 中心化后的数据对向量也容易描述，因为是以原点为基准的 PCA需要找出矩阵的特征向量，也就是主成分（PC）。比如说找到的第一个特征向量是a = [1, 2]，a在坐标平面上就是从原点出发到点（1，2）的一个向量。 如果没有对数据做中心化，那算出来的第一主成分的方向对于数据的“描述”就没有这样的意义了.
|
PCA 对象还提供了 PCA 的概率解释(probabilistic interpretation)，(即评估降维效果的好坏的评分).
该解释可以根据其方差量(amount of variance)提供数据的可能性(likelihood)。因此，它实现了可用于交叉验证的评分方法. 这样可以用来选择模型, 如维度的自动选择.
Sklearn的decomposition.PCA类源码注释
Principal component analysis (PCA)
使用数据的奇异值分解将其投影到较低维空间的线性降维。
在应用 SVD 之前，输入数据居中但未针对每个特征进行缩放。它使用 LAPACK 实现通过 Halko et al. 2009 提出的方法提取完整 SVD 或随机截断 SVD，具体取决于输入数据的形状和要提取的组件数量。
它也可以使用截断的 scipy.sparse.linalg ARPACK 实现SVD。
属性Attributes
components_ :
形状为(n_components, n_features) 的ndarray
特征空间中的主轴(Principal axes)，表示数据中最大方差的方向。等效地，经过中心化的输入数据的右奇异向量(right singular vectors)与其特征向量平行。组件通过“explained_variance_”排序。
explain_variance_ :
形状为 (n_components,) 的ndarray
通过每个选定组件解释的方差量
The amount of variance explained by each of the selected components.
方差估计使用“n_samples - 1”自由度。
等于X 的协方差矩阵的n_components最大特征值
singular_values_：
ndarray, 形状（n_components，）对应于每个选定组件的奇异值。奇异值等于低维空间中“n_components”变量的 2 范数。
mean_ : ndarray of shape (n_features,)
每个特征的经验平均值，从训练集估计。等于`X.mean(axis=0)`。
noise_variance_ : float 遵循 Tipping 和 Bishop 1999 的概率 PCA 模型的估计噪声协方差。参见 C. Bishop 的“模式识别和机器学习”，第 12.2.1 页。 574 或 http://www.miketipping.com/papers/met-mppca.pdf。需要计算估计的数据协方差和分数样本。等于 X 的协方差矩阵的 (min(n_features, n_samples) - n_components) 个最小特征值的平均值。
还提到一些参考资料
See Also
--------
KernelPCA : Kernel Principal Component Analysis.
SparsePCA : Sparse Principal Component Analysis.
TruncatedSVD : Dimensionality reduction using truncated SVD.
IncrementalPCA : Incremental Principal Component Analysis.
For n_components == 'mle', this class uses the method from:
`Minka, T. P.. "Automatic choice of dimensionality for PCA".
In NIPS, pp. 598-604 <https://tminka.github.io/papers/pca/minka-pca.pdf>`_
Implements the probabilistic PCA model from:
`Tipping, M. E., and Bishop, C. M. (1999). "Probabilistic principal
component analysis". Journal of the Royal Statistical Society:
Series B (Statistical Methodology), 61(3), 611-622.
<http://www.miketipping.com/papers/met-mppca.pdf>`_
via the score and score_samples methods.
手刻代码与包内源码对照
Sklearn代码
使用时主要用到fit(X), transform(X)两个函数
_fit_full
通过计算 X 上的full SVD 来拟合模型
| def _fit_full(self, X, n_components): n_samples, n_features = X.shape |
数据中心化
| self.mean_ = np.mean(X, axis=0) X -= self.mean_ |
SVD分解
| U, S, Vt = linalg.svd(X, full_matrices=False) U, Vt = svd_flip(U, Vt)# 翻转特征向量的符号以强制执行确定性输出 components_ = Vt |
其中svd_flip
| max_abs_cols = np.argmax(np.abs(u), axis=0) signs = np.sign(u[max_abs_cols, range(u.shape[1])]) u *= signs v *= signs[:, np.newaxis] |
svd_flip保证svd分解结果的一致性(这就是之前实验得到的结果(PCA的正常解法)与sklearn(工业界常用的SVD解法)不同的原因)
| 参考: https://zhuanlan.zhihu.com/p/85701151 SVD奇异值分解的结果是唯一的，但是分解出来的U矩阵和V矩阵的正负可以不是唯一，只要保证它们乘起来是一致的就行。因此，sklearn为了保证svd分解结果的一致性，它们的方案是：保证U矩阵的每一行(u_i)中，绝对值最大的元素一定是正数，否则将u_i转成-u_i,并将相应的v_i转成-v_i已保证结果的一致。 这又是数学与工程的问题了。在数学上，几种结果都是正确的。但是在工程上，有个很重要的特性叫幂等性(Idempotence)。 高质量的服务，一次请求和多次请求，其副作用（结果）应当是一致的。Scikit Learn正是通过svd_flip这个函数，把一个数学上并不幂等的操作，转化成了幂等的服务(讲究人儿) |
截取所需组件
| self.components_ = components_[:n_components] |
| return U, S, Vt |
手刻代码
与sklearn不一致的地方主要在于:
除了翻转符号部分手刻代码没有(上节提过, 依然符合正常PCA)
还有svd分解部分, 这里没有调用linalg的包, 而是通过scatter matrix做特征分解, 应该是属于PCA两类求解方法中的第一类(协方差的特征向量), 而sklearn中使用的是第二类求解方法(奇异值分解). 这两类方法可参考下篇笔记, 下篇看懂之后也尝试做代码与公式对照.
| # scatter matrix scatter_matrix = np.dot(np.transpose(norm_X), norm_X) # Calculate the eigenvectors and eigenvalues eig_val, eig_vec = np.linalg.eig(scatter_matrix) |
其中散布矩阵scatter_matrix等价于协方差矩阵
散布矩阵可参考:
https://blog.csdn.net/guyuealian/article/details/68922981
最后将原特征X向Y线性转换部分的手刻代码有两种实现方法:
| V = np.array(V_T).T#V从特征向量取出来已经经过一次转置 Y = np.matmul(Vk.T, X) #计算主成分矩阵，将m*n的样本矩阵X转换成k*n的样本主成分矩阵 |
| new_data = np.dot(norm_X, np.transpose(feature)) |
第一种更接近统计学习方法中给出的公式. 第二种不需要进行多次转置.
此处的Vk(每列为一个主成分)
$$ V_k=\left[\begin{array}{llll} v_{1} & v_{2} & \cdots & v_{k} \end{array}\right]=\left[\begin{array}{cccc} v_{11} & v_{12} & \cdots & v_{1 k} \\ v_{21} & v_{22} & \cdots & v_{2 k} \\ \vdots & \vdots & & \vdots \\ v_{m 1} & v_{m 2} & \cdots & v_{m k} \end{array}\right] $$
此处的feature(每行为一个主成分)
$$ A_k=\left[\begin{array}{l} a_{1}^T \\ a_{2}^T \\ \cdots \\ a_{k}^T \end{array}\right] $$
两种方法的矩阵相乘是等价的
$$ V_k^TX= \left[\begin{array}{llll} v_{1} & v_{2} & \cdots & v_{k} \end{array}\right]\left[\begin{array}{l} x_{1}^T \\ x_{2}^T \\ \cdots \\ x_{k}^T \end{array}\right]\\ \Leftrightarrow XA_k^T=\left[\begin{array}{l} x_{1}^T \\ x_{2}^T \\ \cdots \\ x_{k}^T \end{array}\right]\left[\begin{array}{llll} a_{1} & a_{2} & \cdots & a_{k} \end{array}\right]$$
特征重要性在原始空间含义
标准化
真正使用PCA之前, 还应该进行标准化. 一方面是为了解释含义无误, 一方面是PCA涉及到距离计算.
| 需要标准化的原因参考: https://blog.csdn.net/weixin_36604953/article/details/102652160 归一化的最通用模式Normalization(或称rescaling) Mean normalization, 范围[-1,1]: $$ X_{new}=\frac{X_{i}-X_{min}}{X_{max}-X_{min}} $$ 标准化(Standardization)，也叫标准差标准化 $$ X_{new}=\frac{X_{i}-\mu}{\sigma} $$ 归一化与标准化的作用: 1.统计建模中，如回归模型，自变量X的量纲不一致导致了回归系数无法直接解释或者错误解读；需要将X都处理到统一量纲下，这样才可比； 2.机器学习任务和统计学任务中有很多地方要用到“距离”的计算，比如PCA，比如KNN，比如kmeans等等，假使算欧式距离，不同维度量纲不同可能会导致距离的计算依赖于量纲较大的那些特征而得到不合理的结果； 3.参数估计时使用梯度下降，在使用梯度下降的方法求解最优化问题时， 归一化/标准化后可以加快梯度下降的求解速度，即提升模型的收敛速度。 归一化与标准化的选择 标准化是ML中更通用的手段，如果你无从下手，可以直接使用标准化；如果数据不为稳定，存在极端的最大最小值，不要用归一化。在分类、聚类算法中，需要使用距离来度量相似性的时候、或者使用PCA技术进行降维的时候，标准化表现更好；在不涉及距离度量、协方差计算的时候，可以使用归一化方法。 PCA中标准化表现更好的原因可以参考: https://blog.csdn.net/young951023/article/details/78389445 线性变换后，其协方差产生了倍数值的缩放，因此这种方式无法消除量纲对方差、协方差的影响，对PCA分析影响巨大；同时，由于量纲的存在，使用不同的量纲、距离的计算结果会不同。 标准化方式中，新的数据由于对方差进行了归一化，这时候每个维度的量纲其实已经等价了，每个维度都服从均值为0、方差1的正态分布，在计算距离的时候，每个维度都是去量纲化的，避免了不同量纲的选取对距离计算产生的巨大影响。 |
每个维度都服从均值为0、方差1的正态分布, PCA还如何根据方差选取特征呢? PCA最优化的目标是投影之后描述特征之间的散布(斜跨各特征)的协方差越大越好, 不是原空间(一个特征垂直方向)的同一个特征内的方差.
标准化后数据拟合的变换矩阵, 用来逆转换降维后特征的重要性到原空间, 应该是可以解释原空间特征重要度的.
PCA论文
PCA理论参考:
https://www.cs.princeton.edu/picasso/mats/PCA-Tutorial-Intuition_jp.pdf
文章翻译见下篇笔记, 这里简单摘取总结.
目标总结
找正交矩阵 P，其中 Y = PX 使得 \( \mathbf{S}_{\mathbf{Y}} \equiv \frac{1}{n-1} \mathbf{Y} \mathbf{Y}^{T} \)对角化。 P 的行是 X 的主成分。
总结 PCA 的结果
1. 输入数据X 的主成分(PC)是 XXT 的特征向量；或变换矩阵P 的行
2. 协方差矩阵S_Y 的对角线第 i 个值是 X 沿 pi 的方差
在实践中执行 PCA 非常简单:
1. 将数据集组织为 m × n 矩阵，其中 m 是测量类型的数量，n 是试验的数量。
2. 减去每个测量类型或第 xi 行的平均值。
3. 计算协方差的 SVD 或特征向量。
在多个文献领域中，许多作者将单个测量类型 xi 称为源(sources)。投影到主成分 Y = PX 中的数据称为信号(signals)，因为投影数据可能代表“真实”的潜在概率分布。
重构误差
分析PCA映射回原空间的时候，乘转换矩阵得到的结果和原始数据有误差
转换原空间时候的误差似乎能表达成这个问题:
https://www.zhihu.com/question/428735076
但我还是没想通数学上的原因.
(简化问题, 不考虑中心化)
转换时是
| X_pca = np.dot(X, np.transpose(X_components)) |
$$ \mathbf Y=\mathbf P\mathbf X $$
逆转换时是
| X_origin=np.dot(X_pca, X_components) |
$$ \mathbf X=\mathbf P^{T}\mathbf Y $$
得到结果X_origin就是和X不一样
此时得到的转换矩阵已经不是真正的P矩阵了, 而是P矩阵的前k行, 剩下组件的被抛弃了
(参考矩阵论的换基函数, 当P为可逆矩阵, 也就是方阵的时候, 才能进行完整的坐标系变换, 且各个坐标系的坐标轴数量一致.)
此处的P如果完整的话, 根据正交矩阵的性质, P的转置与P的逆等价, 所以可以用P的转置的前几个主成分代替P的逆的前几个主成分. 但如果单考虑这个不完整的转换矩阵P, 它根本不是方阵, 也不可逆.
可以尝试根据《统计学习方法》314页的例16.1手算
每个组件内的各个值含义:
就是做坐标系旋转时的线性变换系数.
2111
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
