偏最小二乘法（NIPALS经典实现--未简化）

billy145533

已于 2024-04-14 12:27:17 修改

阅读量1.1w

点赞数 18

分类专栏：偏最小二乘法文章标签：最小二乘法算法机器学习

于 2019-05-20 21:00:07 首次发布

本文链接：https://blog.csdn.net/billy145533/article/details/90342460

版权

偏最小二乘法专栏收录该内容

13 篇文章 30 订阅

订阅专栏

前言

偏最小二乘法(Partial Least Squares，PLS)的实现方式很多，在理解优化版的PLS代码之前，有必要理解最原始的PLS思想。网上关于PLS的资料已是浩如烟海，但看完总觉得在理解上差一些意思，没有比较好的代码可以一步一步跟踪。本文的目的是力求实现一遍最原始的PLS算法NIPALS，推导其中的过程并给出MATLAB实现。看PLS代码时，总是头痛于代码表述和原理有些出入，我希望自己的代码能将两者统一起来，方便大家理解。

PLS虽小，却是小算法里的集大成者。在这里，我们可以看到PCA，SVD，CCA，LS等算法的影子。“今晚是首次独家大揭露。为了让大家能够看得清楚，我们一连五晚不断地重复、重复再重复整个的解剖过程，务求让大家能够重复、重复再重复地看个仔细。”

偏最小二乘法原理

偏最小二乘法简单的说是潜在空间下的最小二乘法。由于高维数据常常导致其协方差矩阵为奇异，无法使用最小二乘进行回归，实际中常常利用降维的方法提取数据中的重要部分进行建模。这种数据提取的方式构成了PLS与其他基于潜在空间方法的差别。下图展示PLS的基本原理，T和U分别代表从原始数据X和Y中抽取用于建模的部分，即构成了PLS的潜在空间。

对于给定Y，X，采用偏最小二乘法建立其Y，X之间的关系，其实现过程如下

1. 将X，Y 标准化，分别得到E，F

2. 计算投影轴w，v,使得两个成分的协方差最大，即Cov(Ew,Fv)最大

$\max \boldsymbol{Cov}(Ew,Fv)$

由于E和F已经标准化，其协方差为

$\boldsymbol{Cov}(Ew,Fv) = v^TF^TEw \\ \mathbf{s.t.} \ ||v||=1,||w||=1$

利用拉格朗日乘子法

$\small L =v^TF^TEw-\frac{\lambda}{2}(w^Tw-1)-\frac{\theta}{2}(v^Tv-1)$

对v，w求偏导数

$\small \\ \frac{\partial L}{\partial w} =E^TFv-\lambda w = 0 \\ \frac{\partial L}{\partial v} =F^TEw-\theta v = 0$

上面两式分别乘以 $w^T$ ， $v^T$ 得到

$\small \\w^TE^TFv-w^T\lambda w =0\Rightarrow w^TE^TFv = \lambda \\v^TF^TEw-v^T\theta v = 0\Rightarrow v^TF^TEw = \theta \\(v^TF^TEw)^T = w^TE^TFv\Rightarrow \theta = \lambda$

将 $\small v =F^TEw /\lambda$ 带入

$\small E^TFv-\lambda w = E^TFF^TEw/\lambda- \lambda w\Rightarrow E^TFF^TEw=\lambda ^2 w$

同理可以得到

$(F^TEE^TF)v=\lambda ^2 v$

在奇异分解SVD和PCA那两章我们已经讨论过这种情况，w 是矩阵 $\small E^TFF^TE$ 最大特征值对应的向量， $\small v =F^TEw /\lambda$ ，因此只需要计算w就足够了。

3. 计算个投影得分向量t = E*w ,u =F*v

$\small \begin{array}{lcr} E =tp^t+E^*\\ F =uq^t+F^*\\ F =tr^t+F^{**} \\ p = \frac{E^Tt}{\left \| t \right \|^2}\\ q= \frac{F^Tt}{\left \| u \right \|^2}\\ r = \frac{F^Tt}{\left \| t \right \|^2}\\ \\ \end{array}$

有些郁闷，t居然不是由w展开，作为E中提取的第一个矩阵形式的成分。显然w并非E的特征向量，它并不能保证由它展开矩阵提取到最多的成分。但是，w和p确实存在一定的关系

$w^Tp = w^T\frac{E^Tt}{\left \| t \right \|^2} = 1$

再观察一下p的意义

$p = \frac{E^Tt}{\left \| t \right \|^2} = \frac{E^TEw}{\left \| t \right \|^2}$

w是F'E最大特征值对应的特征向量，将自相关协方差矩阵E'E投到w上，确实不好理解。姑且就看看p是怎么求得的吧。p将t粘到E列空间中，希望两者差距越小越好，根据最小二乘法，E将各列分别投影到 t的单位向量上，使得误差最小，各列的投影长度构成了p。w和p的具体差别可以参考另一篇博文。

4. 计算残差

$\small \\ E^* =E-tp^t\\ F^* =F-tr^t$

用E*，F*代替 E，F，跳到第2步，继续求t，p，r。一般而言，迭代次数不多于X的秩。PLS有一个根据残差自动停止的算法，这里为了简化程序，暂时不研究，后续再展开。

5. 建模

$\small \\ E = t_1p_1^T+t_2p_2^T+...+t_rp_r^T=T*P^T \\ F = t_1r_1^T+t_2r_2^T+...+t_rr_r^T+F^* = T*R^T \\w_1^T*w_2=w_1^TE_2^TF_2v_2 \\E_2w_1 = (E_1- t_1p_1^T)*w_1 = t_1-t_1p_1^Tw_1 = 0 \\t_1^T*t_2=(t_1)^T((E_1-t1*p_1^T)w_2) = (t_1^T*E_1-t_1^T*t_1*p_1^t)w_2 = 0$

可以证明，w和t都两两正交,这里暂时不展开。从上式，我们可以知道，要建立模型，只需要求出t，p，r,其余无关的不需要计算。

最终我们要建立X，y的模型，对上式进行整理

$\small \\ E_1 = E_1 \\ E_2 = E_1-t_1p_1^T = E_1-E_1w_1p_1^T=E_1(I-w_1p_1^T) \\ E_n = E_{n-1}-t_{n-1}p_{n-1}^T= E_{n-1}(I-w_{n-1}*p_{n-1}^T) \\=E_1\prod_i^{n-1}(I-w_ip_i^T) \\t_n = E_nw_n=E_1\prod_i^{n-1}(I-w_ip_i^T)w_n$

令

$z_n=\prod_i^{n-1}(I-w_ip_i^T)w_n\Rightarrow t_n = E_1z_n\Rightarrow F_1 = E_1ZR^T$

由上式可以得到 $E_1$ 关于 $F_1$ 的回归系数 $B$ ，即

$B = ZR^T$

这里，已经建立了标准化后的F和E的关系，可以通过简单的转换，变成标准化前的参数,参考下面的MATLAB代码

代码实战

代码改自参考条目1中的代码，但基本重新实现了一遍

下面，用一个实例去解释这段代码运行过程。

function [beta] = myplsregress1(X,Y,ncomp)
%{
    Input：
    X，Y --------输入
    ncomp-------潜在因子，提取的成分数量
    output:
    beta---------模型系数
    中间变量解释
    E,F----------X,Y标准后的数据
    w  ----------E的投影轴
    t  ----------E的得分
    p------------E的回归系数
    r -----------F的回归系数
    chg----------残差变换矩阵
    z = chg*w-----用于计算得分t = E*z
   大写代表矩阵，小写代表向量
   同时也约定，W代表由w的矩阵
%}
[E,MU_x,SIGMA_x]=zscore(X); [F,MU_y,SIGMA_y] = zscore(Y); %数据标准化,变量记做 E和 F
n=size(X,2);m=size(Y,2); %n 是自变量的个数,m 是因变量的个数
num=size(X,1);%求样本点的个数
chg=eye(n); %残差转换矩阵
for i=1:ncomp
    matrix=E'*F*F'*E;
    [U1,S1,V1]=svd(matrix);
    w=V1(:,1); %提出最大特征值对应的特征向量
    Z(:,i)=chg*w; %计算z
    t=E*w; %计算成分 t
    p=E'*t/(t'*t); 
    r= F'*t/(t'*t); 
    R(:,i) = r;
    chg=chg*(eye(n)-w*p'); %计算残差矩阵的变换矩阵
    E=E-t*p';              %可以直接用chg变换阵
    F = F-t*r';
end
B = Z*R';
%系数变换
for i=1:m
    b =    MU_y(1,i)-sum(B(:,i)'.*MU_x./SIGMA_x)*SIGMA_y(1,i);
    coeff = B(:,i)'./SIGMA_x*SIGMA_y(1,i);
    beta(:,i) =[b coeff]';
end
end

除掉注释，没几条代码，枉费前面那一大堆推理。

测试

load spectra
[xl,yl,xs,ys,beta,pctvar,mse] = plsregress(NIR,octane,10);
plot(octane,'g');hold on; 
plsfit = [ones(size(NIR,1),1) NIR]*beta;
plot(plsfit ,'r');
beta = myplsregress1(NIR,octane,10);
myplsfit= [ones(size(NIR,1),1) NIR]*beta;
plot(myplsfit ,'b');