主成分分析（PCA)是一种比较基础的数据降维方法，也是多元统计中的重要部分，在数据分析、机器学习等方面具有广泛应用。主成分分析目的是用较少的变量来代替原来较多的变量，并可以反映原来多个变量的大部分信息。也就是说，在一堆杂乱的数据往往存在相关性，而这些相关性便意味着可以进行数据的降维。对于多维度数据，我们在使用时无法直观的看出数据之间的差距，使用PCA降维后可以更少的维度表示，从而使得表示结果更加直观，减少数据量。在很多情况下数据的降维是十分必要的，一方面有利于问题的简化，另一方面便于计算机的计算：数据降维后变量的减少，会使计算机处理的数据大大减少，从而缩短数据处理时间。

主成分分析原理（PCA）

主成分分析的直观理解，可以认为是旋转坐标轴，使得在旋转坐标轴后这些点在新的坐标系下在各个坐标轴（变量）方向投影的方差变大。其中如果在某坐标上的方差最大，那么这个坐标轴对应的这些散点的坐标就是第一主成分，其次就是第二主成分，依此类推。

对于下图的情况，我们发现这些数据都几乎排列在一条直线上，并且在x轴方向和y轴方向的方差都比较大。但是如果把坐标轴旋转一定角度，使得这些数据在某个坐标轴的投影的方差比较大，便可以用新坐标系下方差较大的一个坐标轴坐标作为主成分。

20241121080033

对于左图，数据为（1，2）、（2，4）……旋转坐标轴后，坐标为$（\sqrt{5},0)$、$（2\sqrt{5},0）$……这样主成分就是新坐标系下变量x的数值：$\sqrt{5}$, $2\sqrt{5}$,$3\sqrt{5}$,…。

PCA求解步骤

输入$m$个样本，特征数为$n$的数据集合：$X = {x_1, x_2, …, x_m}$降维到$k$维。

记样本集为矩阵$X$：
$$
X = \begin{bmatrix}
x_{11} & x_{12} & \dots & x_{1n}\
x_{21} & x_{22} & \dots & x_{2n}\
\dots & \dots & \dots & \dots\
x_{m1} & x_{m2} & \dots & x_{mn}\
\end{bmatrix}
$$
其中每一行代表一个样本，每一列代表一个特征，列号表示特征的维度，共$n$维。

对矩阵去中心化得到新矩阵$X$，即每一列进行零均值化，也即减去这一列的均值$\bar{x_i}$:
$$
\bar{x_i} = \frac{1}{m}\sum_{j=1}^m x_{ji} (i = 1,2,\dots,n)
$$
所求矩阵$X$仍为$m \times n$阶矩阵：
$$
X = \begin{bmatrix}
x_{11}-\bar{x_1} & x_{12}-\bar{x_2} & \dots & x_{1n}-\bar{x_n}\
x_{21}-\bar{x_1} & x_{22}-\bar{x_2} & \dots & x_{2n}-\bar{x_n}\
\dots & \dots & \dots & \dots\
x_{m1}-\bar{x_1} & x_{m2}-\bar{x_2} & \dots & x_{mn}-\bar{x_n}\
\end{bmatrix}
$$
计算去中心化的矩阵X的协方差矩阵:
$$
C = \frac{1}{m-1} X^TX
$$
该矩阵为$n \times n$阶矩阵。
对协方差矩阵$C$进行特征分解，求出协方差矩阵的特征值$\lambda_k$，及对应的特征向量$v_k$:
$$
Cv_k=\lambda_k v_k
$$
将特征向量按对应特征值从左到右按列降序排列成矩阵，取前$k$列组成矩阵$W$，即$n \times k$阶矩阵。
通过$Y = XW$计算降维到$k$维后的样本特征，即$m \times k$阶矩阵。

输出降维后的样本集：$Y = {y_1, y_2, y_3, \dots, y_m}$

PCA python实现

import numpy as np
import matplotlib.pyplot as plt
 
X=np.empty((100,2))
X[:,0]=np.random.uniform(0,100,size=100)
X[:,1]=0.6*X[:,0]+3+np.random.normal(0,10,size=100)
plt.scatter(X[:,0],X[:,1])
 
#对数据进行中心化处理：减去每个特征的均值
def demean(X):
    return X-np.mean(X,axis=0)
X_demean=demean(X)   #保存中心化后的数据
plt.figure(2)
plt.scatter(X_demean[:,0],X_demean[:,1])
 
#定义目标函数：数据在投影到方向 w 上时的方差
def f(w,X):
    return np.sum((X.dot(w)**2))/len(X)
 
#目标函数关于参数 w 的梯度
def df_math(w,X):
    return X.T.dot(X.dot(w))*2/len(X)
 
def direction(w):
    return w / np.linalg.norm(w)
    
def gradient_ascent(df, X, initial_w, eta, n_iters = 1e4, epsilon=1e-8):
    
    w = direction(initial_w) 
    cur_iter = 0
 
    while cur_iter < n_iters:
        gradient = df(w, X)
        last_w = w
        w = w + eta * gradient
        w = direction(w) # 注意1：每次求一个单位方向
        if(abs(f(w, X) - f(last_w, X)) < epsilon):
            break
            
        cur_iter += 1
 
    return w
 
initial_w = np.random.random(X.shape[1]) # 注意2：不能用0向量开始
eta = 0.001
w = gradient_ascent(df_math, X_demean, initial_w, eta)
plt.figure(3)
plt.scatter(X_demean[:,0], X_demean[:,1])
plt.plot([0, w[0]*50], [0 , w[1]*50], color='r')