python pca/svd原理及应用

PCA的定义

• 主成分分析，即Principal Component Analysis（PCA），是多元统计中的重要内容，也广泛应用于机器学习和其它领域。它的主要作用是对高维数据进行降维。PCA把原先的n个特征用数目更少的k个特征取代，新特征是旧特征的线性组合，这些线性组合最大化样本方差，尽量使新的k个特征互不相关。

PCA的实现步骤

1. 计算样本每个特征的平均值；

2. 每个样本数据减去该特征的平均值；

3. 对样本组成的矩阵求协方差矩阵；

4. 对协方差矩阵做分解得到特征值和特征向量，并将特征值和特征向量重新按照降序排列；

5. 对特征值求取累计贡献率；

6. 通过累计贡献率选取适合的特征向量子集合；

7. 对原始数据进行转换

---

注意：假设有m个样本，每个样本有n个特征，即组成m行n列的矩阵A,如果是图像的话，由于图像是二维向量，所以先要将图像转换为一个行向量。求均值是对A矩阵的每一列求均值，此处协方差矩阵应该是A.T * A而不是A*A.T,最后得到的协方差矩阵的形状应为n*n。

协方差矩阵分解

1. 原始算法（硬算，特征数很大时计算量非常大）

2. numpy的svd算法来分解（使用svd分解出来的右奇异矩阵和直接求出来的特征向量矩阵是一个意思，但是svd并不是通过协方差矩阵来计算的，而是通过某种强大的算法近似求解）

3. scikit-learn实现的PCA

PCA实现

mat = [[-1,-1,0,2,1],[2,0,0,-1,-1],[2,0,1,1,0]]
mat = np.array(mat)
#原始算法
def origin_algorithm_pca(mat, 2):
    #对列求均值
    mean = np.mean(mat, 0)
    #去中心化
    A = mat - mean
    #求方差矩阵
    cov = np.dot(A.T, A)/2
    #求解协方差矩阵的特征值和特征向量
    s, vector = np.linalg.eig(cov)
    #对特征进行压缩
    compression = np.dot(A, vector[:2,:].T)
    return compression

def numpy_svd(mat, 2):
    #对列求均值
    mean = np.mean(mat, 0)
    #去中心化
    A = mat - mean
    #
    u,s,VT = np.linalg.svd(A)
    #压缩
    compression = np.dot(A, VT[:2,:])
    return compression

PCA实现人脸识别

1. 将每张图像拉直为一维行向量

2. 将所有样本存储为矩阵形式，矩阵的每一列表示图像的一个特征也是一个像素

3. 按列求取均值，并用原始矩阵减去均值

4. 利用A.T * A得到协方差矩阵(这一步代码中直接利用svd方法得到U,S,VT)

5. 将原始矩阵乘以VT得到压缩后的矩阵，将VT存储起来用于测试测试样本

6. train_A*VT和test_A*VT，将test_A*VT中的每个样本与train_A*VT计算平方差，并按行求和后排序，得到距离最小的对应的train_A的label，就是预测结果

7. 根据预测结果和真实标签计算准确率

数据集地址：

链接：https://pan.baidu.com/s/1T1DQxHBTLa37rJ8TIbtOlQ

提取码：0qvx

# coding:utf-8
import os
from numpy import *
import numpy as np
import cv2
import matplotlib.pyplot as plt
from pylab import mpl

mpl.rcParams['font.sans-serif'] = ['SimHei']


# 图片矢量化
def img2vector(image):
    img = cv2.imread(image, 0)  # 读取图片
    rows, cols = img.shape  #获取图片的像素
    imgVector = np.zeros((1, rows * cols))
    imgVector = np.reshape(img, (1, rows * cols))#使用imgVector变量作为一个向量存储图片矢量化信息，初始值均设置为0
    return imgVector


orlpath = "./ORL"


# 读入人脸库,每个人随机选择k张作为训练集,其余构成测试集
def load_orl(k):#参数K代表选择K张图片作为训练图片使用
    '''
    对训练数据集进行数组初始化，用0填充，每张图片尺寸都定为112*92,
    现在共有40个人，每个人都选择k张，则整个训练集大小为40*k,112*92
    '''
    train_face = np.zeros((40 * k, 112 * 92))
    train_label = np.zeros(40 * k)  # [0,0,.....0](共40*k个0)
    test_face = np.zeros((40 * (10 - k), 112 * 92))
    test_label = np.zeros(40 * (10 - k))
    # sample=random.sample(range(10),k)#每个人都有的10张照片中，随机选取k张作为训练样本(10个里面随机选取K个成为一个列表)
    sample = random.permutation(10) + 1  # 随机排序1-10 (0-9）+1
    for i in range(40):  # 共有40个人
        people_num = i + 1
        for j in range(10):  # 每个人都有10张照片
            image = orlpath + '/s' + str(people_num) + '/' + str(sample[j]) + '.jpg'
            # 读取图片并进行矢量化
            img = img2vector(image)
            if j < k:
                # 构成训练集
                train_face[i * k + j, :] = img
                train_label[i * k + j] = people_num
            else:
                # 构成测试集
                test_face[i * (10 - k) + (j - k), :] = img
                test_label[i * (10 - k) + (j - k)] = people_num

    return train_face, train_label, test_face, test_label


# 定义PCA算法
def PCA(data, r):#降低到r维
    data = np.float32(np.mat(data))
    rows, cols = np.shape(data)
    # print(rows, cols)
    data_mean = np.mean(data, 0)  # 对列求平均值
    A = data - np.tile(data_mean, (rows, 1))  # 将所有样例减去对应均值得到A
    u, s, VT = np.linalg.svd(A) #利用svd求解右奇异向量即需要将原始数组映射的向量空间矩阵
    V_r = VT[:, 0:r]  # 按列取前r个特征向量
    #将原始数据乘上新的空间向量得到降维后的矩阵
    final_data = A * V_r
    return final_data, data_mean, V_r

def face_recongize():
    for r in range(10, 81, 10):  # 最多降到40维,即选取前40个主成分（因为当k=1时，只有40维)
        print("当降维到%d时" % (r))
        x_value = []
        y_value = []
        train_face, train_label, test_face, test_label = load_orl(7)  # 得到数据集

        # 利用PCA算法进行训练
        data_train_new, data_mean, V_r = PCA(train_face, r)
        # print(data_train_new.shape)
        num_train = data_train_new.shape[0]  # 训练脸总数
        num_test = test_face.shape[0]  # 测试脸总数
        temp_face = test_face - np.tile(data_mean, (num_test, 1))
        data_test_new = temp_face * V_r  # 得到测试脸在特征向量下的数据
        data_test_new = np.array(data_test_new)  # mat change to array
        # print(data_test_new.shape)
        data_train_new = np.array(data_train_new)

        # 测试准确度
        true_num = 0
        for i in range(num_test):
            testFace = data_test_new[i, :]
            # print(testFace.shape)
            diffMat = data_train_new - np.tile(testFace, (num_train, 1))  # 训练数据与测试脸之间距离
            print(diffMat.shape)
            sqDiffMat = diffMat ** 2
            sqDistances = sqDiffMat.sum(axis=1)  # 按行求和
            sortedDistIndicies = sqDistances.argsort()  # 对向量从小到大排序，使用的是索引值,得到一个向量
            indexMin = sortedDistIndicies[0]  # 距离最近的索引
            if train_label[indexMin] == test_label[i]:
                true_num += 1
            else:
                pass

            accuracy = float(true_num) / num_test
            x_value.append(7)
            y_value.append(round(accuracy, 2))

        print('当每个人选择%d张照片进行训练时，The classify accuracy is: %.2f%%' % (7, accuracy * 100))

if __name__ == '__main__':
    face_recongize()

PCA的优缺点

1. 优点

• 仅仅需要以方差衡量信息量，不受数据集以外的因素影响。

• 各主成分之间正交，可消除原始数据成分间的相互影响的因素

• 计算方法简单，主要运算是特征值分解，易于实现。

2. 缺点

• 主成分各个特征维度的含义具有一定的模糊性，不如原始样本特征的解释性强。

• 方差小的非主成分也可能含有对样本差异的重要信息，因降维丢弃可能对后续数据处理有影响。

SVD原理

任意形状的矩阵A经过svd分解会得到U，Σ，VT（V的转置）三个矩阵，其中U是一个M×M的方阵，被称为左奇异向量，方阵里面的向量是正交的；Σ是一个M×N的对角矩阵，除了对角线的元素其他都是0，对角线上的值称为奇异值；VT(V的转置)是一个N×N的矩阵，被称为右奇异向量，方阵里面的向量也都是正交的。

可以计算ΣΣ前x个奇异值的平方和占所有奇异值的平方和的比例，如果大于90%，我们就选这x个奇异值重构矩阵（剩余的数据代表的可能是噪声，无用数据）。

那么重构的A'为：

SVD的应用

• svd实现图像压缩

import scipy.misc
from scipy.linalg import svd
import matplotlib.pyplot as plt
import numpy as np
import numpy
img = scipy.misc.face()[:,:,0]
img = np.array(img)
U,s,Vh=svd(img)
plt.gray()
#将特征压缩到十维，并还原图像
A = numpy.dot(U[:,0:10],numpy.dot(numpy.diag(s[0:10]),Vh[0:10,:]))
plt.subplot(121,aspect='equal')
plt.title("img")
plt.imshow(img)

plt.subplot(111,aspect='equal')
plt.title(":10")
plt.imshow(A)

• svd实现信号去噪

import matplotlib.pylab as plt
import numpy as np
import numpy
from scipy.linalg import svd
x = np.linspace(-np.pi, np.pi, 400)
mu, sigma = 0, 0.05
ss = np.random.normal(mu, sigma, 400)
a = (np.sin(x) + ss).reshape(20, 20)

U,s,Vh = svd(a)
print s
A = numpy.dot(U[:,0:2],numpy.dot(numpy.diag(s[0:2]),Vh[0:2,:]))
aa = A.reshape((400,))

plt.subplot(221,aspect='equal')
plt.plot(x, np.sin(x))
plt.title("sin")

plt.subplot(222,aspect='equal')
plt.plot(x, np.sin(x) + ss)
plt.title("sin with noise")

plt.subplot(223,aspect='equal')
plt.title("sin remove noise x = 2")
plt.plot(x, aa)

A = numpy.dot(U[:,0:1],numpy.dot(numpy.diag(s[0:1]),Vh[0:1,:]))
aa = A.reshape((400,))
plt.subplot(224,aspect='equal')
plt.title("sin remove noise x = 1")
plt.plot(x, aa)

plt.show()

第一副图是纯sin波形，第2副图是加了噪声的波形，第三副选择2个奇异值重构A结果很接近标准sin波形实现了去噪，第四副图是选了1个奇异值重构A的结果。

参考文献

http://liao.cpython.org/scipy06.html

https://www.cnblogs.com/jclian91/p/8024101.html

https://blog.csdn.net/m0_50572604/article/details/121018988

https://www.cnblogs.com/pinard/p/6239403.html

https://www.cnblogs.com/pinard/p/6251584.html

https://github.com/Gaoshiguo/PCA-Principal-Components-Analysis