Dimensionality Reduction —降维

预备知识

1. Data:
   $X = (x_1,x_2,...,x_N)^T_{N*M} = \left(\begin{matrix} x_1^T\\x_2^T\\...\\ x_N^T\end{matrix}\right) = \left(\begin{matrix} x_{11}x_{12}...x_{1M}\\x_{21}x_{22}...x_{2M}\\...\\x_{N1}x_{N2}...x_{NM} \end{matrix}\right)$X=(x1​,x2​,...,xN​)N∗MT​=⎝⎜⎜⎛​x1T​x2T​...xNT​​⎠⎟⎟⎞​=⎝⎜⎜⎛​x11​x12​...x1M​x21​x22​...x2M​...xN1​xN2​...xNM​​⎠⎟⎟⎞​

2. Sample Mean:
   $\overline{X}_{M×1} = \frac{1}{N}\sum_{i = 1}^{N}x_i$XM×1​=N1​i=1∑N​xi​

3. Sample Covariance:
   $S_{M×M} = \frac{1}{N}\sum_{i=1}^{N}(x_i-\overline{x})(x_i-\overline{x})^T$SM×M​=N1​i=1∑N​(xi​−x)(xi​−x)T
   注解：N为特征个数，或样本个数。M表示维度数。 样本均值表示每个维度求均值，因此结果为M×1维。方差转置放在后一项因为假设Xi为列向量（从Data开始）。

S表示各个维度之间的协方差。即因此S矩阵为M×M维

下面将均值和方差矩阵化：
$\overline{X}_{M×1} = \frac{1}{N}\sum_{i = 1}^{N}x_i = \frac{1}{N}(\begin{matrix} x_1 &x_2&x_3&...&x_N\end{matrix})\left(\begin{matrix}1\\1\\1\\.\\.\\.\\1\end{matrix}\right) =\frac{1}{N}X^T1_N$XM×1​=N1​i=1∑N​xi​=N1​(x1​​x2​​x3​​...​xN​​)⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛​111...1​⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞​=N1​XT1N​

$S_{M×M} = \frac{1}{N}\sum_{i=1}^{N}(x_i-\overline{x})(x_i-\overline{x})^T = \frac{1}{N}X^THH^TX = \frac{1}{N}X^THX$SM×M​=N1​i=1∑N​(xi​−x)(xi​−x)T=N1​XTHHTX=N1​XTHX

$其中：H_N = I_N - \frac{1}{N}1_N1_N^T$其中：HN​=IN​−N1​1N​1NT​

写程序就用这个，少用几个for

PCA(Principal Component Analysis) 主成分分析

一个中心，两个基本点:

一个中心：对原始特征空间的重构 （把一组线性相关的变量变换为线性无关的变量）

我的理解：基向量的变换。

两个基本点：
（1）最大投影方差

（2）最小重构距离

两个方法都是达到一个中心的目的，本质相同

最大投影方差

向量投影公式 ：
$令投影方向|\vec{b}| = 1,则\vec{a}在\vec{b}上的投影可以表示为\vec{a}\vec{b},矩阵表示为a^Tb$令投影方向∣b∣=1,则a在b上的投影可以表示为ab,矩阵表示为aTb
我们在计算之前，一般要中心化 即
$对于每一个x_i,计算x_i-\overline{x}$对于每一个xi​,计算xi​−x
假设我们已经知道主轴（主成分）为u1，且|u1| = 1，那么将所有样本投影到主轴上 ，即
$（x_i-\overline{x}）^Tu_1$（xi​−x）Tu1​
因为我们已经中心化，因此投影之后的方差（最大投影方差）为：
$\begin{aligned}J &= \frac{1}{N}\sum_{i=1}^N((x_i-\overline{x})^Tu_1)^2 \\&= \sum_{i=1}^N\frac{1}{N}u_1^T(x_i-\overline{x})(x_i-\overline{x})^Tu_1 \\&= u_1^T(\sum_{i=1}^N\frac{1}{N}(x_i-\overline{x})(x_i-\overline{x})^T)u_1 \\&= u_1^TSu_1\end{aligned}$J​=N1​i=1∑N​((xi​−x)Tu1​)2=i=1∑N​N1​u1T​(xi​−x)(xi​−x)Tu1​=u1T​(i=1∑N​N1​(xi​−x)(xi​−x)T)u1​=u1T​Su1​​
因此，求最大投影方差就变换为求解最优化问题:
$\begin{aligned} \hat{u_1} &= argmax(u_1^TSu_1)\\ s.t.&\space\space u_1^Tu_1 = 1 \end{aligned}$u1​^​s.t.​=argmax(u1T​Su1​)  u1T​u1​=1​
通过拉格朗日方程求解。

很显然，求值的结果是协方差矩阵S的特征值分解。

这个过程我们可以理解为

1.特征空间变换
2.通过特征值 筛选特征向量
3 .降维

最小重构代价

假设ui为变换之后的特征向量，则无损失重构为：
$x_i = \sum_{k=1}^M(x_i^Tu_k)u_k$xi​=k=1∑M​(xiT​uk​)uk​
如果对变换后的特征向量进行降维处理(降为q维)，则重构之后的损失为:
$\hat{x}_i=\sum_{k=1}^q(x_i^Tu_k)u_k$x^i​=k=1∑q​(xiT​uk​)uk​

$J = \frac{1}{N}\sum_{i=1}^N||x_i-\hat{x}_i||^2= \frac{1}{N}\sum_{i=1}^N\sum_{k=q+1}^M(x_i^Tu_k)^2$J=N1​i=1∑N​∣∣xi​−x^i​∣∣2=N1​i=1∑N​k=q+1∑M​(xiT​uk​)2

因为实际中我们要对xi进行中心化，因此上式中
$x_i = x_i- \overline{x_i}$xi​=xi​−xi​​
替换上式子中的Xi，则有
$\frac{1}{N}\sum_{i=1}^N\sum_{k=q+1}^M(x_i^Tu_k)^2 = \frac{1}{N}\sum_{i=1}^N\sum_{k=q+1}^M((x_i- \overline{x_i})^Tu_k)^2 =\sum_{k=q+1}^M\sum_{i=1}^N\frac{1}{N}((x_i- \overline{x_i})^Tu_k)^2 =\sum_{k=q+1}^Mu_k^TSu_k$N1​i=1∑N​k=q+1∑M​(xiT​uk​)2=N1​i=1∑N​k=q+1∑M​((xi​−xi​​)Tuk​)2=k=q+1∑M​i=1∑N​N1​((xi​−xi​​)Tuk​)2=k=q+1∑M​ukT​Suk​
问题变为：
$u_k=argmin\sum_{k=q+1}^Mu_k^TSu_k\\ s.t.\space u_k^Tu_k=1$uk​=argmink=q+1∑M​ukT​Suk​s.t. ukT​uk​=1
解拉格朗日方程，就可看出与最大投影方差相同。

从SVD（Singular Value Decomposition）奇异值分解的角度分析PCA

我们从样本的角度分析

首先对样本进行中心化处理：
$X-\overline{X}=HX=U\sum V^T\\ S=X^THX=X^TH^THX=(HX)^T(HX)=V\sum U^TU\sum V^T=V\sum\sum V^T$X−X=HX=U∑VTS=XTHX=XTHTHX=(HX)T(HX)=V∑UTU∑VT=V∑∑VT
因此，我们可以不求协方差矩阵，直接求X的奇异值分解即可
$T = HXX^TH=U\sum \sum U^T$T=HXXTH=U∑∑UT
若我们直接用S做，可以直接得到坐标和基向量，若我们用T做，则不能得到基向量，但是可以得到坐标。用T做我们叫它主坐标分析（Principle Coordinate Analysis）PCoA

S一般为M×M，T一般为N×N,选择合适的角度进行分解，可以简化计算量。

从概率的角度分析PCA，P-PCA

https://www.bilibili.com/video/av32709936/?p=6

代码实现

PCA

import numpy as np
from sklearn.datasets import load_iris
import matplotlib.pyplot as plt
import matplotlib
plt.rcParams['savefig.dpi'] = 150 #图片像素
plt.rcParams['figure.dpi'] = 150 #分辨率
#初始化鸢尾花数据
iris = load_iris()
X,y = iris.data , iris.target
class1_X = X[:100,:2]       # ['sepal length (cm)', 'sepal width (cm)']
class2_X = X[:100,2:]       # ['petal length (cm)', 'petal width (cm)']

X = class1_X[:50,:]
x_overline = np.sum(X,axis = 0)/X.shape[0]
a = X - x_overline
b = a.T.dot(a)
V,sigma,VT = np.linalg.svd(b,1,1)
mark = np.argmax(sigma)
V = V[:,mark] 
print(V)
result = X.dot(V.T).reshape([-1,1])*V.reshape([1,-1])
print(X.dot(V.T))

plt.title("test")
plt.xlabel("sepal length")
plt.ylabel("sepal width")
plt.plot(result[:50,0],result[:50,1],'bo',color = 'blue',label = '0')
plt.plot(X[:50,0],X[:50,1],'bo',color = 'orange',label = '1') 

#下面是fisher判别，因为想到这个了，就写一些
X1 = class1_X[:50,:]
X2 =class1_X[50:,:]
x1_overline = np.sum(X1,axis = 0)/X1.shape[0]
x2_overline = np.sum(X2,axis = 0)/X2.shape[0] 


In = np.identity(X1.shape[0])
Hn = In - 1/X1.shape[0]*np.ones([X1.shape[0],X1.shape[0]])
S1 = 1 / X1.shape[0] * X1.T.dot(Hn).dot(X1)
S2 = 1 / X2.shape[0] * X2.T.dot(Hn).dot(X2)
Sw = S1 + S2
w = np.linalg.inv(Sw).dot(x1_overline-x2_overline)
m = (x1_overline + x2_overline) / 2
m2 = m.dot(w.T).reshape([-1,1])*w
print(m2[:,1])
result1 = X1.dot(w.T).reshape([-1,1])*w
result2 = X2.dot(w.T).reshape([-1,1])*w

plt.title("PCA")
plt.xlabel("sepal length")
plt.ylabel("sepal width")
plt.plot(class1_X[:50,0],class1_X[:50,1],'bo',color = 'blue',label = '0')
plt.plot(class1_X[50:,0],class1_X[50:,1],'bo',color = 'orange',label = '1') 
plt.plot(result1[:50,0],result1[:50,1],'bo',color = 'blue',label = '0')
plt.plot(result2[:50,0],result2[:50,1],'bo',color = 'orange',label = '0')
plt.plot(m2[:,0],m2[:,1],'bo',color = 'red',label = '1') 

参考视频

https://www.bilibili.com/video/av32709936/?p=6