PCA - 主成分分析（Python实现图像压缩）

PCA(Principal Component Analysis)

PCA的本质就是找一些投影方向，使得数据在这些投影方向上的方差最大，而且这些投影方向是相互正交的。
这其实就是找新的正交基的过程，计算原始数据在这些正交基上投影的方差，方差越大，就说明在对应正交基上包含了更多的信息量。
而原始数据协方差矩阵的特征值越大，对应的方差越大，在对应的特征向量上投影的信息量就越大。反之，如果特征值较小，
则说明数据在这些特征向量上投影的信息量很小，可以将小特征值对应方向的数据删除，从而达到了降维的目的。

算法思想

将一组P维向量降为K维（K<P），其目标是选择K个单位正交基，使得原始数据变换到这组基上后，
各字段两两间协方差为0，而字段的方差则尽可能大——这不就是对角矩阵吗？

Y = XP，其中X为原始数据（每个维度按列组成），Y为X对P做基变换后的数据，P的基行列组成。
设X、Y的协方差矩阵（对称矩阵）分别为C、D，则：D=PCP’。

我们希望D为对角矩阵，则D为C对P做相似对角化所得。
由于协方差矩阵为对称矩阵，则构成相似变换矩阵P的特征向量是正交的，实现逆变换。

实验步骤

（1）分别求每个维度的平均值，然后对于所有的样例，都减去对应维度的均值，得到去中心化的数据；

（2）求协方差矩阵C：用去中心化的数据矩阵乘上它的转置，然后除以(N-1)即可，N为样本数量；

（3）求协方差的特征值和特征向量；

（4）将特征值按照从大到小排序，选择前k个，然后将其对应的k个特征向量分别作为列向量组成特征向量矩阵；

（5）将样本点从原来维度投影到选取的k个特征向量，得到低维数据；

（6）通过逆变换，重构低维数据，进行复原。

import numpy as np

创建 eigValPct(eigVals, percentage)

通过方差的百分比来计算将数据降到多少维。

函数传入的参数是特征值eigVals和百分比percentage，返回需要降到的维度数num

def eigValPct(eigVals, percentage):
    sortArray=np.sort(eigVals)[::-1] # 特征值从大到小排序
    pct = np.sum(sortArray)*percentage
    tmp = 0
    num = 0
    for eigVal in sortArray:
        tmp += eigVal
        num += 1
        if tmp>=pct:
            return num

创建 im_PCA(dataMat, percentage=0.9)

函数有两个参数，其中dataMat是已经转换成矩阵matrix形式的数据集，每列表示一个维度；

其中的percentage表示取前多少个特征需要达到的方差占比，默认为0.9。

注意：

np.cov(dataMat, rowvar=False)，按照rowvar的默认值，会把一行当成一个特征，一列当成一个样本。

def im_PCA(dataMat, percentage=0.9):
    meanVals = np.mean(dataMat, axis=0)
    meanRemoved = dataMat - meanVals
    # 这里不管是对去中心化数据 or 原始数据计算协方差矩阵，结果都一样，特征值大小会变，但相对大小不会改变
    # covMat = np.cov(dataMat, rowvar=False)
    # 标准的计算需要除以(dataMat.shape[0]-1)，不算也不会影响结果，理由同上
    covMat = np.dot(np.transpose(meanRemoved), meanRemoved)
    
    eigVals, eigVects = np.linalg.eig(np.mat(covMat))
    k = eigValPct(eigVals,percentage) #要达到方差的百分比percentage，需要前k个向量
    print('K =', k)
    
    eigValInd = np.argsort(eigVals)[::-1]  #对特征值eigVals从大到小排序
    eigValInd = eigValInd[:k]
    redEigVects = eigVects[:,eigValInd] #主成分
    lowDDataMat = meanRemoved*redEigVects #将原始数据投影到主成分上得到新的低维数据lowDDataMat
    reconMat = (lowDDataMat*redEigVects.T)+meanVals   #得到重构数据reconMat
    return lowDDataMat, reconMat

注意：

图像Matrix格式必须转换为uint8格式，否则cv2.imshow时图像不能正常显示！np.array(reconMat, dtype=‘uint8’)。

但是，强制类型转化后会丢失信息，比如将 6. 变成 5 ，因为强制类型转化是直接采用截断二进制位的方式。。。

import cv2

img = cv2.imread('ava.jpg')
blue = img[:,:,0]
dataMat = np.mat(blue)
lowDDataMat, reconMat = im_PCA(dataMat, 1)
print('原始数据', blue.shape, '降维数据', lowDDataMat.shape)
# print(dataMat)
# print(reconMat)
# 格式必须转换为uint8格式，这里丢失了很多信息！！！
# reconMat = np.array(reconMat, dtype='uint8')

cv2.imshow('blue', blue)
cv2.imshow('reconMat', np.array(reconMat, dtype='uint8'))
cv2.waitKey(0)

K = 426
原始数据 (640, 426) 降维数据 (640, 426)

计算降维误差率

根据PCA的误差计算方式：信息丢失率 =

注意：np.dot(A, B)

（1）如果处理的是一维数组，则得到的是两数组的内积；

（2）如果是二维数组/矩阵之间的运算，得到的是矩阵积；

def PrintError(data, recdata):
    sum1 = 0
    sum2 = 0
    D_value = data - recdata # 计算两幅图像之间的差值矩阵
    # 计算两幅图像之间的误差率，即信息丢失率
    for i in range(data.shape[0]):
        sum1 += np.dot(data[i],data[i])
        sum2 += np.dot(D_value[i], D_value[i])
    print('丢失信息量：', sum2)
    print('原始信息量：', sum1)
    print('信息丢失率：', sum2/sum1)

print(dataMat)
print(reconMat)
PrintError(np.array(blue, dtype='double'), np.array(reconMat, dtype='double'))

[[ 6  6  6 ... 11 13 13]
 [ 6  6  6 ... 12 15 15]
 [ 7  7  7 ... 13 17 17]
 ...
 [16 16 16 ... 29 29 29]
 [16 16 16 ... 29 29 29]
 [16 16 16 ... 29 29 29]]
[[ 6.  6.  6. ... 11. 13. 13.]
 [ 6.  6.  6. ... 12. 15. 15.]
 [ 7.  7.  7. ... 13. 17. 17.]
 ...
 [16. 16. 16. ... 29. 29. 29.]
 [16. 16. 16. ... 29. 29. 29.]
 [16. 16. 16. ... 29. 29. 29.]]
丢失信息量： 7.140138346734739e-16
原始信息量： 784219850.0
信息丢失率： 9.104766152928595e-25

使用sklearn的PCA方法

参数说明：

sklearn.decomposition.PCA(n_components=None, copy=True, whiten=False);

n_components: PCA算法中所要保留的主成分个数n，也即保留下来的特征个数n。赋值为int，比如n_components=1，将把原始数据降到一个维度。赋值为string，比如n_components=‘mle’，将自动选取特征个数n，使得满足所要求的方差百分比。

copy: 是否在原始数据的副本上进行运算。

whiten: 白化，使得每个特征具有相同的方差。

from sklearn.decomposition import PCA

pca = PCA(n_components=426).fit(blue)
# 降维
x_new = pca.transform(blue)
# 还原降维后的数据到原空间
recdata = pca.inverse_transform(x_new)
print(recdata)
# 计算误差
PrintError(np.array(blue, dtype='double'), np.array(reconMat, dtype='double'))

cv2.imshow('sklearn-recdata', np.array(recdata, dtype='uint8'))
cv2.waitKey(0)

[[ 6.  6.  6. ... 11. 13. 13.]
 [ 6.  6.  6. ... 12. 15. 15.]
 [ 7.  7.  7. ... 13. 17. 17.]
 ...
 [16. 16. 16. ... 29. 29. 29.]
 [16. 16. 16. ... 29. 29. 29.]
 [16. 16. 16. ... 29. 29. 29.]]
丢失信息量： 7.140138346734739e-16
原始信息量： 784219850.0
信息丢失率： 9.104766152928595e-25

Reference

计算过程

https://www.cnblogs.com/lzllovesyl/p/5235137.html

PCA详解，原理及步骤，包括数学基础

https://my.oschina.net/gujianhan/blog/225241#OSC_h2_1

PCA，图像压缩还原

https://blog.csdn.net/Vincent_zbt/article/details/88648739

scikit-learn中PCA的使用方法

https://blog.csdn.net/u012162613/article/details/42192293/