【python】双向二维PCA（2D-2D PCA）算法实现

最新推荐文章于 2025-02-25 13:35:23 发布

yaotouge

最新推荐文章于 2025-02-25 13:35:23 发布

阅读量7.8k

点赞数 12

分类专栏： python CV 文章标签： python 算法图片

本文链接：https://blog.youkuaiyun.com/w450468524/article/details/54895477

版权

python 同时被 2 个专栏收录

2 篇文章

订阅专栏

2 篇文章

订阅专栏

本文深入讲解了2D-2DPCA的原理及其在图像处理中的应用，对比了传统PCA与2DPCA的方法，介绍了如何通过特征向量降低图像数据的维度，并给出了Python实现。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

原理介绍

PCA

这篇教程讲的非常详细。

以图片为例，将每张图片（假设宽高为 m n）转化为m*n的一维向量，然后减去每张图片像素的均值（使得每张图片均值为零），然后数据集中的所有图片一起合成一个大的矩阵X。

计算图片矩阵的方差-协方差矩阵，我们的优化目标就是将这个矩阵对角化，使得字段（代表一张图片）两两之间的协方差为零，也就是相关性最小，字段与自身的方差尽可能大。如下：
$D\space =\frac{1}{m}YY^T \\ \quad =\frac{1}{m}(PX)(PX)^T \\ \quad =\frac{1}{m}PXX^TP^T \\ \quad =P(\frac{1}{m}XX^T)P^T \\ \quad = PCP^T$

Y=PX，P为一组基，Y为在P上映射后的数据，D为映射后数据的协方差矩阵，也就是我们的优化目标矩阵，X为数据矩阵。 $C=\frac{1}{m} X X^T$ ，也就是原数据的协方差矩阵。我们的目标现在是找到能将矩阵C对角化，使得D成为一个对角矩阵。

数学上矩阵对角化只需计算出C的特征值和特征向量，并取前k大的特征值对应的特征向量，也就是投影之后方差最大的那些基，即可得到P矩阵，k就是降维后的维数。

2D PCA

相较于传统PCA，2D PCA不用将图片转换为1D向量，极大的减少了数据的维度。基本思路和PCA相同，也是将数据投影到某组基上，然后使得投影之后数据的协方差矩阵成为对角阵。投影后数据的协方差矩阵为：

$S_x=E[(Y-EY)(Y-EY)^T]=E[AX-E(AX)][AX-E(AX)]^T \\ \quad = E[(A-EA)X][(A-EA)X]^T$

Y为投影之后的数据，Y=AX，A是原数据，X是一组正交基，Sx也就是协方差矩阵。上述公式变形之后之后可得（利用tr(AB)=tr(BA)，只要保证迹不变即可，这里还有点没搞懂为什么）：

$S_x=X^T[E(A-EA)^T(A-EA)]X$

然后再定义矩阵Gt为：

$G_t=E(A-EA)^T(A-EA)$
也就是原数据的协方差矩阵。

类似于传统PCA，我们优化目标就是将使矩阵S_x成为对角阵，并找出对应的特征向量和特征值，也就是将G_t对角化，取前K大的特征向量。

2D-2DPCA

双向2D PCA的改进之处在于，普通2D PCA中的变换只提取了数据矩阵行内的特征，对行进行了变换，而双向PCA则加上了对列进行的变换。

$S_x=E[(B-EB)(B-EB)^T]=E[Z^TA-E(Z^TX)][Z^TA-E(Z^TA)]^T \\ \quad = Z^TE[(A-EA)(A-EA)^T]Z$

列变换对应的Gt矩阵为：

$G_t=E[(A-EA)(A-EA)^T]$

对列进行2D PCA也就是对上述矩阵提取特征向量，得到一组基Z，然后对原始数据进行投影。

为了将列和行的变换合并，最后进行如下变换：

$C=Z^TAX$

X为行方向上2D PCA得到的变换矩阵，Z为列方向上2D PCA得到的变换矩阵，A为原数据。要还原图片数据也很简单，已知C矩阵：

$A=ZCX^T$

因为特征向量相互正交，Z，X均为正交阵，正交阵的逆矩阵就是其转置。

python 实现

用python实现了一下算法：

# a implementation of 2D^2 PCA algorithm

import numpy as np
from PIL import Image

def PCA2D_2D(samples, row_top, col_top):
    '''samples are 2d matrices'''
    size = samples[0].shape
    # m*n matrix
    mean = np.zeros(size)

    for s in samples:
        mean = mean + s

    # get the mean of all samples
    mean /= float(len(samples))

    # n*n matrix
    cov_row = np.zeros((size[1],size[1]))
    for s in samples:
        diff = s - mean;
        cov_row = cov_row + np.dot(diff.T, diff)
    cov_row /= float(len(samples))
    row_eval, row_evec = np.linalg.eig(cov_row)
    # select the top t evals
    sorted_index = np.argsort(row_eval)
    # using slice operation to reverse
    X = row_evec[:,sorted_index[:-row_top-1 : -1]]

    # m*m matrix
    cov_col = np.zeros((size[0], size[0]))
    for s in samples:
        diff = s - mean;
        cov_col += np.dot(diff,diff.T)
    cov_col /= float(len(samples))
    col_eval, col_evec = np.linalg.eig(cov_col)
    sorted_index = np.argsort(col_eval)
    Z = col_evec[:,sorted_index[:-col_top-1 : -1]]

    return X, Z


samples = []
for i in range(1,6):
    im = Image.open('image/'+str(i)+'.png')
    im_data  = np.empty((im.size[1], im.size[0]))
    for j in range(im.size[1]):
        for k in range(im.size[0]):
            R = im.getpixel((k, j))
            im_data[j,k] = R/255.0
    samples.append(im_data)

X, Z = PCA2D_2D(samples, 90, 90)

res = np.dot(Z.T, np.dot(samples[0], X))
res = np.dot(Z, np.dot(res, X.T))

row_im = Image.new('L', (res.shape[1], res.shape[0]))
y=res.reshape(1, res.shape[0]*res.shape[1])

row_im.putdata([int(t*255) for t in y[0].tolist()])
row_im.save('X.png')