PCA原理及代码实现

Posted on Edited on

0 PCA的过程

1. 数据准备与中心化

原始表达矩阵 A (n×p)：n个样本，p个基因
对每个基因（特征）进行中心化：B = A - A.mean(axis=0)
中心化后矩阵 B：每列均值为0

2. 奇异值分解（SVD）

这里用的是紧凑型SVD，为零的奇异值直接省略了。
对中心化矩阵 B 进行SVD分解：B = UΣVᵀ
其中：
U：左奇异向量矩阵 (n×r)
Σ：奇异值对角矩阵 (r×r)
Vᵀ：右奇异向量矩阵转置 (r×p)

3. 主成分得分计算

主成分方向：V 的列向量（特征向量）
样本在主成分空间的坐标：scores = BV
等价计算：scores = UΣ
每个样本在PC1, PC2, …上的坐标即为PCA降维结果

4. 解释方差

各主成分解释方差比例：explained_variance_ratio = s² / sum(s²)
其中 s 为奇异值向量

1 一些问题

1 BV的意义是什么？

我们如果只考虑 bi*vi 的话，bi为p维空间中的一点（是一个p维的向量），vi是一个p维的特征向量。现在我们考虑吧一下两个向量点乘的几何意义：a·b = |a|·|b|·cosx,如果b为单位向量的话，a/b向量点乘的几何意义为a向量在b向量方向上的投影（坐标）。

所以BV可以理解成，样本点在PCi（特征向量）上的投影（坐标）。我们希望每一个样本点在PCi上的投影坐标的方差尽量大，考虑一个极端情况，如果所有的样本点在PC1上的坐标相差很大，实际上我们完全可以将PC1作为一个轴，并将样本散布在上面，也就是实现了降维。

2 为什么vi的方向是B投影坐标方差最大的方向？

下面来进行数学推导：

2 PCA的代码实现及其应用

手动PCA

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris

# 1. 加载Iris数据集
iris = load_iris()
X = iris.data  
y = iris.target
target_names = iris.target_names

print("=== Iris Dataset Positive Control ===")
print(f"Data shape: {X.shape}")
print(f"Samples: {X.shape[0]}, Features: {X.shape[1]}")
print(f"Classes: {target_names}")
print(f"Feature names: {iris.feature_names}")

# 2. 数据处理
iris_df = pd.DataFrame(X, columns=iris.feature_names)
iris_df['target'] = y
iris_df['target_name'] = [target_names[i] for i in y]

print("\n=== Manual PCA (Your Method) ===")

# 行为样本 ，列为特征
A = iris_df[iris.feature_names]  # (150, 4)
print(f"Final A shape: {A.shape}")
print(f"the rank of A is {np.linalg.matrix_rank(A)}")

# 中心化
B = A - A.mean()
print(f"Centered B shape: {B.shape}")
print(f"the rank of B is {np.linalg.matrix_rank(B)}")

# 3. 进行SVD分解
U, s, Vt = np.linalg.svd(B, full_matrices=False)
scores_manual = U @ np.diag(s)

print(f"SVD Results:")
print(f"U shape: {U.shape}, s length: {len(s)}, Vt shape: {Vt.shape}")
print(f"Principal component scores shape: {scores_manual.shape}")

# 4. 解释方差计算
explained_variance = s**2 / (B.shape[0] - 1)
total_variance = np.sum(explained_variance)
explained_variance_ratio = explained_variance / total_variance

print(f"\nExplained Variance Ratio (Manual PCA):")
for i, ratio in enumerate(explained_variance_ratio):
    print(f"PC{i+1}: {ratio:.3%}")

# 5. 绘制PCA结果
plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
colors = ['red', 'blue', 'green']
for i, target in enumerate(target_names):
    idx = y == i
    plt.scatter(scores_manual[idx, 0], scores_manual[idx, 1], 
               alpha=0.7, label=target, s=50)

plt.title(f"Manual PCA: Iris Dataset\n(Should show clear class separation)")
plt.xlabel(f"PC1 ({explained_variance_ratio[0]:.2%} variance)")
plt.ylabel(f"PC2 ({explained_variance_ratio[1]:.2%} variance)")
plt.legend()
plt.grid(False)



# 6. 显示主成分中每一个特征的贡献（系数）,系数也有正负
pc_columns = [f"PC{i+1}" for i in range(Vt.shape[0])]
pc_df = pd.DataFrame(Vt.T,  
                    index=iris.feature_names,
                    columns=pc_columns) 
print(pc_df)

# 可视化前两个主成分
plt.subplot(1, 2, 1)
pc_df[['PC1', 'PC2']].plot(kind='bar', figsize=(5, 5))
plt.title("Feature Contributions to PC1 and PC2 (Manual PCA)")
plt.xlabel("Features")
plt.ylabel("Contribution")
plt.grid(True, axis='y', alpha=0.3)
plt.tight_layout()
plt.xticks(rotation=45, ha='right') 
plt.show()

手动PCA运行结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

=== Iris Dataset Positive Control ===
Data shape: (150, 4)
Samples: 150, Features: 4
Classes: ['setosa' 'versicolor' 'virginica']
Feature names: ['sepal length (cm)', 'sepal width (cm)', 'petal length (cm)', 'petal width (cm)']

=== Manual PCA (Your Method) ===
Final A shape: (150, 4)
the rank of A is 4
Centered B shape: (150, 4)
the rank of B is 4
SVD Results:
U shape: (150, 4), s length: 4, Vt shape: (4, 4)
Principal component scores shape: (150, 4)

Explained Variance Ratio (Manual PCA):
PC1: 92.462%
PC2: 5.307%
PC3: 1.710%
PC4: 0.521%
                        PC1       PC2       PC3       PC4
sepal length (cm)  0.361387 -0.656589  0.582030  0.315487
sepal width (cm)  -0.084523 -0.730161 -0.597911 -0.319723
petal length (cm)  0.856671  0.173373 -0.076236 -0.479839
petal width (cm)   0.358289  0.075481 -0.545831  0.753657

sklearn PCA

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
# 6. 使用sklearn的PCA进行对比
from sklearn.decomposition import PCA

print("\n=== sklearn PCA Comparison ===")

pca_sklearn = PCA()
scores_sklearn = pca_sklearn.fit_transform(A)

print(f"sklearn PCA scores shape: {scores_sklearn.shape}")
print(f"Explained variance ratio (sklearn):")
for i, ratio in enumerate(pca_sklearn.explained_variance_ratio_[:4]):
    print(f"PC{i+1}: {ratio:.3%}")

# 绘制sklearn PCA结果
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
for i, target in enumerate(target_names):
    idx = y == i
    plt.scatter(scores_sklearn[idx, 0], scores_sklearn[idx, 1], 
               alpha=0.7, label=target, s=50)

plt.title("sklearn PCA: Iris Dataset\n(Reference comparison)")
plt.xlabel(f"PC1 ({pca_sklearn.explained_variance_ratio_[0]:.2%} variance)")
plt.ylabel(f"PC2 ({pca_sklearn.explained_variance_ratio_[1]:.2%} variance)")
plt.legend()
plt.grid(False)
plt.tight_layout()
plt.show()

correlation = np.corrcoef(scores_manual[:, 0], scores_sklearn[:, 0])[0, 1]

print(f"Correlation between Manual PCA PC1 and sklearn PCA PC1: {correlation:.4f}")

sklearn PCA 运行结果

1
2
3
4
5
6
7
8
9
10
=== sklearn PCA Comparison ===
sklearn PCA scores shape: (150, 4)
Explained variance ratio (sklearn):
PC1: 92.462%
PC2: 5.307%
PC3: 1.710%
PC4: 0.521%

Correlation between Manual PCA PC1 and sklearn PCA PC1: 1.0000

3 写在最后

谁能想到在生物学中最经常使用的PCA分析涉及的知识会这么广泛，可能有人会说，PCA已经包装成函数了，没有必要去探究背后的数学原理，真的是这样吗？我对此也没有答案。不过我相信无用为大用，看似现在没用，未来可能有用，即使未来没用，我们也获得了探寻的快乐。PCA的核心知识就是线性代数中的奇异值分解，注意PCA是一种线性降维的方法，还有非线性降维的方法，这就触及到笔者的知识盲区了。不求进步多快，只要每天进步一点点，and happy everyday!

其他学习资料

瑞利商定理：
https://www.cnblogs.com/xinyy-HilbertSpace/articles/16013387.html