声明：本文基于在校课程及吴恩达ML教程，代码参考自多份博客（已在参考链接中表明），如需转载请标明出处。
源代码、实验数据、实验指导书： https://pan.baidu.com/s/1yTTI0_w2bZ7o8uuxravLyA 提取码: spbp
PCA 是最常见的降维算法。
比如说，现在想把
2
2
2 维 降到
1
1
1 维。PCA 的做法是找到一个方向向量（Vector direction），当把所有的数据都投射到该向量上时，投射均方误差能尽可能地小。方向向量是一个经过原点的向量，而投射误差是从特征向量向该方向向量作垂线的长度。
现在将问题一般化，即将 n 维数据降至 k 维，对于 PCA 来说就是找到高维向量
U
=
(
u
1
,
u
2
,
.
.
.
,
u
k
)
U = (u^1, u^2, ... , u^k)
U=(u1,u2,...,uk) 使得总的投射误差最小。
z-score 标准化(zero-mean normalization)，也叫标准差标准化。这种方法给予原始数据的均值（mean）和标准差（standard deviation）进行数据的标准化。经过处理的数据符合标准正态分布，即均值为 0 0 0，标准差为 1 1 1。注意，一般来说 z-score 不是归一化，而是标准化，归一化只是标准化的一种。from 数据标准化/归一化normalization
f o r e a c h c o l u m n , ( X − m e a n ) / s t d for \ each \ column, (X - mean) / std for each column,(X−mean)/std
其中 mean 为所有样本数据的均值，std 为所有样本数据的标准差。
x ( i ) x^{(i)} x(i) 代表第 i i i 个样本，是一个 n ∗ 1 n *1 n∗1 的向量， n n n 是特征维数， m m m 是样本个数。
U U U 是 n ∗ n n*n n∗n 维的主成分（Principal Component）集
- U r e d u c e U_{reduce} Ureduce 是 n ∗ k n*k n∗k 维的前 k k k 个主成分
- U r e d u c e T U_{reduce}^T UreduceT 是前 k k k 个主成分的转置（ k ∗ n k*n k∗n 维）， x ( i ) x^{(i)} x(i) 是 n ∗ 1 n*1 n∗1 维的均值归一化后第 i i i 个样本数据，因此结果 z ( i ) z^{(i)} z(i) 为 k ∗ 1 k*1 k∗1 维的投影向量
- 整体使用矩阵计算即 m ∗ n m*n m∗n 维的 X X X 【矩阵乘】 n ∗ k n*k n∗k 维的 U r e d u c e U_{reduce} Ureduce 即可。 Z = X U r e d u c e Z = XU_{reduce} Z=XUreduce
若想使用压缩后的数据近似的获得原有的特征，可以使用如下公式：
因为：
z
(
i
)
=
U
r
e
d
u
c
e
T
x
(
i
)
z^{(i)} = U_{reduce}^Tx^{(i)}
z(i)=UreduceTx(i)
所以相反的方程：
x
a
p
p
r
o
x
(
i
)
=
U
r
e
d
u
c
e
z
(
i
)
x
a
p
p
r
o
x
(
i
)
≈
x
(
i
)
X
a
p
p
r
o
x
=
Z
U
r
e
d
u
c
e
T
x^{(i)}_{approx} = U_{reduce}z^{(i)} \\ x^{(i)}_{approx} \approx x^{(i)} \\ X_{approx} = Z U_{reduce}^T
xapprox(i)=Ureducez(i)xapprox(i)≈x(i)Xapprox=ZUreduceT
- U r e d u c e U_{reduce} Ureduce n ∗ k n*k n∗k 维， z ( i ) z^{(i)} z(i) 为 k ∗ 1 k*1 k∗1 维，可得 n ∗ 1 n*1 n∗1 维的 x a p p r o x ( i ) x^{(i)}_{approx} xapprox(i)
- 整体使用矩阵计算， Z Z Z 是 m ∗ k m*k m∗k 维， U r e d u c e T U_{reduce}^T UreduceT 是 k ∗ n k*n k∗n 维，可得 m ∗ n m*n m∗n 维的 X a p p r o x X_{approx} Xapprox
x
a
p
p
r
o
x
(
i
)
=
U
r
e
d
u
c
e
z
(
i
)
x
a
p
p
r
o
x
(
i
)
≈
x
(
i
)
x^{(i)}_{approx} = U_{reduce}z^{(i)} \\ x^{(i)}_{approx} \approx x^{(i)}
xapprox(i)=Ureducez(i)xapprox(i)≈x(i)
如图所知，这是一个漂亮的重构，它们与原始数据相当相似。上图直观的展现了从低维表示
Z
Z
Z 回到未压缩的表示。
import numpy as np
class PCA():
def normalize(self, X):
"""
1. 均值归一化
1. 计算出每一维特征的均值𝜇_𝑗 ，令 𝑥_𝑗=𝑥_𝑗−𝜇_𝑗。
2. 如果特征是在不同的数量级上，还需要将其除以标准差 。
for each column, (X - mean) / std
"""
means = X.mean(axis=0)
stds = X.std(axis=0, ddof=1) # numpy.std() 求标准差的时候默认是除以 n 的，即是有偏的，np.std无偏样本标准差方式为加入参数 ddof = 1；
X_norm = (X - means) / stds
return X_norm
def covariance_matrix(self, X_norm):
"""
2. 计算协方差矩阵（covariance matrix）𝛴
∑=1/𝑚 ∑1_(𝑖=1)^m▒〖(𝑥^((𝑖)) ) (𝑥^((𝑖)) )^𝑇 〗= 1/𝑚 𝑋𝑋^𝑇
∑ = 1/𝑚〖𝑋^𝑇 𝑋〗
"""
m = X_norm.shape[0]
return (X_norm.T @ X_norm) / m
def dimensional_reduction(self, X, keep_dims=None):
if not keep_dims:
keep_dims = X.shape[1] - 1
# 1. 均值、归一化
normalize_x = self.normalize(X)
# 2. 计算协方差矩阵
cov_x = self.covariance_matrix(normalize_x)
# 3. 计算协方差矩阵𝜮的特征向量, 使用奇异值分解（SVD分解）
U, S, V = np.linalg.svd(cov_x) # U: principle components (n, n)
# 4. 将数据从𝑛维降至𝑘维，从𝑈中选取前𝑘个向量
U_reduce = U[:, :keep_dims] # U_reduce.shape : (n, k)
# 5. 将二维数据投影到主成分方向（二维数据降维到一维）得到降维的结果
Z = normalize_x @ U_reduce # Z.shape : (m, k)
return normalize_x, U_reduce, Z
def recover(self, Z, U_reduce):
"""
6. 数据恢复：将降维后一维数据再投影回二维空间。
𝒙_𝒂𝒑𝒑𝒓𝒐𝒙^((𝒊) )=𝑼_𝒓𝒆𝒅𝒖𝒄𝒆 𝒛^((𝒊))
"""
return Z @ U_reduce.T
import scipy.io as scio
import numpy as np
import matplotlib.pyplot as plt
from PCA import PCA
# 读取图片数据
data = scio.loadmat('ex7data1.mat')
X = data['X'] # X.shape : (50, 2)
if __name__ == '__main__':
plt.figure(0, (7, 7))
plt.title("raw data")
plt.xlim(xmin=0, xmax=7)
plt.ylim(ymin=2, ymax=8)
plt.scatter(X[:, 0], X[:, 1], c='', marker='o', facecolors='none', edgecolors='b')
pca = PCA()
X_norm, U_reduce, Z = pca.dimensional_reduction(X)
print("U_reduce : ", U_reduce)
print("Z[0] : ", Z[0])
X_approx = pca.recover(Z, U_reduce)
print("X_approx[0] : ", X_approx[0])
plt.figure(1, (7, 7))
plt.scatter(X_norm[:, 0], X_norm[:, 1], c='', marker='o', facecolors='none', edgecolors='b')
plt.scatter(X_approx[:, 0], X_approx[:, 1], c='', marker='o', facecolors='none', edgecolors='r')
for i in range(X.shape[0]):
x = [X_norm[i, 0], X_approx[i, 0]]
y = [X_norm[i, 1], X_approx[i, 1]]
plt.plot(x, y, '--k', linewidth=1)
plt.title("PCA dimensional reduction and recover")
plt.xlim(xmin=-4, xmax=3)
plt.ylim(ymin=-4, ymax=3)
plt.show()
计算出的第一主成分：U_reduce : [[-0.70710678] [-0.70710678]]
第一个样本的投影值：Z[0] : [1.48127391]
第一个样本投影回二维数据空间位置：X_approx[0] : [-1.04741883 -1.04741883]
import scipy.io as scio
import numpy as np
import matplotlib.pyplot as plt
from PCA import PCA
from displayData import display_data
# 读取图片数据
data = scio.loadmat('ex7faces.mat')
X = data['X'] # X.shape : (5000, 1024)
if __name__ == '__main__':
# 前100个人脸
display_data(X, "Original faces")
pca = PCA()
X_norm, U_reduce, Z = pca.dimensional_reduction(X, 100)
# 前100维主成分（转换成人脸形式）
display_data(U_reduce.T, "Principle Components", 6, 6)
X_approx = pca.recover(Z, U_reduce)
# 100个人脸仅使用前100个主成分表示再恢复得到的人脸
display_data(X_approx, "Recovered faces")
展示人脸数据的函数（displayData.py）
# displayData.py
import matplotlib.pyplot as plt
import numpy as np
# 函数说明：把输入的图像数据X进行重新排列，显示在一个面板figurePane中，
# 面板中有多个小imge用来显示每一行数据
def display_data(x, title = "show top 100", rows = 10, cols = 10):
(m,n) = x.shape
# 设置每个小图例的宽度和高度
width = np.round(np.sqrt(n)).astype(int)
height = (n / width).astype(int)
print("width : ", width, " height : ", height)
# 设置图片的行数和列数
# rows = 10 # np.floor(np.sqrt(m)).astype(int)
# cols = 10 # np.ceil(m / rows).astype(int)
print("rows : ", rows, " cols : ", cols)
# 设置图例之间的间隔
pad = 1
# 初始化图像数据
display_array = -np.ones((pad + rows*(height+pad),
pad + cols*(width + pad)))
print(display_array.shape)
# 把数据按行和列复制进图像中
# current_image = np.random.randint(0, m)
current_image = 0
for j in range(rows):
for i in range(cols):
if current_image >= m:
break
# [:,np.newaxis]可以让指定的那一列数据以列的形式返回和指定
# 否则返回的是行的形式
max_val = np.max(np.abs(x[current_image,:]))
display_array[pad + j*(height + pad) + np.arange(height),
pad + i*(width + pad) + np.arange(width)[:,np.newaxis]] = \
x[current_image,:].reshape((height,width)) / max_val
# current_image = np.random.randint(0, m)
current_image += 1
if current_image >= m :
break
# 显示图像
plt.figure(figsize=(8, 8))
# 设置图像色彩为灰度值，指定图像坐标范围
plt.imshow(display_array,cmap = 'gray',extent =[-1,1,-1,1])
plt.axis('off')
plt.title(title, fontsize=20)
plt.show()
