EM算法对人脸数据降维(机器学习作业06)

第一题

第二题

代码如下：

import numpy as np
import os
from PIL import Image
from scipy.linalg import sqrtm


def loadFile(filepath):
    sample_list = np.zeros((0, 112 * 92))
    for root, dirs, files in os.walk(filepath):
        for file in files:
            if os.path.splitext(file)[1] == '.pgm':
                # print(os.path.join(root, file))
                im = Image.open(os.path.join(root, file))
                im = np.array(im).flatten()  # 展开数组
                sample_list = np.vstack((sample_list, im))  # 将数组水平拼接
    X = sample_list.T

    return X


def SVD_PCA(X, k=8):
    # 数据中心化
    x_mean = np.sum(X, axis=1) / X.shape[1]
    X = X - x_mean[:, np.newaxis]
    p = X.shape[0]  # 原始特征数量
    m = X.shape[1]  # 样本个数
    # 先求解XTX的协方差矩阵
    c = np.dot(X.T, X)  # 协方差矩阵
    # 求解协方差矩阵的特征向量和特征值
    eigenvalue, featurevector = np.linalg.eig(c)
    # 对特征值索引排序 从大到小
    aso = np.argsort(eigenvalue)
    indexs = aso[::-1]
    # print("特征值:", eigenvalue)
    # print("特征向量:", featurevector)
    # print("降为", k, "维")
    eigenvalue_sum = np.sum(eigenvalue)
    W = []
    for i in range(k):
        # print("第", indexs[i], "特征的解释率为:", (eigenvalue[indexs[i]] / eigenvalue_sum))
        W.append(np.dot(X, featurevector[:, indexs[i]]) / np.sqrt(eigenvalue[indexs[i]]))  # 取前k个特征值大的特征向量作为基向量
    W = np.array(W).T
    # print('W.shape:',W.shape)
    X_trans = np.dot(W.T, X)

    return X_trans


def MLE_PCA(X, k=8):
    # 数据中心化
    x_mean = np.sum(X, axis=1) / X.shape[1]
    X = X - x_mean[:, np.newaxis]
    p = X.shape[0]  # 原始特征数量
    m = X.shape[1]  # 样本个数
    # 先求解XXT的协方差矩阵特征向量与特征值
    c = (1 / m) * np.dot(X, X.T)  # 协方差矩阵
    # 求解协方差矩阵的特征向量和特征值
    eigenvalue, featurevector = np.linalg.eig(c)
    # 对特征值索引排序 从大到小
    aso = np.argsort(eigenvalue)
    indexs = aso[::-1]  # 特征值从大到小的索引列表
    eigenvalue_sum = np.sum(eigenvalue)
    U = []
    A = []
    for i in range(k):
        # print("第", indexs[i], "特征的解释率为:", (eigenvalue[indexs[i]] / eigenvalue_sum))
        U.append(featurevector[:, indexs[i]])  # 取前k个特征值大的特征向量作为基向量
        A.append(eigenvalue[indexs[i]])  # 保存对应特征值
    U = np.array(U).T
    A = np.diag(A)  # 将特征值列表变为对应对角矩阵
    # 计算σ sigma2
    sigma2 = 0
    for j in indexs[k:]:
        sigma2 = j + sigma2
    sigma2 = 1 / (p - k) * sigma2
    # 计算Wml
    W = np.dot(U, np.sqrt(A - sigma2 * np.eye(A.shape[0])))
    # 计算z
    Z = []
    for i in range(m):
        zi = np.dot(np.dot(np.linalg.inv(np.dot(W.T, W) + sigma2 * np.eye(k)), W.T), X[:, i])
        Z.append(zi)
    Z = np.array(Z).T

    return Z


def EM_PCA(X, k=8, iter_num=20):
    # 数据中心化
    x_mean = np.sum(X, axis=1) / X.shape[1]
    X = X - x_mean[:, np.newaxis]
    W = np.random.random([X.shape[0], k])
    for i in range(iter_num):
        # E步
        Z = np.dot(np.dot(np.linalg.inv(np.dot(W.T, W)), W.T), X)
        # M步
        W = np.dot(np.dot(X, Z.T), np.linalg.inv(np.dot(Z, Z.T)))

    return Z


if __name__ == '__main__':
    filepath = r"orl_faces"
    X = loadFile(filepath)
    X_trans = SVD_PCA(X, 8)
    print("SVD_PCA_shape:", X_trans.shape)
    Z = MLE_PCA(X, 8)
    print("MLE_PCA_shape:", Z.shape)
    Z1 = EM_PCA(X)
    print("EM_PCA_shape:", Z1.shape)

结果如下：

MLE_PCA(X, 8)
print("MLE_PCA_shape:", Z.shape)
Z1 = EM_PCA(X)
print("EM_PCA_shape:", Z1.shape)

```

结果如下：