线性与非线性数据降维方法汇总（Python代码实现）

作者简介：大家好，我是车神哥，府学路18号的车神磊

个人主页：应无所住而生其心的博客_府学路18号车神_CSDN博客
点赞➕评论➕收藏 == 养成习惯（一键三连）

⚡希望大家多多支持珞~一起加油 

专栏《Fault Diagnosis》

其他专栏：

《LeetCode天梯》
《Neural Network》
《Python》
《Algorithm》

数据降维方法及Python实现

前言
大纲
- 线性数据降维
- - 主元分析（PCA）
  - 偏最小二乘法（PLS）
  - 独立成分分析（ICA）
  - 线性判别分析（LDA）
  - 典型相关分析（CCA）
  - 慢特征分析（SFA）
- 非线性数据降维
- - 神经网络非线性方法
  - 核（Kernel）方法
  - - 线性核函数
    - 高斯核函数
  - 流行学习方法
  - - LLE(Locally Linear Embedding) 局部线性嵌入
    - LE (Laplacian Eigenmaps) 拉普拉斯特征映射
    - LPP(Locality Preserving Projection) 局部保留投影

前言

11月份了，目前再准备开题，依旧搞故障诊断，在工业过程中的过程监控现状大多是从数据驱动的角度进行监控和诊断，然鹅，我们遇到的很多工况数据确实非线性的，在实际生产中，大部分都是非线性的数据，但也有一些线性的，由此我们可以从全局和局部出发，来对数据进行挖掘和分析，对此，近期我做了很多关于数据降维的一些工作用于毕业论文，现在在下面讲讲我的一些见解和Python代码的复现。

Python代码在文末哟！~

大纲

线性数据降维

线性降维方法，如主元分析（PCA）＂、部分最小二乘化（PLS）及独立元分析（ICA），线性辨别分析（LDA）等。

主元分析（PCA）

之前有过一篇blog介绍过这一部分，-> 看这里

过程监控中的应用可看 -> 这里

还是在介绍一遍吧~

PCA是一种统计方法，广泛应用于工程和科学应用中，与傅里叶分析相比，尤其适用于质量监测。

设 x ∈ R m boldsymbol{x} in mathfrak{R}^{m} x∈Rm表示 m m m个传感器矢量的样本测量值。

假设每个传感器有 N N N个样本，数据矩阵 X = [ x 1 x 2 ⋯ x N ] T ∈ R N × m mathbf{X}=left[begin{array}{llll} boldsymbol{x}_{1} & boldsymbol{x}_{2} & cdots & boldsymbol{x}_{N} end{array}right]^{T} in mathfrak{R}^{N times m} X=[x1x2⋯xN]T∈RN×m，由代表样本 x i T x^T_i xiT的每一行组成。

正常数据矩阵 X X X的一个重要要求是，它应具有丰富的正常变化，以代表过程的共同原因变化。矩阵 X X X被缩放为零均值，通常为PCA建模的单位方差。矩阵 X X X通过奇异值分解（SVD）分解为得分矩阵 T T T和加载矩阵 P P P，
X = T P T + X ~ (1) mathbf{X}=mathbf{T P}^{T}+tilde{mathbf{X}}tag{1} X=TPT+X~(1)

其中 T = X P T=XP T=XP包含 l l l 个左前导奇异向量和奇异值，P 包含 l l l个右前导奇异向量， X ~ tilde{mathbf{X}} X~ 是残差矩阵。因此，T 的列是正交的，P 的列是正交的。将样本协方差矩阵表示为
S = 1 N − 1 X T X (2) mathbf{S}=frac{1}{N-1} mathbf{X}^{T} mathbf{X}tag{2} S=N−11XTX(2)

作为SVD的替代方法，可以对 S 进行特征分解，以获得 P 作为 S 的 l l l 个前导特征向量，特征值表示为

Λ = diag ⁡ { λ 1 , λ 2 , … , λ l } (3) mathbf{Lambda}=operatorname{diag}left{lambda_{1}, lambda_{2}, ldots, lambda_{l}right}tag{3} Λ=diag{λ1,λ2,…,λl}(3)
第 i i i 个特征值可与得分矩阵 T 的第 i i i 列相关，如下所示：
λ i = 1 N − 1 t i T t i ≈ var ⁡ { t i } (4) lambda_{i}=frac{1}{N-1} mathbf{t}_{i}^{T} mathbf{t}_{i} approx operatorname{var}left{mathbf{t}_{i}right}tag{4} λi=N−11tiTti≈var{ti}(4)

这是第 i i i个得分向量 t i ∈ R N mathbf{t}_{i} in mathfrak{R}^{N} ti∈RN的样本方差。主成分子空间（PCS）是 S p = span ⁡ { P } mathcal{S}_{p}=operatorname{span}{mathbf{P}} Sp=span{P}，剩余子空间（RS） S r S_r Sr是 S p S_p Sp的正交补。将测量空间划分为PCS和RS，使得RS仅包含微小的奇异值，这些奇异值对应于通常具有较小变化的子空间，或者主要是噪声的子空间。因此，残差类似于根据质量平衡和能量平衡建立的数学模型中的方程误差。

样本向量 x ∈ R m mathbf{x} in mathfrak{R}^{m} x∈Rm可以分别投影到PCS和RS上，

x ^ = P t = P P T x ∈ S p (5) hat{boldsymbol{x}}=mathbf{P} boldsymbol{t}=mathbf{P P}^{T} boldsymbol{x} in mathcal{S}_{p}tag{5} x^=Pt=PPTx∈Sp(5)

其中，

t = P T x ∈ R l (6) boldsymbol{t}=mathbf{P}^{T} boldsymbol{x} in mathfrak{R}^{l}tag{6} t=PTx∈Rl(6)

为 l l l 个潜在变量得分的向量。

残差向量：

x ~ = x − x ^ = ( I − P P T ) x ∈ S r (7) tilde{boldsymbol{x}}=boldsymbol{x}-hat{boldsymbol{x}}=left(mathbf{I}-mathbf{P P}^{T}right) boldsymbol{x} in mathcal{S}_{r}tag{7} x~=x−x^=(I−PPT)x∈Sr(7)

因为 S p S_p Sp 和 S r S_r Sr 是正交的，

x ^ T x ~ = 0 (8) hat{boldsymbol{x}}^{T} tilde{boldsymbol{x}}=0tag{8} x^Tx~=0(8)

且

x = x ^ + x ~ (9) boldsymbol{x}=hat{boldsymbol{x}}+tilde{boldsymbol{x}}tag{9} x=x^+x~(9)

其中，一个重要的概念是，数据的PCA模型， x ^ hat{boldsymbol{x}} x^由潜变量 t ∈ R m mathbf{t} in mathfrak{R}^{m} t∈Rm 参数化。

偏最小二乘法（PLS）

PLS具体算法过程如下：

将X和Y经过标准化（包括减均值、除标准差等）（此操作可参考这里的CCA算法）。
设X的第一个主成分为 p 1 p_1 p1，Y的第一个主成分为 q 1 q_1 q1，两者都经过单位化。（注意：这里的主成分并不是通过PCA得到的主成分）（~~菜鸡解释：~~这里主成分可简单的看成CCA系数矩阵中的第一系数成分，如 a 1 a_1 a1）。
u 1 = X p 1 ， v 1 = Y q 1 u_1=Xp_1，v_1=Yq_1 u1=Xp1，v1=Yq1，这和CCA几乎一样，可得下面的期望的约束条件。
V a r ( u 1 ) → M a x ， V a r ( v 1 ) → M a x Var(u_1)rightarrow Max，Var(v_1)rightarrow Max Var(u1)→Max，Var(v1)→Max，代表在主成分分量上的投影，得到了期望的方差最大化值（~~貌似这样解释有点奇怪~~）。
R ( u 1 , v 1 ) → M a x R_(u_1,v_1)rightarrow Max R(u1,v1)→Max，和CCA一样。
综合以上条件，可得 C o v ( u 1 ， v 1 ) = V a r ( u 1 ) V a r ( v 1 ) R ( u 1 , v 1 ) → M a x Cov(u_1，v_1)=sqrt{Var(u_1)Var(v_1)}R_(u_1,v_1) rightarrow Max Cov(u1，v1)=Var(u1)Var(v1) R(u1,v1)→Max。

简而言之，为了实现偏最小二乘回归的基本思想，要求p1和q1的协方差最大，即求解下面优化后的目标函数：

M a x : < X p 1 , Y q 1 > S . t . : ∣ ∣ p 1 ∣ ∣ = 1 , ∣ ∣ q 1 ∣ ∣ = 1 Max: \ S.t.:||p_1||=1,||q_1||=1 Max:S.t.:∣∣p1∣∣=1,∣∣q1∣∣=1

看似比CCA的简单，这里附上CCA的目标函数：

M a x : R ( U , V ) = Cov ⁡ ( U , V ) Var ⁡ [ U ] Var ⁡ [ V ] = C o v ( U , V ) = t k T C o v ( A , B ) h k = t k T Σ 12 h k S . t . : V a r ( U k ) = V a r ( t k T A ) = t k T Σ 11 t k = 1 , V a r ( V k ) = V a r ( h k T A ) = h k T Σ 22 h k = 1 Max:R_{(U,V)}=frac{operatorname{Cov}(U, V)}{sqrt{operatorname{Var}[U] operatorname{Var}[V]}}=Cov(U,V)={t_k}^TCov(A,B)h_k={t_k}^TSigma_{12} h_k\ S.t.:Var(U_k)=Var({t_k^T}{A})={t_k^T}Sigma_{11}t_k=1, Var(V_k)=Var({h_k^T}{A})={h_k^T}Sigma_{22}h_k=1 Max:R(U,V)=Var[U]Var[V] Cov(U,V)=Cov(U,V)=tkTCov(A,B)hk=tkTΣ12hkS.t.:Var(Uk)=Var(tkTA)=tkTΣ11tk=1,Var(Vk)=Var(hkTA)=hkTΣ22hk=1

上面CCA是一次求解的过程，而我们的PLS回归只是对目前的第一主成分做了优化计算，剩下的主成分还得再计算。

关于优化的目标求解的办法，和CCA一样，也是引用了拉格朗日乘数法来求解。（下面给出详细计算步骤）

首先，引入拉格朗日乘子：
L = p 1 T X T Y q 1 − λ 2 ( p 1 T p 1 − 1 ) − θ 2 ( q 1 T q 1 − 1 ) mathcal{L}=p_{1}^{T} X^{T} Y q_{1}-frac{lambda}{2}left(p_{1}^{T} p_{1}-1right)-frac{theta}{2}left(q_{1}^{T} q_{1}-1right) L=p1TXTYq1−2λ(p1Tp1−1)−2θ(q1Tq1−1)
分别对 p 1 p_1 p1和 q 1 q_1 q1求偏导,
∂ L ∂ p 1 = X τ Y q 1 − λ p 1 = 0 ∂ L ∂ q 1 = Y τ X p 1 − θ q 1 = 0 begin{array}{l} frac{partial mathcal{L}}{partial p_{1}}=X^{tau} Y q_{1}-lambda p_{1}=0 \\ frac{partial mathcal{L}}{partial q_{1}}=Y^{tau} X p_{1}-theta q_{1}=0 end{array} ∂p1∂L=XτYq1−λp1=0∂q1∂L=YτXp1−θq1=0
和CCA一样，可求得 λ lambda λ和 θ theta θ相等。

将 λ − 1 X τ Y q 1 = p 1 lambda^{-1}X^{tau} Y q_{1}= p_{1} λ−1XτYq1=p1代入上面第二式子，可得
Y τ X X τ Y q 1 = λ 2 q 1 Y^{tau}XX^{tau} Y q_{1}= lambda^{2} q_{1} YτXXτYq1=λ2q1
两边均乘以 p 1 p_1 p1或 q 1 q_1 q1，再利用约束条件 ∣ ∣ p 1 ∣ ∣ = 1 , ∣ ∣ q 1 ∣ ∣ = 1 ||p_1||=1,||q_1||=1 ∣∣p1∣∣=1,∣∣q1∣∣=1，可得：
X τ Y Y τ X p 1 = λ 2 p 1 X^{tau}YY^{tau} X p_{1}= lambda^{2} p_{1} XτYYτXp1=λ2p1
故上式 λ 2 lambda^{2} λ2则为 X τ Y Y τ X p 1 X^{tau}YY^{tau} X p_{1} XτYYτXp1的特征值， p 1 p_1 p1为相应的单位特征向量， q 1 q_1 q1一样。

求得 p 1 p_1 p1、 q 1 q_1 q1这样 ⟨ X p 1 , Y q 1 ⟩ → p 1 τ X τ Y q 1 → p 1 τ ( λ p 1 ) → λ leftlangle X p_{1}, Y q_{1}rightrangle rightarrow p_{1}^{tau} X^{tau} Y q_{1} rightarrow p_{1}^{tau}left(lambda p_{1}right) rightarrow lambda ⟨Xp1,Yq1⟩→p1τXτYq1→p1τ(λp1)→λ可得到最优解。

可见 p 1 p_1 p1和 q 1 q_1 q1是投影方差最大和两者相关性最大上的权衡，而CCA只是相关性上最大化。

到此，我们可以得到 u 1 、 v 1 u_1、v_1 u1、v1值，这里的 u 1 、 v 1 u_1、v_1 u1、v1在图上面只是表示为绿色点，如果这样就完成求解，那和CCA的过程一样，得不到X到Y的映射。

具体整个算法的详细原理(模型+回归)可看 -> 这里

独立成分分析（ICA）

独立成分分析 ICA(Independent Component Correlation Algorithm)是一种函数，X为n维观测信号矢量，S为独立的m（m<=n)维未知源信号矢量，矩阵A被称为混合矩阵。ICA的目的就是寻找解混矩阵W（A的逆矩阵），然后对X进行线性变换，得到输出向量U。

这里使用最大似然估计来解释算法，我们假定每个 s i s_i si有概率密度 p s p_s ps，那么给定时刻原信号的联合分布就是
p ( s ) = ∏ i = 1 n p s ( s i ) mathrm{p}(mathrm{s})=prod_{i=1}^{n} p_{s}left(s_{i}right) p(s)=i=1∏nps(si)
此公式代表一个假设前提：每个人发出的声音信号各自独立。

有了 p ( s ) p(s) p(s)，我们可以求得 p ( x ) p(x) p(x)
p ( x ) = p s ( H x ) ∣ H ∣ = ∣ H ∣ ∏ i = 1 n p s ( h i T x ) mathrm{p}(mathrm{x})=mathrm{p}_{s}(H x)|mathrm{H}|=|mathrm{H}| prod_{i=1}^{n} p_{s}left(h_{i}{ }^{T} xright) p(x)=ps(Hx)∣H∣=∣H∣i=1∏nps(hiTx)
左边是每个采样信号 x x x的概率，右边是每个原信号概率的乘积的 ∣ H ∣ |H| ∣H∣倍。

若没有先验知识，我们无法求得 H H H和 s s s。

因此我们需要知道 p s ( s i ) p_s(s_i) ps(si)，我们打算选取一个概率密度函数赋给 s s s，但是我们不能选取高斯分布的密度函数。在概率论里我们知道密度函数p(x)由累计分布函数（cdf）F(x)求导得到。F(x)要满足两个性质是：单调递增和在[0,1]。我们发现sigmoid函数很适合，定义域负无穷到正无穷，值域0到1，缓慢递增。我们假定 s s s的累积分布函数符合sigmoid函数
g ( s ) = 1 1 + e − s g(s)=frac{1}{1+e^{-s}} g(s)=1+e−s1
求导可得，
p s ( s ) = g ′ ( s ) = e s ( 1 + e s ) 2 p_{s}(s)=g^{prime}(s)=frac{e^{s}}{left(1+e^{s}right)^{2}} ps(s)=g′(s)=(1+es)2es
这就是 s s s的密度函数。此时的 s s s是实数。

要是我们预先知道 s s s的分布函数，那就不用假设了，但在未知的情况下，sigmoid函数能够在大多数问题上取得不错的效果。

由于上式中 p s ( s ) p_s(s) ps(s)是个对称函数，因此E[s]=0（s的均值为0），那么E[x]=E[As]=0，x的均值也是0。

现在我们知道了 p s ( s ) p_s(s) ps(s)，下面开始求 H H H。

采样后的训练样本为 X ( i ) = ( x 1 ( i ) , x 2 ( i ) , … , x n ( i ) ) ; ( i = 1 , … , m ) mathrm{X}^{(i)}=left(x_{1}^{(i)}, x_{2}^{(i)}, ldots, x_{n}^{(i)}right) ;( i=1, ldots, m) X(i)=(x1(i),x2(i),…,xn(i));(i=1,…,m)，使用前面得到的 x x x的概率密度函数，得其样本对数似然估计：
ℓ ( H ) = ∑ i = 1 m ( ∑ j = 1 n log ⁡ g ′ ( h j T x ( i ) ) + log ⁡ ∣ H ∣ ) ell(H)=sum_{i=1}^{m}left(sum_{j=1}^{n} log g^{prime}left(h_{j}^{T} x^{(i)}right)+log |H|right) ℓ(H)=i=1∑m(j=1∑nlogg′(hjTx(i))+log∣H∣)
其中，括号里的一大堆为 p ( x ( i ) ) p(x^{(i)}) p(x(i))，然后再对 H H H进行求导操作。在上式中包含有行列式，对行列式|W|进行求导的方法可参考这里。

最终得到的求导结果公式（~~很复杂很繁琐~~–心情）：
H : = H + α ( [ 1 − 2 g ( h 1 T x ( i ) ) 1 − 2 g ( h 2 T x ( i ) ) ⋮ 1 − 2 g ( h n T x ( i ) ) ] x ( i ) T + ( H T ) − 1 ) H:=H+alphaleft(left[begin{array}{c} 1-2 gleft(h_{1}^{T} x^{(i)}right) \ 1-2 gleft(h_{2}^{T} x^{(i)}right) \ vdots \ 1-2 gleft(h_{n}^{T} x^{(i)}right) end{array}right] x^{(i)^{T}}+left(H^{T}right)^{-1}right) H:=H+α⎝⎜⎜⎜⎛⎣⎢⎢⎢⎡1−2g(h1Tx(i))1−2g(h2Tx(i))⋮1−2g(hnTx(i))⎦⎥⎥⎥⎤x(i)T+(HT)−1⎠⎟⎟⎟⎞
其中 α alpha α表示的是梯度上升速率，可自定义。

当通过多次迭代后，可求出 H H H，便可得到 s ( i ) = H x ( i ) s^{(i)}=Hx^{(i)} s(i)=Hx(i)来还原出原始信号。

具体独立成分分析ICA原理及应用可看 —> 这里

线性判别分析（LDA）

LDA的思想：由所给定的数据集，设法将样例数据投影在一条直线上，使得同类数据的投影点尽可能的接近、而异类数据的投影点之间将可能间隔更远。在我们做新样本数据的分类时，将其投影到同样的直线上，再根据投影点的位置来确定新样本的类别。如下图（源自周志华《机器学习》）所示：

这里的投影直线也用到了最小二乘的思想，所有数据样本垂直投影在直线上，只是我们的约束条件表为了不同数据样本之间的投影在直线上的距离度量。
我们需要寻找到在投影方向 w w w上，使得数据样本满足两个条件：1) 相同数据之间投影距离最小；2）不同数据之间投影点位置最大（可通过求其不同数据的投影中心点来判别）

LDA二维示意图

图中，“+”和“-”代表的是两种不同的数据簇，而椭圆表示数据簇的外轮廓线，虚线表示其投影，红色实心圆●和红色实心三角形△分别代表的两类数据簇投影到 w w w方向上的中心点。

对于上面投影方向 y = w T x y=mathbf{w^Tx} y=wTx，有博主认为描述的不够准确，书中并未提及关于 y y y的解释，但是对于 y y y其实是有所提及的。

但我认为，这里的 y y y，仅仅是为了体现投影的一个方向，将数据 x x x投影在方向为 w w w的直线上，而不是代表的这根投影直线为 y = w T x y=mathbf{w^Tx} y=wTx，或许会被人误认为是投影后的值0。
—~~菜鸡理解~~（如有不对，请批评指正）

已知给定的数据集为
D = { ( x i , y i ) } i = 1 m ， y i ∈ { 0 , 1 } D={(x_i,y_i)}_{i=1}^{m}，y_iin {0,1} D={(xi,yi)}i=1m，yi∈{0,1}
假设 X i 、 μ i 、 Σ i X_i、mu_i、Sigma_i Xi、μi、Σi分别表示第 i ∈ { 0 , 1 } iin{0,1} i∈{0,1}类（注意：这里的 i i i 指代有多少个不同的类别数据集，图中只有两类，故为0和1）示例的集合、均值向量、协方差矩阵。

假如将所有的样本数据点都投影到直线 w w w上来，那么两类不同的样本数据的中心点在直线上的投影可表示为 w T μ 0 、 w T μ 1 w^{T}mu_0、w^{T}mu_1 wTμ0、wTμ1；同样，所有样本投影到直线上后，我们得到的两类样本的协方差分别为 w T Σ 0 w w^{T}Sigma_0w wTΣ0w和 w T Σ 1 w w^{T}Sigma_1w wTΣ1w.

由于我们只是在一维平面上的直线，故为一维空间，由此 w T μ 0 、 w T μ 1 、 w T Σ 0 w 、 w T Σ 1 w w^{T}mu_0、w^{T}mu_1、w^{T}Sigma_0w、w^{T}Sigma_1w wTμ0、wTμ1、wTΣ0w、wTΣ1w都是实数。

为什么说这里是一维空间呢？可以看上图，假设每个样本都是d维向量（上图为二维 x 1 、 x 2 x_1、x_2 x1、x2坐标系）。现在就简单一点，想用一条直线 w w w表示这些样本，称之为样本集合的一维表达。所以这里说的一维讲的是投影到一条直线上以后的数据，在直线上是属于一维空间表达的。

下面思考另一个问题，如何让同类的数据样本投影点尽可能的靠近，而使得不同样本投影点离得更远呢？

这里需要引入协方差的概念，小小复习一下协方差及样本方差的知识（~~因为本菜鸡数学基础差~~）

协方差（Covariance）在概率论和统计学中用于衡量两个变量的总体误差。而方差是协方差的一种特殊情况，即当两个变量是相同的情况。
上面的 Σ 0 、 Σ 1 Sigma_0、Sigma_1 Σ0、Σ1因为是自协方差也就是代表方差（也即为样本方差）。

方差：当数据分布比较分散（即数据在平均数附近波动较大）时，各个数据与平均数的差的平方和较大，方差就较大；当数据分布比较集中时，各个数据与平均数的差的平方和较小。
总的说来：方差越大，数据的波动越大；方差越小，数据的波动就越小。

协方差表示的是两个变量的总体的误差，这与只表示一个变量误差的方差不同。如果两个变量的变化趋势一致，也就是说如果其中一个大于自身的期望值，另外一个也大于自身的期望值，那么两个变量之间的协方差就是正值。如果两个变量的变化趋势相反，即其中一个大于自身的期望值，另外一个却小于自身的期望值，那么两个变量之间的协方差就是负值。
简而言之：两个变量之间差距越大，协方差就越小；相反，两个变量越相似变化趋势一致，则协方差越大。

复习完协方差、样本方差的知识后，解决上面的问题应该不难。

按照我们的需求，让同类的样本投影点尽可能的靠近，换句话说就是让同类样本投影的协方差尽可能的小（注意：这里由于是自协方差==样本方差，也就满足上面大字第一条），即 w T Σ 0 w + w T Σ 1 w w^{T}Sigma_0 w+w^{T}Sigma_1 w wTΣ0w+wTΣ1w尽可能的小，这样数据的波动就小，之间的距离就更小更靠近。

关于不同数据样本投影点之间的操作，使其更加的远离。我们可以通过不同数据集投影的中心点来判别，不同中心点之间的距离越大，那么表示他们之间离得更远，则 ∣ ∣ w T μ 0 + w T μ 1 ∣ ∣ 2 2 ||w^{T}mu_0+w^{T}mu_1||_2^2 ∣∣wTμ0+wTμ1∣∣22（欧式距离）更大。

好！现在我们同时考虑两者的情况，则可以使得得到最大化的目标，建立我们的模型：
J = ∥ w T μ 0 − w T μ 1 ∥ 2 2 w T Σ 0 w + w T Σ 1 w = w T ( μ 0 − μ 1 ) ( μ 0 − μ 1 ) T w w T ( Σ 0 + Σ 1 ) w begin{aligned} J &=frac{left|boldsymbol{w}^{mathrm{T}} boldsymbol{mu}_{0}-boldsymbol{w}^{mathrm{T}} boldsymbol{mu}_{1}right|_{2}^{2}}{boldsymbol{w}^{mathrm{T}} boldsymbol{Sigma}_{0} boldsymbol{w}+boldsymbol{w}^{mathrm{T}} boldsymbol{Sigma}_{1} boldsymbol{w}} \ &=frac{boldsymbol{w}^{mathrm{T}}left(boldsymbol{mu}_{0}-boldsymbol{mu}_{1}right)left(boldsymbol{mu}_{0}-boldsymbol{mu}_{1}right)^{mathrm{T}} boldsymbol{w}}{boldsymbol{w}^{mathrm{T}}left(boldsymbol{Sigma}_{0}+boldsymbol{Sigma}_{1}right) boldsymbol{w}} end{aligned} J=wTΣ0w+wTΣ1w∥∥wTμ0−wTμ1∥∥22=wT(Σ0+Σ1)wwT(μ0−μ1)(μ0−μ1)Tw
观察上式目标函数，当我们的 J → M a x Jrightarrow Max J→Max则是我们想要的结果。式子太复杂，那我们再优化一下吧。

引入一下类内和间散度矩阵的知识：

类间散度矩阵用于表示各样本点围绕均值的散布情况。
类内散度矩阵用于表示样本点围绕均值的散步情况，关于特征选择和提取的结果，类内散布矩阵的积越小越好。
具体可参考这里，还有这里。

首先，我们来定义“类内散度矩阵”（within-class scatter matrix）
S w = Σ 0 + Σ 1 = ∑ x ∈ X 0 ( x − μ 0 ) ( x − μ 0 ) T + ∑ x ∈ X 1 ( x − μ 1 ) ( x − μ 1 ) T begin{aligned} mathbf{S}_{w} &=boldsymbol{Sigma}_{0}+boldsymbol{Sigma}_{1} \ &=sum_{boldsymbol{x} in X_{0}}left(boldsymbol{x}-boldsymbol{mu}_{0}right)left(boldsymbol{x}-boldsymbol{mu}_{0}right)^{mathrm{T}}+sum_{boldsymbol{x} in X_{1}}left(boldsymbol{x}-boldsymbol{mu}_{1}right)left(boldsymbol{x}-boldsymbol{mu}_{1}right)^{mathrm{T}} end{aligned} Sw=Σ0+Σ1=x∈X0∑(x−μ0)(x−μ0)T+x∈X1∑(x−μ1)(x−μ1)T
“类间散度矩阵”（between-class scatter matrix）：
S b = ( μ 0 − μ 1 ) ( μ 0 − μ 1 ) T mathbf{S}_{b}=left(boldsymbol{mu}_{0}-boldsymbol{mu}_{1}right)left(boldsymbol{mu}_{0}-boldsymbol{mu}_{1}right)^{mathrm{T}} Sb=(μ0−μ1)(μ0−μ1)T
然后我们的 J J J可以表示为
J = w T S b w w T S w w J=frac{boldsymbol{w}^{mathrm{T}} mathbf{S}_{b} boldsymbol{w}}{boldsymbol{w}^{mathrm{T}} mathbf{S}_{w} boldsymbol{w}} J=wTSwwwTSbw
这样看起来简单多了，这就是我们的LDA想要最大化的目标函数，比较专业的说法为， S b S_b Sb和 S w S_w Sw的“广义瑞利商”（generalizad Rayleigh quotient）。

关于“广义瑞利商”（generalizad Rayleigh quotient）的解释，可以参考这里和这里。
瑞利商经常出现在降维和聚类任务中，因为降维聚类任务往往能导出最大化最小化瑞利熵的式子，进而通过特征值分解的方式找到降维空间。
大体内容如下：

下面开始构建我们的函数及约束条件。

首先得确定我们的 w w w，由于 J J J的分母分子都是关于 w w w的二项式子，则与 w w w的长度无关，且只与方向有关。故我们令 w T S w w = 1 boldsymbol{w}^{mathrm{T}} mathbf{S}_{w} boldsymbol{w}=1 wTSww=1，则：
min ⁡ w − w T S b w s.t. w T S w w = 1 begin{array}{ll} min _{boldsymbol{w}} & -boldsymbol{w}^{mathrm{T}} mathbf{S}_{b} boldsymbol{w} \ text { s.t. } & boldsymbol{w}^{mathrm{T}} mathbf{S}_{w} boldsymbol{w}=1 end{array} minw s.t. −wTSbwwTSww=1
由拉格朗日乘数法（具体可参考CCA中Lagrange的应用）可得，
S b w = λ S w w (*) boldsymbol{S_bw=lambda S_w w}tag{*} Sbw=λSww(*)

其中， λ lambda λ为拉格朗日乘子。

由上“类间散度矩阵”可知， S b w boldsymbol{S_bw} Sbw为 μ 0 − μ 1 mu_0-mu_1 μ0−μ1的平方，故 S b w boldsymbol{S_bw} Sbw的方向则恒为 μ 0 − μ 1 mu_0-mu_1 μ0−μ1，向量的方向可以确定了，我们再令
S b w = λ ( μ 0 − μ 1 ) boldsymbol{S_bw=lambda (mu_0-mu_1)} Sbw=λ(μ0−μ1)

向量方向确定， λ lambda λ只是代表方向向量的长度，所以 S b w boldsymbol{S_bw} Sbw可由上式表达。可能会有人疑惑了，这里的 λ lambda λ和（*）式的 λ lambda λ是一个 λ lambda λ吗？
答案是肯定的。

将上式代入(*)式，可得关于 S w boldsymbol{S_w} Sw的式子：
w = S w − 1 ( μ 0 − μ 1 ) boldsymbol{w=S_w^{-1}(mu_0-mu_1)} w=Sw−1(μ0−μ1)
这里需要对 S w boldsymbol{S_w} Sw求逆，考虑到数值解的稳定性，常规实践操作中，需要对 S w boldsymbol{S_w} Sw进行奇异值分解（也就是我们在矩阵理论中学到的SVD方法），原理很简单，此处，即为 S w = U Σ V T boldsymbol{S_w=USigma V^{T}} Sw=UΣVT，其中 Σ Sigma Σ是一个实对角矩阵，对角线上的元素也就是所谓的“迹”是 S w boldsymbol{S_w} Sw的奇异值。我们需要求解的是 S w boldsymbol{S_w} Sw的逆，故式子变为了这样，
S w − 1 = V Σ − 1 U T boldsymbol{S_w^{-1}=VSigma^{-1} U^{T}} Sw−1=VΣ−1UT
至此，我们得到了 S w − 1 boldsymbol{S_w^{-1}} Sw−1，从而可求得直线向量 w w w，找到使得 J J J 最大的 w w w.

LDA还可从贝叶斯决策理论的角度来描述（关于贝叶斯可参考这里），可证明，当两类数据同先验、满足高斯分布（正态分布）且协方差相等时，LDA可以达到最优的分类效果。

上述讲了这么多都是二分类问题，那么关于多分类任务。

具体线性判别分析LDA和Fisher判别分析原理及推广可看 —> 这里

典型相关分析（CCA）

从字面意义上理解CCA，我们可以知道，简单说来就是对不同变量之间做相关分析。较为专业的说就是，一种度量两组变量之间相关程度的多元统计方法。

关于相似性度量距离问题，在这里有一篇Blog可以参考参考。

首先，从基本的入手。

当我们需要对两个变量 X ， Y X，Y X，Y进行相关关系分析时，则常常会用到相关系数来反映。学过概率统计的小伙伴应该都知道的吧。还是解释一下。

相关系数：是一种用以反映变量之间相关关系密切程度的统计指标。相关系数是按积差方法计算，同样以两变量与各自平均值的离差为基础，通过两个离差相乘来反映两变量之间相关程度；着重研究线性的单相关系数。
R ( X , Y ) = Cov ⁡ ( X , Y ) Var ⁡ [ X ] Var ⁡ [ Y ] R(X, Y)=frac{operatorname{Cov}(X, Y)}{sqrt{operatorname{Var}[X] operatorname{Var}[Y]}} R(X,Y)=Var[X]Var[Y] Cov(X,Y)
其中， C o v ( X , Y ) Cov(X,Y) Cov(X,Y)表示 X , Y X,Y X,Y的协方差矩阵， V a r [ X ] Var[X] Var[X]为 X X X的方差， V a r [ Y ] Var[Y] Var[Y]为 Y Y Y的方差.

复习了一下大学本科概率统计知识，那么，如果我们需要分析的对象是两组或者多组向量，又该怎么做呢？

CCA的数学表达：

这里举例两组变量 A ( a 1 , a 2 , . . . , a n ) , B ( b 1 , b 2 , . . . , b m ) A(a_1,a_2,...,a_n),B(b_1,b_2,...,b_m) A(a1,a2,...,an),B(b1,b2,...,bm)，那么我们的公式会是这样：
R ( X i , Y j ) = ∑ i = 1 , j = 1 n , m C o v ( X i , Y j ) V a r [ X i ] V a r [ Y j ] R(X_i,Y_j)=sum_{i=1,j=1}^{n,m} frac{Cov(X_i,Y_j)}{sqrt{Var[X_i]Var[Y_j]}} R(Xi,Yj)=i=1,j=1∑n,mVar[Xi]Var[Yj] Cov(Xi,Yj)
我们会得到一个这样的矩阵：
[ R ( X 1 , Y 1 ) . . . R ( X 1 , Y m − 1 ) R ( X 1 , Y m ) R ( X 2 , Y 1 ) . . . R ( X 2 , Y m − 1 ) R ( X 2 , Y m ) . . . . . . . . . . . . R ( X n , Y 1 ) . . . . . . R ( X n , Y m ) ] begin{bmatrix} R(X_1,Y_1) &... & R(X_1,Y_{m-1}) & R(X_1,Y_m)\R(X_2,Y_1) & ...& R(X_2,Y_{m-1})& R(X_2,Y_m)\ ...& ...& ...&... \ R(X_n,Y_1) & ...& ...&R(X_n,Y_m) end{bmatrix} ⎣⎢⎢⎡R(X1,Y1)R(X2,Y1)...R(Xn,Y1)............R(X1,Ym−1)R(X2,Ym−1)......R(X1,Ym)R(X2,Ym)...R(Xn,Ym)⎦⎥⎥⎤

这样的话，我们把每个变量的相关系数都求了出来，不知道会不会和我一样觉得这样很繁琐呢。如果我们能找到两组变量之间的各自的线性组合，那么我们就只分析讨论线性组合之间的相关分析。

典型相关系数：是先对原来各组变量进行主成分分析，得到新的线性关系的综合指标，再通过综合指标之间的线性相关系数来研究原各组变量间相关关系。

现在我们利用主成分分析（PCA）的思想，可以把多个变量与多个变量之间的相关转化成两个变量之间的相关。

先得到两组变量 ( A T , B T ) (A^T,B^T) (AT,BT)的协方差矩阵
Σ = [ Σ 11 Σ 12 Σ 21 Σ 22 ] Sigma=left[begin{array}{l} Sigma_{11} Sigma_{12} \ Sigma_{21} Sigma_{22} end{array}right] Σ=[Σ11 Σ12Σ21 Σ22]
其中， Σ 11 = C o v ( A ) , Σ 22 = C o v ( B ) , Σ 12 = Σ 12 T = C o v ( A , B ) Sigma_{11} = Cov(A),Sigma_{22} = Cov(B),Sigma_{12}=Sigma_{12}^T = Cov(A,B) Σ11=Cov(A),Σ22=Cov(B),Σ12=Σ12T=Cov(A,B).

把上面两组变量 A ( a 1 , a 2 , . . . , a n ) , B ( b 1 , b 2 , . . . , b m ) A(a_1,a_2,...,a_n),B(b_1,b_2,...,b_m) A(a1,a2,...,an),B(b1,b2,...,bm)分别组合成两个变量U、V，则用线性表示
U = t 1 a 1 + t 2 a 2 + . . . + t n a n , V = h 1 b 1 + h 2 b 2 + . . . + h m b m begin{matrix} U=t_1a_1+t_2a_2+...+t_na_n,\ \V=h_1b_1+h_2b_2+...+h_mb_m end{matrix} U=t1a1+t2a2+...+tnan,V=h1b1+h2b2+...+hmbm

然后，找出最大可能的相关系数 t k = ( t 1 , t 2 , . . . , t n ) T , h k = ( h 1 , h 2 , . . . , h m ) T {t_k}=(t_1,t_2,...,t_n)^T,{h_k}=(h_1,h_2,...,h_m)^T tk=(t1,t2,...,tn)T,hk=(h1,h2,...,hm)T,

使得， R ( U , V ) ⟶ M a x R(U,V)longrightarrow Max R(U,V)⟶Max，这样，就得到了典型相关系数；而其中的 U , V U,V U,V 为典型相关变量。

典型相关分析最朴素的思想：首先分别在每组变量中找出第一对典型变量，使其具有最大相关性，然后在每组变量中找出第二对典型变量，使其分别与本组内的第一对典型变量不相关，第二对本身具有次大的相关性。如此下去，直到进行到K步，两组变量的相关系被提取完为止，可以得到K组变量。

So，

注意：此时的 ( U , V ) (U,V) (U,V)若不能反映两组变量之间的相关关系，我们需要继续构造下一组关系变量来表示，具体可构造 K K K对这样的关系

直到 R ( U , V ) ⟶ M a x R(U,V)longrightarrow Max R(U,V)⟶Max为止
U k = t k T A = t 1 k a 1 + t 2 k a 2 + . . . + t n k a n V k = h k T B = h 1 k b 1 + h 2 k b 2 + . . . + h m k b m begin{matrix} U_k={t_k^T}{A}=t_{1k}a_1+t_{2k}a_2+...+t_{nk}a_n\ \ V_k={h_k^T}{B}=h_{1k}b_1+h_{2k}b_2+...+h_{mk}b_m end{matrix} Uk=tkTA=t1ka1+t2ka2+...+tnkanVk=hkTB=h1kb1+h2kb2+...+hmkbm

其中，我们需要一个约束条件满足，使得 R ( U , V ) ⟶ M a x R(U,V)longrightarrow Max R(U,V)⟶Max

V a r ( U k ) = V a r ( t k T A ) = t k T Σ 11 t k = 1 V a r ( V k ) = V a r ( h k T A ) = h k T Σ 22 h k = 1 C o v （ U k , U i ） = C o v ( U k , V i ) = C o v ( V i , U k ) = C o v ( V k , V i ) = 0 ( 1 < = i < k ) begin{matrix} Var(U_k)=Var({t_k^T}{A})={t_k^T}Sigma_{11}t_k=1\ \ Var(V_k)=Var({h_k^T}{A})={h_k^T}Sigma_{22}h_k=1\ \ Cov（U_k,U_i）=Cov(U_k,V_i)=Cov(V_i,U_k)=Cov(V_k,V_i)=0(1<=i 典型相关系数公式 R ( U , V ) R_{(U,V)} R(U,V)
R ( U , V ) = Cov ⁡ ( U , V ) Var ⁡ [ U ] Var ⁡ [ V ] = C o v ( U , V ) = t k T C o v ( A , B ) h k = t k T Σ 12 h k R_{(U,V)}=frac{operatorname{Cov}(U, V)}{sqrt{operatorname{Var}[U] operatorname{Var}[V]}}=Cov(U,V)={t_k}^TCov(A,B)h_k={t_k}^TSigma_{12} h_k R(U,V)=Var[U]Var[V] Cov(U,V)=Cov(U,V)=tkTCov(A,B)hk=tkTΣ12hk

在此约束条件下， t k , h k t_k,h_k tk,hk系数得到最大，则使得 R ( U , V ) R_{(U,V)} R(U,V)最大。

具体典型相关分析CCA算法及故障诊断应用可看 —> 这里

慢特征分析（SFA）

给定一个 i 维输入信号 x ( t ) = [ x 1 ( t ) , … , x I ( t ) ] T mathbf{x}(t)=left[x_{1}(t), ldots, x_{I}(t)right]^{T} x(t)=[x1(t),…,xI(t)]T，考虑一个输入-输出函数 g ( x ) = [ g 1 ( x ) , … , g J ( x ) ] T mathbf{g}(mathbf{x})=left[g_{1}(mathbf{x}), ldots, g_{J}(mathbf{x})right]^{T} g(x)=[g1(x),…,gJ(x)]T，每个分量都是一组K个非线性函数的加权和 h k ( x ) : g j ( x ) : = ∑ k = 1 K w j k h k ( x ) h_{k}(mathbf{x}): g_{j}(mathbf{x}):=sum_{k=1}^{K}w_{j k} h_{k}(mathbf{x}) hk(x):gj(x):=∑k=1Kwjkhk(x)。通常K > max(I, J)。应用 h = [ h 1 , … , h K ] T mathbf{h}=left[h_{1}, ldots, h_{K}right]^{T} h=[h1,…,hK]T对输入信号产生非线性扩展信号 z ( t ) : = h ( x ( t ) ) mathbf{z}(t):=mathbf{h}(mathbf{x}(t)) z(t):=h(x(t))。经过这种非线性展开后，可以将该问题在扩展的信号分量 z k ( t ) z_k(t) zk(t)中视为线性问题。这是将非线性问题转化为线性问题的常用方法。一个众所周知的例子是支持向量机(Vapnik, 1995)。权向量 w j = [ w j 1 , … , w j K ] T mathbf{w}_{j}=left[w_{j 1}, ldots, w_{j K}right]^{T} wj=[wj1,…,wjK]T进行学习，第j个输出信号分量由 y j ( t ) = g j ( x ( t ) ) = w j T h ( x ( t ) ) = w j T z ( t ) y_{j}(t)=g_{j}(mathbf{x}(t))=mathbf{w}_{j}^{T} mathbf{h}(mathbf{x}(t))=mathbf{w}_{j}^{T} mathbf{z}(t) yj(t)=gj(x(t))=wjTh(x(t))=wjTz(t)给出。

目标(见方程1)是优化输入-输出函数，从而使权值达到
Δ ( y j ) = ⟨ y ˙ j 2 ⟩ = w j T ⟨ z ˙ z ˙ T ⟩ w j Deltaleft(y_{j}right)=leftlangledot{y}_{j}^{2}rightrangle=mathbf{w}_{j}^{T}leftlangledot{mathbf{z}} dot{mathbf{z}}^{T}rightrangle mathbf{w}_{j} Δ(yj)=⟨y˙j2⟩=wjT⟨z˙z˙T⟩wj (3.1)

是最小的。

假设所选择的非线性函数hkare使扩展信号z(t)具有零均值和单位协方差矩阵。这样一组非线性函数 h k h_k hk可以很容易地通过一个球阶从任意集合 h k ′ h_{k}^{prime} hk′中导出，如下所述。然后我们发现约束条件(见公式2-4)

⟨ y j ⟩ = w j T ⟨ z ⟩ ⏟ = 0 = 0 , （ 3.2 ） ⟨ y j 2 ⟩ = w j T ⟨ z z T ⟩ ⏟ = I w j = w j T w j = 1 （ 3.3 ） , ∀ j ′ < j : ⟨ y j ′ y j ⟩ = w j ′ T ⟨ z z T ⟩ ⏟ = I w j = w j ′ T w j = 0 , （ 3.4 ） begin{array}{r} leftlangle y_{j}rightrangle = mathbf{w}_{j}^{T} underbrace{langlemathbf{z}rangle}_{ = 0} = 0, （3.2）\ leftlangle y_{j}^{2}rightrangle = mathbf{w}_{j}^{T} underbrace{leftlanglemathbf{z z}^{T}rightrangle}_{ = mathbf{I}} mathbf{w}_{j} = mathbf{w}_{j}^{T} mathbf{w}_{j} = 1 （3.3）, \ forall j^{prime}

当且仅当我们约束权重向量为向量的标准正交集时，自动满足。

因此，对于输入输出函数的第一个分量，优化问题简化为寻找使方程(3.1)中的 Δ ( y 1 ) Deltaleft(y_{1}right) Δ(y1)最小的赋范权向量。解是矩阵 ⟨ z ˙ z ˙ T ⟩ leftlangledot{mathbf{z}} dot{mathbf{z}}^{T}rightrangle ⟨z˙z˙T⟩的赋范特征向量，它对应于最小的特征值(cf. Mitchison, 1991)。下一个更高的特征值的特征向量产生输入-输出函数的下一个分量与下一个更高的1个值。这就引出了解决上述优化问题的算法。

明确区分原始信号、来自训练数据的精确归一化信号和来自测试数据的近似归一化信号是有用的。设 x ~ ( t ) tilde{mathbf{x}}(t) x~(t)是一个可以有任意均值和方差的原始输入信号。为了计算方便和显示的目的，信号被归一化为零的平均值和单位方差。这种归一化对训练数据x(t)是精确的。用相同的偏移量和因子来校正测试数据，通常会得到一个近似归一化的输入信号 x ′ ( t ) mathbf{x}^{prime}(t) x′(t)，因为每个数据样本的均值和方差略有不同，而归一化总是使用从训练数据中确定的偏移量和因子来完成。在下面，原始信号有一个波浪线，测试数据有一个破折号;没有波浪线或破折号的符号通常(但不总是)指的是规范化训练数据。

该算法现在有如下形式(参见下图):

输入信号。用于训练，i维输入信号为 x ~ ( t ) tilde{mathbf{x}}(t) x~(t)。
输入信号归一化。对输入信号进行归一化得到

x ( t ) : = [ x 1 ( t ) , … , x I ( t ) ] T （ 3.5 ） with x i ( t ) : = x ~ i ( t ) − ⟨ x ~ i ⟩ ⟨ ( x ~ i − ⟨ x ~ i ⟩ ) 2 ⟩ , （ 3.6 ） so that ⟨ x i ⟩ = 0 （ 3.7 ） and ⟨ x i 2 ⟩ = 1 （ 3.8 ） begin{aligned} mathbf{x}(t) &:=left[x_{1}(t), ldots, x_{I}(t)right]^{T} （3.5）\ text { with } quad x_{i}(t) &:=frac{tilde{x}_{i}(t)-leftlangletilde{x}_{i}rightrangle}{sqrt{leftlangleleft(tilde{x}_{i}-leftlangletilde{x}_{i}rightrangleright)^{2}rightrangle}}, （3.6） \ text { so that } quadleftlangle x_{i}rightrangle &=0 （3.7）\ text { and } quadleftlangle x_{i}^{2}rightrangle &=1 （3.8） end{aligned} x(t) with xi(t) so that ⟨xi⟩ and ⟨xi2⟩:=[x1(t),…,xI(t)]T（3.5）:=⟨(x~i−⟨x~i⟩)2⟩ x~i(t)−⟨x~i⟩,（3.6）=0（3.7）=1（3.8）

非线性扩张。应用一组非线性函数 h ~ ( x ) tilde{mathbf{h}}(mathbf{x}) h~(x)，生成扩展信号 z ~ ( t ) tilde{mathbf{z}}(t) z~(t)。这里使用了一次的所有单项式(导致线性 S F A SFA SFA有时用 S F A 1 SFA_1 SFA1表示)或一次和两次的单项式，包括混合项，如 x 1 x 2 x_1x_2 x1x2(导致二次SFA有时用 S F A 2 SFA_2 SFA2表示)，但也可以使用任何其他的函数集。因此对于二次 S F A SFA SFA，

h ~ ( x ) : = [ x 1 , … , x I , x 1 x 1 , x 1 x 2 , … , x I x I ] T （ 3.9 ） z ~ ( t ) : = h ~ ( x ( t ) ) = [ x 1 ( t ) , … , x I ( t ) , x 1 ( t ) x 1 ( t ) , x 1 ( t ) x 2 ( t ) , … , x I ( t ) x I ( t ) ] T （ 3.10 ） begin{aligned} tilde{mathbf{h}}(mathbf{x}):=&left[x_{1}, ldots, x_{I}, x_{1} x_{1}, x_{1} x_{2}, ldots, x_{I} x_{I}right]^{T} （3.9） \ tilde{mathbf{z}}(t):=tilde{mathbf{h}}(mathbf{x}(t))=&left[x_{1}(t), ldots, x_{I}(t), x_{1}(t) x_{1}(t),right.\ &left.x_{1}(t) x_{2}(t), ldots, x_{I}(t) x_{I}(t)right]^{T} （3.10） end{aligned} h~(x):=z~(t):=h~(x(t))=[x1,…,xI,x1x1,x1x2,…,xIxI]T（3.9）[x1(t),…,xI(t),x1(t)x1(t),x1(t)x2(t),…,xI(t)xI(t)]T（3.10）

属于 h ~ ( x ) tilde{mathbf{h}}(mathbf{x}) h~(x)和属于 z ~ ( t ) tilde{mathbf{z}}(mathbf{t}) z~(t)的一、二度分量词的维数为 K = I + I ( I + 1 ) / 2 K = I + I(I + 1)/2 K=I+I(I+1)/2。

球体。将扩展的信号属于 z ~ ( t ) tilde{mathbf{z}}(mathbf{t}) z~(t)通过仿射变换归一化，以生成具有零均值和单位协方差矩阵 I I I的 z ( t ) z(t) z(t)，

z ( t ) : = S ( z ~ ( t ) − ⟨ z ~ ⟩ ) , （ 3.11 ） with ⟨ z ⟩ = 0 （ 3.12 ） and ⟨ z z T ⟩ = I （ 3.13 ） begin{aligned} mathbf{z}(t) &:=mathbf{S}(tilde{mathbf{z}}(t)-langletilde{mathbf{z}}rangle), （3.11） \ text { with } quadlanglemathbf{z}rangle &=mathbf{0} （3.12） \ text { and }leftlanglemathbf{z} mathbf{z}^{T}rightrangle &=mathbf{I} （3.13） end{aligned} z(t) with ⟨z⟩ and ⟨zzT⟩:=S(z~(t)−⟨z~⟩),（3.11）=0（3.12）=I（3.13）

这种归一化称为球面化(或白化)。矩阵 S S S是球面矩阵，可以通过对矩阵的主成分分析(PCA)确定 ( Z ~ ( t ) − ⟨ z ~ ⟩ ) (tilde{mathbf{Z}}(t)-langletilde{mathbf{z}}rangle) (Z~(t)−⟨z~⟩)。因此，它依赖于特定的训练数据集。这也定义了

h ( x ) : = S ( h ~ ( x ) − ⟨ z ~ ⟩ ) mathbf{h}(mathbf{x}):=mathbf{S}(tilde{mathbf{h}}(mathbf{x})-langletilde{mathbf{z}}rangle) h(x):=S(h~(x)−⟨z~⟩) （3.14）

它是一个标准化函数，而 z ( t ) z(t) z(t)是球化数据。

主成分分析。将主成分分析法应用于矩阵 ⟨ z ˙ z ˙ T ⟩ leftlangledot{mathbf{z}} dot{mathbf{z}}^{T}rightrangle ⟨z˙z˙T⟩。具有最小特征值的 J J J个特征向量 λ j λ_j λj得到归一化权向量

w j : ⟨ z ˙ z ˙ T ⟩ w j = λ j w j （ 3.15 ） with λ 1 ≤ λ 2 ≤ ⋯ ≤ λ J （ 3.16 ） begin{array}{ll} mathbf{w}_{j}: & leftlangledot{mathbf{z}} dot{mathbf{z}}^{T}rightrangle mathbf{w}_{j}=lambda_{j} mathbf{w}_{j} （3.15） \ text { with } & lambda_{1} leq lambda_{2} leq cdots leq lambda_{J} （3.16） end{array} wj: with ⟨z˙z˙T⟩wj=λjwj（3.15）λ1≤λ2≤⋯≤λJ（3.16）

提供输入-输出功能

g ( x ) : = [ g 1 ( x ) , … , g J ( x ) ] T （ 3.17 ） with g j ( x ) : = w j T h ( x ) （ 3.18 ） begin{aligned} mathbf{g}(mathbf{x}) &:=left[g_{1}(mathbf{x}), ldots, g_{J}(mathbf{x})right]^{T} （3.17）\ text { with } quad g_{j}(mathbf{x}) &:=mathbf{w}_{j}^{T} mathbf{h}(mathbf{x}) （3.18） end{aligned} g(x) with gj(x):=[g1(x),…,gJ(x)]T（3.17）:=wjTh(x)（3.18）

输出信号

y ( t ) : = g ( x ( t ) ) （ 3.19 ） with ⟨ y ⟩ = 0 （ 3.20 ） ⟨ y y T ⟩ = I （ 3.21 ） and Δ ( y j ) = ⟨ y ˙ j 2 ⟩ = λ j . （ 3.22 ） begin{aligned} mathbf{y}(t) &:=mathbf{g}(mathbf{x}(t)) （3.19）\ text { with } &langlemathbf{y}rangle &=mathbf{0} （3.20）\ leftlanglemathbf{y} mathbf{y}^{T}rightrangle &=mathbf{I} （3.21）\ text { and } Deltaleft(y_{j}right)=leftlangledot{y}_{j}^{2}rightrangle &=lambda_{j} . （3.22） end{aligned} y(t) with ⟨yyT⟩ and Δ(yj)=⟨y˙j2⟩:=g(x(t))（3.19）⟨y⟩=I（3.21）=λj.（3.22）=0（3.20）

输出信号的各分量的均值和单位方差都为零，并且是不相关的。

重复。如果需要，将输出信号 y ( t ) y(t) y(t)(或它的前几个分量或不同输出信号的组合)作为输入信号 x ( t ) x(t) x(t)，用于学习算法的下一个应用。继续步骤3。
测试。为了对测试信号进行系统测试，将步骤2至6中推导的归一化和输入-输出函数应用于新的输入信号(属于x0(T))。注意，这个测试信号需要使用与训练信号相同的偏移量和因子进行归一化，以准确地再现学习的输入-输出关系。因此，训练信号只是近似地归一化降低

x ′ ( t ) : = [ x 1 ′ ( t ) , … , x I ′ ( t ) ] T （ 3.23 ） with x i ′ ( t ) : = x ~ i ′ ( t ) − ⟨ x ~ i ⟩ ⟨ ( x ~ i − ⟨ x ~ i ⟩ ) 2 ⟩ （ 3.24 ） so that ⟨ x i ′ ⟩ ≈ 0 （ 3.25 ） and ⟨ x i ′ 2 ⟩ ≈ 1 （ 3.26 ） begin{aligned} & mathbf{x}^{prime}(t):=left[x_{1}^{prime}(t), ldots, x_{I}^{prime}(t)right]^{T} （3.23）\ text { with } quad x_{i}^{prime}(t) &:=frac{tilde{x}_{i}^{prime}(t)-leftlangletilde{x}_{i}rightrangle}{sqrt{leftlangleleft(tilde{x}_{i}-leftlangletilde{x}_{i}rightrangleright)^{2}rightrangle}}（3.24） \ text { so that } quadleftlangle x_{i}^{prime}rightrangle & approx 0 （3.25）\ text { and }leftlangle x_{i}^{prime 2}rightrangle & approx 1 （3.26） end{aligned} with xi′(t) so that ⟨xi′⟩ and ⟨xi′2⟩x′(t):=[x1′(t),…,xI′(t)]T（3.23）:=⟨(x~i−⟨x~i⟩)2⟩ x~i′(t)−⟨x~i⟩（3.24）≈0（3.25）≈1（3.26）

只有在测试信号对训练信号具有代表性的情况下，归一化才是准确的。输出信号也是如此

y ′ ( t ) : = g ( x ′ ( t ) ) （ 3.27 ） with ⟨ y ′ ⟩ ≈ 0 , （ 3.28 ） and ⟨ y ′ y ′ T ⟩ ≈ I . （ 3.29 ） begin{aligned} mathbf{y}^{prime}(t) &:=mathbf{g}left(mathbf{x}^{prime}(t)right)（3.27） \ text { with } quadleftlanglemathbf{y}^{prime}rightrangle & approx mathbf{0}, （3.28）\ text { and }leftlanglemathbf{y}^{prime} mathbf{y}^{prime T}rightrangle & approx mathbf{I} . （3.29） end{aligned} y′(t) with ⟨y′⟩ and ⟨y′y′T⟩:=g(x′(t))（3.27）≈0,（3.28）≈I.（3.29）

由于实际原因，在步骤4和步骤5中使用了奇异值分解而不是PCA。对于一些特征值非常接近于零的退化数据，奇异值分解是一种较好的分析方法，然后在步骤4中丢弃这些特征值。非线性展开有时会导致退化数据，因为它产生高度冗余的表示，其中一些成分可能具有线性关系。一般来说，特征值接近于零的信号分量通常包含噪声，如舍入误差，在归一化后会非常迅速地波动和在任何情况下，SFA都不会在步骤5中选择。因此，决定哪些小组件应该被丢弃并不重要。

【图3】执行SFA的两种可能的网络结构。(上)用一组单位来解释，在树(粗线)上进行复杂的计算，如sigma-pi单位。(下)解释为简单单元与隐含层中固定的非线性单元的分层网络，如径向基函数网络具有不可适应的隐含单元。在这两种情况下，输入-输出函数分量由 g j ( x ) = w j T h ( x ) = w ~ j 0 + w ~ j T h ~ ( x ) g_{j}(mathbf{x})=mathbf{w}_{j}^{T} mathbf{h}(mathbf{x})=tilde{w}_{j 0}+tilde{mathbf{w}}_{j}^{T} tilde{mathbf{h}}(mathbf{x}) gj(x)=wjTh(x)=w~j0+w~jTh~(x)给出，并带有适当的原始 w ~ ( J ) tilde{mathbf{w}}(mathbf{J}) w~(J)权重向量。这里假设输入信号分量是标准化的。

不是直接用 Δ Delta Δ的值，而是用一种更直观的解释来测量信号的不变性是有用的。

可以是定义为的指标η

g j ( x ) = w j T h ( x ) = w ~ j 0 + w ~ j T h ~ ( x ) η ( y ) : = T 2 π Δ ( y ) g_{j}(mathbf{x})=mathbf{w}_{j}^{T} mathbf{h}(mathbf{x})=tilde{w}_{j 0}+tilde{mathbf{w}}_{j}^{T} tilde{mathbf{h}}(mathbf{x})eta(y):=frac{T}{2 pi} sqrt{Delta(y)} gj(x)=wjTh(x)=w~j0+w~jTh~(x)η(y):=2πTΔ(y) （3.30）

对 t ∈ [ t 0 , t 0 + T ] t ∈ [t_0,t_0+T] t∈[t0,t0+T], 对于一个纯正弦波 y ( t ) : = 2 sin ⁡ ( n 2 π t / T ) y(t):=sqrt{2} sin (n 2 pi t / T) y(t):=2 sin(n2πt/T)整数数目的振荡，η(y)就是振荡的数目，即η(y) = n。因此，任意信号的指数η表示相同值的纯正弦波的振荡次数，至少对η的整数值是这样的。低η值表示信号慢。由于从测试数据导出的输出信号只是近似归一化， η ( y 0 ) η(y_0) η(y0)意味着包含 y 0 y_0 y0到零的平均值和单位方差的精确归一化，以使 η eta η 指标独立于一个偶然的比例因子。

具体慢特征分析SFA原理及应用可看 —> 这里

非线性数据降维

实际工业过程往往表现出复杂的特性，过程的非线性特性尤为突出，因此非线性数据降维方法在故障诊断领域得到了广泛的关注和应用。针对工业过程中存在的大量非线性过程，尤其是当输入的取值范围较大时，很多过程信息及这种非线性关系无法再被PCA描述。针对非线性恃征的故障诊断研究从九十年代发展至今，主要提出了以下几类方法：神经网络非线性方法；（Kernel）方法；流形学习（Maniford Learning）方法。

神经网络非线性方法

神经网络非线性PCA方法最早是有Kramer提出的，其模型是一个五层神经网络，输入和输出采用同一组数据。由于中间层神经元个数小于输入层的个数，实现了降维的效果。有作者指出该网络结构在线性情况下并不能保证与线性PCA等价，提出了一种改进的方法，即先计算出每个训练数据点的具有物理意义的主曲线的低维表示后，再使用两个三层神经网络分别实现降维和重反构映射；有学者引入核概率密度函数研究了非线性PCA中主元不服从正态分布的问题；有学者研究了非线性过程数据含有的动态性问题。

其中神经网络还可结合LDA、ICA、SFA、CCA等上述线性方法，这里就不再展开啦！~~

这里只以神经网络加入PCA方法为例。但是针对神经网络（NN）加入PCA的方法，从实际上来讲目前任然存在很多的问题，首当其冲的就是整个神经网络的训练复杂度，会很高，提高成本，在实际工业生产中的数据并不是完整的，会存在样本的不均衡和故障样本的缺失或样本量极少，所以，整个模型的泛化能力不是很好或者说是难以保证。怎样确定整个神经网络模型的结构及对应生产工况数据的规模，优化算法，效率及在线功能实现，保证学习样本的完整性和价值问题仍未得到解决，故对于此类加入神经网络的降维方法，不是很看好，简单来说就是——太玄了！

核（Kernel）方法

核方法是将数据的原始输入空间非线性地变换到一个高维的核特征空间后、再使用线性技术的一种间接处理方法。

最开始将线性推广至非线性的为支持向量机（SVM），后来又有了Kernel PCA用于人脸识别和语音识别，同时在过程监控领域也有了很大的研究进展，推动了过程监控故障诊断的进步，虽然有很多学者将核kernel应用到了非线性降维方法中来，但均存在着一个很重要的问题，就是核函数究竟该怎么选取。而且，对于加入核函数的方法，仅仅依赖于反映数据中的非线性，在其降维用的目标函数中没有考虑到保持其数据结构的完整性，甚至存在将数据强行变为非线性结构的缺陷，由此，此种方法并不是很好，还需要进一步的研究和改善。

线性核函数

让转换函数 ϕ ( x ) = x phi (x)=x ϕ(x)=x，则得到线性核函数，则两个向量的点积为：
k ( x , x ′ ) = x T x ′ k(x,x^{'})=x^{T}x^{'} k(x,x′)=xTx′
线性核函数的特征空间F的维度与输入空间 χ chi χ的维度一样，每个向量的特征数也一样，其中 x = ( x 1 , x 2 , . . . , x n ) x=(x_{1},x_{2},...,x_{n}) x=(x1,x2,...,xn)叫做特征， x x x代表一个样本。

当不需要在特征空间进行运算时，可以使用线性核函数。

如果原始数据已经是高维且可进行比较，在输入空间中线性可分。

线性内核适用于由大量固定长度的特征数据表示的对象。

高斯核函数

高斯核函数其又称作——Square exponential kernel（SE Kernel ）or radial basis function(RBF，径向基函数)，其形如：
k ( x , x ′ ) = e − 1 2 ( x − x ′ ) T Σ − 1 ( x − x ′ ) k(x,x^{'})=e^{-frac{1}{2}(x-x^{'})^{T}Sigma^{-1} (x-x^{'})} k(x,x′)=e−21(x−x′)TΣ−1(x−x′)
其中， Σ Sigma Σ表示的是输入样本的每个特征数据所对应的协方差，P维数据矩阵。

当 Σ Sigma Σ是对角矩阵时，为
k ( x , x ′ ) = e − 1 2 ∑ j = 1 p ( x j − x j ′ ) 2 σ j 2 k(x,x^{'})=e^{-frac{1}{2}sum_{j=1}^{p} frac{(x_{j}-x^{'}_{j})^2}{sigma^{2}_{j}}} k(x,x′)=e−21∑j=1pσj2(xj−xj′)2
其中， σ j sigma_j σj被定义为特征 j j j 的伸缩尺度。

若 Σ Sigma Σ是球形的时候，则为
k ( x , x ′ ) = e − ∥ x − x ′ ∥ 2 2 σ 2 , σ j = σ , ∀ j k(x,x^{'})=e^{-frac{left | x-x^{'} right |^{2}}{2sigma^{2}}},sigma_j=sigma,forall j k(x,x′)=e−2σ2∥∥x−x′∥∥2,σj=σ,∀j
该核函数的特征空间的维度是无限的。核函数避免转换函数的计算，故可用相对马氏距离来机选n阶的Gram方阵，即使已经将隐式的将对象投影到无限维度的特征空间中去，依旧有效。

流行学习方法

流形学习是非线性降维的重要研究领域，也是近些年来，利用流行学习方法进行特征提取和故障诊断的一种新的探索，并受到了广泛的关注。

流形(manifold)是一般几何对象的总称，包括各种维度的曲线与曲面等，和一般的降维分析一样，流形学习是把一组在高维空间中的数据在低维空间中重新表示。不同之处是，在流形学习中假设：所处理的数据采样与一个潜在的流形上，或者说对于这组数据存在一个潜在的流形。

流形上的点本身是没有坐标的，所以为了表示这些数据点，我们把流形放入到外围空间(ambient space)，用外围空间上的坐标来表示流形上的点，例如三维空间 R 3 R^3 R3 中球面是一个2维曲面，即球面上只有两个自由度，但我们一般采用外围空间R^3 空间中的坐标来表示这个球面。

流形学习 = 微分流形 + 黎曼几何

流行学习可以概括为：在保持流形上点的某些几何性质特征的情况下，找出一组对应的内蕴坐标(intrinsic coordinate)，将流形尽量好的展开在低维平面上，这种低维表示也叫内蕴特征(intrinsic feature)，外围空间的维数叫观察维数，其表示叫自然坐标，在统计上称为observation。(引用自此Blog)

LLE(Locally Linear Embedding) 局部线性嵌入

LLE算法是在2000年被S.T.Roweis和L.K.Saul提出来的非线性降维方法，它的基本思想是将全局非线性转化为局部线性，而互相重叠的局部邻域能够提供全局结构的信息，这样对每个局部进行线性降维后，在按照某种规则将结果组合在一起，就能够得到低维的全局坐标表示。

LLE算法有一个前提假设：采样数据所在的低维流形在局部是线性的，即每个采样点可以用它的近邻点线性表示。这个假设可以由最终结果的有效性上得到证明。该算法的学习目标是：在低维空间中保持每个邻域中的权值不变，即假设嵌入映射在局部是线性的条件下，最小化重构误差。

上图为LLE算法的运行步骤

设有高维空间 R D R^D RD中的N个公维随机向量 { x i } = ( x 1 , x 2 , . . . , x N ) , i ∈ [ 1 ， N ] {x_{i}}=(x_{1},x_{2},...,x_{N}),iin [1，N] {xi}=(x1,x2,...,xN),i∈[1，N]，要通过降维把它们映射到低维空间 R d R^{d} Rd中，需要求出映射后的N个ｄ维向量台 { y i } , i ∈ [ 1 ， N ] {y_{i}},iin [1，N] {yi},i∈[1，N]。我们当然希望能够用显示的映射 F F F来表示，但{ x i x_i xi}柄是非线性结构时，求出显示映射 F F F相当困难，基于流形学习的原理，可Ｗ将非线性结构的局部用线性结构来近似，即可Ｗ将非线性结构的降维分解为局部的线性降维。这就是LLE算法的基本原理。

主要步骤：

寻找每个样本点的 k k k个近邻点；
由每个样本点的近邻点计算出该样本点的局部重建权值矩阵
由该样本点的局部重建权值矩阵和其近邻点计算出该样本点的输出向量

具体的详细计算步骤就不放这里了，有一个大佬的写的不错，这跳转这里

LE (Laplacian Eigenmaps) 拉普拉斯特征映射

拉普拉斯特征映射(Laplacian Eigenmaps，LE)也是一种局部的流行学习方法，和LE的方法有点相似，LE算法也具有直观降维的目标函数，也就是说，通过保持了高维数据之间的近邻关系来得到低维度数据流行结构。

LE算法的主要有以下四步：

建立无向图：对于每个高维样本数据点 { x i } = ( x 1 , x 2 , . . . , x N ) , i ∈ [ 1 ， N ] {x_{i}}=(x_{1},x_{2},...,x_{N}),iin [1，N] {xi}=(x1,x2,...,xN),i∈[1，N]，寻求其近邻数据点（去欧氏距离或最小K近邻的数据点），由此建立邻域关系图（近邻图）
计算每条近邻数据点执念的权值，（不相连的边权值为0），构建权值矩阵 W = { w i j } W={w_{ij} } W={wij}.
其中数据点中权值的计算方法有两种：
（如果xi和xj在邻域关系图中为近邻点，则将它们之间的权值赋值为 w i j = 1 w_{ij}=1 wij=1；相反，则为0.简单说来就是，节点相邻取1，不相邻取0。）
（热核权值：如果xi和xj在邻域中为近邻点，其热核函数为其赋值为： w i , j = e x p ( ∥ x i − x j ∥ 2 σ 2 ) w_{i,j}=exp(frac{left | x_{i}-x_{j} right |^{2}}{sigma^2}) wi,j=exp(σ2∥xi−xj∥2)）,将其 w i , j w_{i,j} wi,j带入到拉普拉斯矩阵中得到拉式矩阵L；相反，则权值赋值为0.）
求解特征向量方差（映射结果）：目标函数—> m i n Σ i j ( y i − y j ) 2 w i j min Sigma_{ij}(y_i-y_j)^2w_{ij} minΣij(yi−yj)2wij ,经过推导，对低维数据的计算转化为特征值分解的方式，也即使球 L y = λ D y Ly=lambda Dy Ly=λDy,将点 x i x_i xi映射到 y i y_i yi中。
使用最小的n和非零特征值对应的特征向量作为降维后的结果输出。

LPP(Locality Preserving Projection) 局部保留投影

局部保留投影LPP算法提出的目的是为了实现非线性流形的学习和分析，LPP可以提取最具有判别性的特征来进行降维，是一种保留了局部信息，降低影响图像识别的诸多因素的降维方法，这种算法本质上是一种线性降维方法，由于其巧妙的结合了拉普拉斯特征映射算法（LE）的思想，从而可以在对高维数据进行降维后有效地保留数据内部的非线性结构。

与其他非线性降维方法相比，LPP方法可以将新增的测试数据点，通过映射在降维后的子空间找到对应的位置，而其他非线性方法只能定义训练数据点，无法评估新的测试数据。LPP方法可以很容易地将新的测试数据点根据特征映射关系（矩阵），投影映射在低维空间中。

将 n 维原数据映射为 l 维数据，l<< n ；实现数据降维，样本个数为m 。（引用此大佬）

转自大佬

或许看到这已经很累了，确实，写这么多也很累，后面都直接转载了我（菜鸡），但是干什么又不累呢，是吧，一会儿还得刷力扣~

加油！~

对了，整理的代码放在这里了—> 数据降维算法代码集合

链接：https://pan.baidu.com/s/1yusgn7VtgSADEu7tiBF3nA （可直接点上面链接跳转的哟~）
提取码：yyds

❤坚持读Paper，坚持做笔记，坚持学习，坚持刷力扣LeetCode❤！！！
坚持刷题！！！打天梯！！！
⚡To Be No.1

⚡⚡哈哈哈哈

⚡创作不易⚡，过路能❤关注、收藏、点个赞❤三连就最好不过了

ღ( ´･ᴗ･` )

❤

『
物质决定意识，意识反作用于物质。
』

线性与非线性数据降维方法汇总（Python代码实现）

Python相关栏目本月热门文章