HIT机器学习lab4

最新推荐文章于 2022-03-26 17:57:03 发布

ln栓

最新推荐文章于 2022-03-26 17:57:03 发布

阅读量266

点赞数

文章标签：机器学习 python opencv

本文链接：https://blog.csdn.net/weixin_44118063/article/details/110927938

版权

1 实验目的
实现一个PCA模型，能够对给定数据进行降维（即找到其中的主成分）
2 实验要求及实验环境
实验要求
首先人工生成一些数据（如三维数据），让它们主要分布在低维空间中，如首先让某个维度的方差远小于其它维度，然后对这些数据旋转。生成这些数据后，用你的PCA方法进行主成分提取。
找一个人脸数据（小点样本量），用你实现PCA方法对该数据降维，找出一些主成分，然后用这些主成分对每一副人脸图像进行重建，比较一些它们与原图像有多大差别（用信噪比衡量）。
实验环境
windows64，VScode，python3.8.6

3算法原理
PCA(主成分分析，Principal Component Analysis)是最常用的一种统计分析、简化数据集的方法。通过PCA可以使要分析的数据的维度降低，且这些维度还会包含原数据集的主要信息。关于PCA的推导有两种方式：最小投影距离和最大投影方差。
最小投影距离：样本点到这个超平面的距离都足够近
最大投影方差：样本点在这个超平面上的投影尽可能分开

4 具体实现
4.1中心化和坐标变化
假设m个n维数据X=(x(1) , x(2),… . , x(m)，我们令x(i)= x(i) - u，其中u=1/m∑x(i)。
这样我们会得到显∑x(i)=0，这就是中心化。中心化的目的是方便求协方差矩阵:中心化后XTX就是样本集的协方差矩阵。我们知道X(i)是一个n维向量，这个n维空间中的基底是(1,0,…,0)”，(0,1,… ,0)T，…，(0,0,…,1)T。如果我们换一组基底，x(i)又是另一种表现形式。
现在假设我们n维空间中的基底时W=(w(1) , w(2) ,-- , w(n)),那么Z= WTX。
假设我们只取W中的前d个向量，即W=(w(1) , w(2) ,-- , w(d))，其中d <n。那么Z=WTX就是d维空间下的一个坐标，即Z是X从n维空间到d维空间的一个投影。
那么我们从d维空间下的坐标Z变化到n维空间，就是X=ZW。总结就是:
某种投影W =(w(1) , w/2) , -… , w(d))投影坐标:Z=wTx。还原坐标： X=ZW显然我们想要
X(i) 这些点的投影越分散越好、距离平面越近越好，就像这样：
在这里插入图片描述

因此我们有两种策略：最小投影距离和最大投影方差。事实上二者相当于优化同一个函数。

4.2 PCA的推导
对于任意一个样本x(i)，在新的坐标系中的投影为WTx(i),在新坐标系中的投影方差WTX(i)X(i)TW。要使所有的样本的投影方差和最大，也就是最大化∑WTX(i)X(i)TW，即
arg max tr(WTXXTW) s.t.wTw =I
4.3参数求解
拉格朗日乘数法:
J(W)= tr(WTXXTW＋λ(WTw -I))
∂J(W)/Wi = 0
xXTw= ( -入) W
可以看出W是特征向量。
tr(WTXXTW)=>入
因此W=(w(1) , w(2),.… , w(d))，其中w(i)是XXT第i大的特征值对应的特征向量。

4.4算法实现
给定样本集和低维空间的维数

对所有的样本进行中心化操作X=(x(1) , x(2),… . , x(m)：
计算样本均值 u=1/m∑x(i)
所有样本减去均值 x(i)= x(i) - u
计算样本的协方差矩阵 XTX
对协方差矩阵XTX 进行特征值分解
取最大的 d个特征值对应的单位特征向量 w(1) , w(2) ,-- , w(n)，构造投影矩阵W=(w(1) , w(2) ,-- , w(d))

6 结论
PCA方法寻找的是用来有效表示同一类样本共同特点的主轴方向，这对于表示同一类数据样本的共同特征是非常有效的。但PCA不适合用于区分不同的样本类。这点和GMM有很大不同。
PCA算法中舍弃了 n-d个最小的特征值对应的特征向量，用交大的 d特征值对应的特征向量来表示数据，称其为主成分。PCA在数据压缩消除冗余和数据噪音消除等领域都有广泛的应用，如图片压缩。

#original_pca
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt


w  = 3

def norm_(x): #去均值，标准化
    rows, cols = data.shape 
    data_mean = np.sum(data, 0) / rows
    c_data = data - data_mean # 中心化
    return c_data,data_mean

dataset =  np.loadtxt('data_pca.csv', delimiter=',')
data = dataset[:, 0:w]
print("原始数据：")
print(data)
data_new,data_mean = norm_(data)
#print("去均值，标准化")
#print(data_new)
ew, ev = np.linalg.eig(np.cov(data_new.T))
ew_order = np.argsort(ew)[::-1]#排序由大到小
ew_sort = ew[ew_order]
ev_sort = ev[:,ew_order]
#print("ew和ev")
print(ew_sort)
print(ev_sort)
V = ev_sort[:,:2]#取前两个最大的，降维到2
X_new = data_new.dot(V)
print("降维后的数据")
print(X_new)

#picture_pca
import os
import matplotlib.image as mpimg
import cv2
from PIL import Image
import math
import matplotlib.pyplot as plt
from cv2 import cv2 
import numpy as np
size = (40, 40) # 由于较大的数据在求解特征值和特征向量时很慢，故统一压缩图像为size大小


def norm_(x): #去均值，标准化
    rows, cols = data.shape 
    data_mean = np.sum(data, 0) / rows
    c_data = data - data_mean # 中心化
    return c_data,data_mean

'''
    对数据data用PCA降至k维
    data.shape = (N, D)
    返回值：
    X_new降维后的矩阵
'''
def PCA(data,k):
    data_new,data_mean = norm_(data)
    ew, ev = np.linalg.eig(np.cov(data_new.T))
    ew_order = np.argsort(ew)[::-1]
    ew_sort = ew[ew_order]
    ev_sort = ev[:,ew_order]
    V = ev_sort[:,:k]#取前两个最大的，降维到2
    return data_new,V,data_mean

'''
    从file_path中读取面部图像数据
'''
def read_faces(file_path):
    file_list = os.listdir(file_path)
    data = []
    i = 1
    plt.figure(figsize=size)
    for file in file_list:
        path = os.path.join(file_path, file)
        plt.subplot(3, 3, i)
        with open(path) as f:
            img = cv2.imread(path) # 读取图像
            img = cv2.resize(img, size) # 压缩图像至size大小
            img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) # 三通道转换为灰度图
            plt.imshow(img_gray) # 预览
            h, w = img_gray.shape
            img_col = img_gray.reshape(h * w) # 对(h,w)的图像数据拉平
            data.append(img_col)
        i += 1
    plt.show()
    return np.array(data)


'''
    计算峰值信噪比psnr
'''
def psnr(img1, img2):
   mse = np.mean((img1 / 255. - img2 / 255.) ** 2 )
   if mse < 1.0e-10:
      return 100
   PIXEL_MAX = 1
   return 20 * math.log10(PIXEL_MAX / math.sqrt(mse))


data = read_faces('test1')
k=1
row, Column = data.shape
print('降维到',k)
data_new, V,data_mean= PCA(data, k) # PCA降维
V = np.real(V) # 一旦降维维度超过某个值，特征向量矩阵将出现复向量，对其保留实部
pca_data = np.dot(data_new, V) # 计算降维后的数据
recon_data = np.dot(pca_data, V.T)+data_mean#重构数据
plt.figure(figsize=size)
for i in range(row):
    plt.subplot(3,3,i+1)
    plt.imshow(recon_data[i].reshape(size))
plt.show()


print("信噪比如下：")
for i in range(row):
    a = psnr(data[i], recon_data[i])
    print('图', i, '的信噪比: ', a)

#Swiss_Roll_pca
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


'''
    生成瑞士卷数据
    n_sample，生成的数据点数量，default=100
    noise，生成数据点的噪声程度，default=0.0
    y_scale，瑞士卷的厚度，default=100
'''
def roll(n_sample, noise, y_scale):
    t = 1.5 * np.pi * (1 + 2 * np.random.rand(1, n_sample))
    x = t * np.cos(t)
    y = y_scale * np.random.rand(1, n_sample)
    z = t * np.sin(t)
    X = np.concatenate((x, y, z))
    X += noise * np.random.randn(3, n_sample)
    X = X.T
    return X

def show_3D(X):
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.view_init(elev=20, azim=80)
    ax.scatter(X[:,0], X[:,1], X[:,2], c=X[:,0])
    ax.legend(loc='best')
    plt.show()


def show_2D(X):
    plt.scatter(X[:,0], X[:,1], c=X[:,0])
    plt.show()


def norm_(data): #去均值，标准化
    rows, cols = data.shape 
    data_mean = np.sum(data, 0) / rows
    data_new = data - data_mean # 中心化
    return data_new,data_mean

'''
    对数据data用PCA降至k维
    data.shape = (N, D)
    返回值：
    X_new降维后的矩阵
'''
def PCA(data,k):
    #print("原始数据：")
    #print(data)
    data_new,data_mean = norm_(data)
    #print("去均值，标准化")
    #print(data_new)
    ew, ev = np.linalg.eig(np.cov(data_new.T))
    ew_order = np.argsort(ew)[::-1]
    ew_sort = ew[ew_order]
    ev_sort = ev[:,ew_order]
    V = ev_sort[:,:k]#取前两个最大的，降维到k
    X_new = data_new.dot(V)
    return X_new

'''
    人工数据PCA实验
'''


X1 = roll(2000, 0, 10)
show_3D(X1)
pca_data= PCA(X1, 2)
#print(pca_data)
show_2D(pca_data)