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

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

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

偏最小二乘法原理

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

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

2.计算投影轴w，v,使得Conv(E*w,F*v)最大

$max <Conv(E*w,F*v)>$

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

$Conv(Ew,Fv) = v^TF^TEw \ s.t. \left \| w \right \| = \left \| v \right \| = 1$

拉格朗日乘子法

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

对v，w求偏导数

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

上面两式分别乘以 wT，vT 得到

$\\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$

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

$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 是矩阵 $E^TFF^TE$ 最大特征值对应的向量， $v =F^TEw /\lambda$ ，因此只需要计算w就足够了。

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

$\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

4 计算残差

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

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

5 建模

$\\ 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的模型，对上式进行整理

$\\ 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$

求得系数

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

代码实战

代码改自参考条目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:ncompmatrix=E'*F*F'*E;[U1,S1,V1]=svd(matrix);w=V1(:,1); %提出最大特征值对应的特征向量Z(:,i)=chg*w; %计算zt=E*w; %计算成分 tp=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:mb =    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');

效果上和matlab自带的plsregress接近，但肯定会有差别，后续的小目标是实现一个和matlab自带plsregress一样的版本，最终还是实现数种简化版本的PLS以及相应过程的推导。

PLS潜在变量构造的几何解释

OLS,PCA,CCA,PLS和CR的关系总结及几何解释

参考：

【建模应用】PLS偏最小二乘回归原理与应用 - pigcv - 博客园

《数学建模算法与引用》