一、主成分分析的数学方法
主成分分析（Principal Components Analysis，PCA）是一种基于统计的数据降维方法，又称主元素分析、主分量分析。
主成分分析只需要特征值分解，就可以对数据进行压缩、去噪，应用非常广泛。众多原始变量之间往往具有一定的相关关系。
这意味着相关变量所反映的信息有一定程度的重叠，因此可以用较少的综合指标聚合、反映众多原始变量所包含的全部信息或主要信息。
主成分分析方法研究特征变量之间的相关性、相似性，将一组相关性高的高维变量转换为一组彼此独立、互不相关的低维变量，从而降低数据的维数。主成分分析方法的思想是，将高维特征（p维）映射到低维空间（k维）上，新的低维特征是在原有的高维特征基础上通过线性组合而重构的，并具有相互正交的特性，称为主成分特性。
通过正交变换构造彼此正交的新的特征向量，这些特征向量组成了新的特征空间。将特征向量按特征值排序后，样本数据集中所包含的全部方差，大部分就包含在前几个特征向量中，其后的特征向量所含的方差很小。因此，可以只保留前 k个特征向量，而忽略其它的特征向量，实现对数据特征的降维处理。主成分分析的基本步骤是：对原始数据归一化处理后求协方差矩阵，再对协方差矩阵求特征向量和特征值；对特征向量按特征值大小排序后，依次选取特征向量，直到选择的特征向量的方差占比满足要求为止。
主成分分析方法得到的主成分变量具有几个特点：
（1）每个主成分变量都是原始变量的线性组合；
（2）主成分的数目大大少于原始变量的数目；
（3）主成分保留了原始变量的绝大多数信息；
（4）各主成分变量之间彼此相互独立。算法的基本流程如下：
（1）归一化处理，数据减去平均值；
（2）通过特征值分解，计算协方差矩阵；
（3）计算协方差矩阵的特征值和特征向量；
（4）将特征值从大到小排序；
（5）依次选取特征值最大的 k个特征向量作为主成分，直到其累计方差贡献率达到要求；
（6）将原始数据映射到选取的主成分空间，得到降维后的数据。主成分分析方法的主要优点是：
（1）仅以方差衡量信息量，不受数据集以外的因素影响；
（2）各主成分之间正交，可消除原始数据各变量之间的相互影响；
（3）方法简单，易于实现。在图像处理中，把每幅二维图像拉伸为一维向量，即展平为一维数组。
一组 m 幅图像就构造为一个 m 维向量，使用 Karhunen-Loève transform（KLT） 变换得到变换矩阵，选取特征值最大的 k个特征向量作为主成分，从而实现特征降维。图像压缩过程是把一组原始图像变换成低维向量的过程，图像重建就是由低维向量变换重建图像组的过程。
使用主成分分析进行图像压缩和重建会有少量信息损失，但可以把损失控制到很小。-
对于一组 P 维向量 X，通过线性变换转换为另一组 K 维向量 Y。
向量 Xi 包含 m 个数据样本，记为：
计算平均值
计算协方差矩阵：
协方差矩阵C_x是实对称矩阵，可以协方差矩阵的特征向量和对应的特征值。
通过变换矩阵 A 将 X 映射为由 Y 表示的向量：
Y=A∗(X−u)
Y 的协方差矩阵用对角阵 A 和 Cx 表示，即：
则 Cy 是一个对角阵：
是Cx 的特征值，将其按降序排列。
选择前 K 个特征向量，作为主成分特征，可以重建 p 维向量 X，并使均方误差最小。
二、例程
- 14.15：特征描述之主分量
本例程的图像来自 R.C.Gonzalez 《数字图像处理（第四版）》P622 例11.16。
原始图像显示了对应于 6 个波段的 6 幅多光谱卫星图像：可见蓝光（450-520nm）、可见绿光（520-600nm）、可见红光（630-690nm）、近红外（760-900nm）、中红外（1550-1750nm）和热红外（10400-12500nm）。
本例的目的是说明如何使用主分量作为图像特征。
# # 14.15 特征描述之主成分分析
# 读取光谱图像组
img = cv2.imread("../images/Fig1138a.tif", flags=0)
height, width = img.shape[:2] # (564, 564)
nBands = 6 # 光谱波段种类
snBands = ['a','b','c','d','e','f'] # Fig1138a~f
imgMulti = np.zeros((height, width, nBands)) # (564, 564, 6)
mbMatrix = np.zeros((img.size, nBands)) # (318096, 6)
print(imgMulti.shape, mbMatrix.shape)
# 显示光谱图像组
fig1 = plt.figure(figsize=(9, 6)) # 原始图像，6 个不同波段
fig1.suptitle("Spectral image of multi bands by NASA")
for i in range(nBands):
path = "../images/Fig1138{}.tif".format(snBands[i])
imgMulti[:,:,i] = cv2.imread(path, flags=0) # 灰度图像
ax1 = fig1.add_subplot(2,3,i+1)
ax1.set_xticks([]), ax1.set_yticks([])
ax1.imshow(imgMulti[:,:,i], 'gray') # 绘制光谱图像 snBands[i]
plt.tight_layout()
# 主成分分析 (principal component analysis)
mbMean = np.zeros((nBands,))
for i in range(nBands):
mbArray = imgMulti[:,:,i].flatten() # 转为一维数组
mbMean[i] = mbArray.mean()
mbMatrix[:,i] = mbArray - mbMean[i] # 数据标准化 (318096, 6)
# cov = np.cov(mbMatrix.transpose(),rowvar=0) # 以列为变量计算协方差矩阵 (6,6)
cov = np.cov(mbMatrix, rowvar=0) # 以列为变量计算协方差矩阵 (6,6)
# eigenvalues: 特征值，一维数组 (M,)，eigenvectors: 特征向量，二维矩阵 (M,M)
eigenValues, eigenVectors = np.linalg.eig(cov) # 计算特征值和特征向量
principal = np.matmul(mbMatrix, eigenVectors) # 主元素变换，投影到特征向量方向 (318096, 6)
print(eigenValues.shape, eigenVectors.shape, principal.shape)
print("Eigenvalues:\n", eigenValues.round(4)) # 特征值，从大到小
# print("Eigenvectors:\n", eigenVectors.round(4))
# 显示主成分变换图像
fig2 = plt.figure(figsize=(9, 6)) # 主元素图像，q=6
fig2.suptitle("Principal component images")
imgPCA = np.zeros((height, width, nBands)) # (564, 564, 6)
for i in range(nBands):
pca = principal[:, i].reshape(-1, img.shape[1]) # 主元素图像 (564, 564)
imgPCA[:,:,i] = np.uint8(cv2.normalize(-pca, (height, width), 0, 255, cv2.NORM_MINMAX))
ax2 = fig2.add_subplot(2,3,i+1)
ax2.set_xticks([]), ax2.set_yticks([])
ax2.imshow(imgPCA[:,:,i], 'gray') # 绘制主成分图像
plt.tight_layout()
# 保留的主成分数量
remainRatio = 0.95 # 设定的主成分累计方差贡献率
remainValue = remainRatio * sum(eigenValues)
IndexDes = np.argsort(-eigenValues) # 特征值降序排列的索引
topKValue = 0.0
for i in range(len(eigenValues)):
topKValue += eigenValues[IndexDes[i]] # 降序排列的第 i 个特征值
print("k={}, topKValue = {:.2f}, topKRatio = {:.4f}"
.format(i, topKValue, topKValue/sum(eigenValues)))
if topKValue > remainValue:
K = i + 1 # （0，1，2）-> K=3
break
print("number of PCA features: K=", K) # 主成分方差贡献率 95% 时的特征维数 K=3
indexPCA = IndexDes[:K] # 选择特征值最大的 K 个特征向量的索引
eigenVectPCA = eigenVectors[:, indexPCA] # 选择 K 个主要特征向量组成降维特征矩阵 (P=6, K=3)
print("PCA eigenvalues:\n", eigenValues[indexPCA].round(4)) # 特征值，从大到小
print("PCA eigenvectors:\n", eigenVectPCA.round(4))
三、资料
youcans_的博客：
https://blog.csdn.net/youcans/article/details/125761655
