偏最小二乘法算法

偏最小二乘法1.1基本原理偏最小二乘法（PLS）是基于因子分析的多变量校正方法，其数学基础为主成分分析。但它相对于主成分回归（PCR）更进了一步，两者的区别在于PLS法将浓度矩阵Y和相应的量测响应矩阵X同时进行主成分分解：X=TP+EY=UQ+F式中T和U分别为X和Y的得分矩阵，而P和Q分别为X和Y的载荷矩阵，E和F分别为运用偏最小二乘法去拟合矩阵X和Y时所引进的误差。偏最小二乘法和主成分回归很相似，其差别在于用于描述变量Y中因子的同时也用于描述变量X。为了实现这一点，数学中是以矩阵Y的列去计算矩阵X的因子。同时，矩阵Y的因子则由矩阵X的列去预测。分解得到的T和U矩阵分别是除去了大部分测量误差的响应和浓度的信息。偏最小二乘法就是利用各列向量相互正交的特征响应矩阵T和特征浓度矩阵U进行回归：U=TB得到回归系数矩阵，又称关联矩阵B：B=(TTT-1)TTU因此，偏最小二乘法的校正步骤包括对矩阵Y和矩阵X的主成分分解以及对关联矩阵B的计算。1.2主成分分析主成分分析的中心目的是将数据降维，以排除众多化学信息共存中相互重叠的信息。他是将原变量进行转换，即把原变量的线性组合成几个新变量。同时这些新变量要尽可能多的表征原变量的数据结构特征而不丢失信息。新变量是一组正交的，即互不相关的变量。这种新变量又称为主成分。如何寻找主成分，在数学上讲，求数据矩阵的主成分就是求解该矩阵的特征值和特征矢量问题。下面以多组分混合物的量测光谱来加以说明。假设有n个样本包含p个组分，在m个波长下测定其光谱数据，根据比尔定律和加和定理有：An×m=Cn×pBp×m如果混合物只有一种组分，则该光谱矢量与纯光谱矢量应该是方向一致，而大小不同。换句话说，光谱A表示在由p个波长构成的p维变量空间的一组点（n个），而这一组点一定在一条通过坐标原点的直线上。这条直线其实就是纯光谱b。因此由m个波长描述的原始数据可以用一条直线，即一个新坐标或新变量来表示。如果一个混合物由2个组分组成，各组分的纯光谱用b1，b2表示，则有：1122TTTiiiacbcb有上式看出，不管混合物如何变化，其光谱总可以用两个新坐标轴b1,b2来表示。因此可以推出，如果混合物由p个组分组成，那么混合物的光谱就可由p个主成分轴的线性组合表示。因而现在的问题就变成了如何求解这些主成分轴。而寻找这些坐标轴的基本原则是使新坐标轴包含原数据的最大方差。即沿着新坐标轴的方向，使方差达到最大。而其他方向，使方差达到最小。从几何角度看，就是变量空间中所有的点到这个新坐标轴的距离最短。以二维空间的为例说明如何寻找主成分坐标轴。变量空间的每一个数据点（一个样本）都可以用通过该点与坐标原点的一个矢量xi表征。主成分轴x2x1v1byicaθ0上图中直角三角形的三个边长分别以a,b,c表示，那么这n个点到第一个主成分轴v1距离的平方和可以通过勾股定理与矢量点积得出：22211()nniiiiibca因为2||||iicx与11||||||||cosTiixvxv，所以222111(||||())nnTiiiiibxxv22111||||()nnTiiiixxv21111||||()()nnTTiiiiixvxxv2111||||nTTiixvXXvmin上式等价于11TTvXXvmax（最大特征值λ）上式中v1表示第一个主成分轴矢量，即第一个特征矢量，所对应的最大值称为特征值，用λ1表示。从上面推导看出，寻找主成分轴就是求X矩阵的协方差矩阵XTX中的最大特征值（λi）和特征向量（vi）。下面考虑变量数为m的一般情况。在m为空间中新变量可以表示为：11111221221122221122iiimimiiimimimmimimmimuvxvxvxuvxvxvxuvxvxvx其中系数矩阵V为V=1112121221112mmmmmvvvvvvvvv用u和x分别表示新变量和原始矢量，则12iiimuuuu，12iiimxxxxuVx上述m维主成分系数必须满足下面两个条件（1）正交条件：任意两个主成分uk、ur，其系数的乘积之和为0。11220krkrkmrmvvvvvv（2）归一化条件：对于任一主成分系数的平方和等于1。222121kkkmvvv满足这两个条件的矩阵，称之为正交矩阵。正交矩阵具有如下性质：1,TTVVIVV1.3矩阵的主成分分解根据特征向量和特征值的定义,1,2,,TTiiivXXvim(*)同时令X的协方差矩阵为TZXX(*)式两边同时左乘vi，有,1,2,,iiiZvvim主成分系数矩阵V也可写为12(,,,)mVvvv因此可得{}iZVVdiag其中{}idiag表示一个对角矩阵，即对角线元素为i，非对角线元素为0的矩阵。上式两边同时左乘VT，得{}TiVZVdiag(){}TTTTiVZVVXXVXVXVdiag令TXV,则上式变为{}TiTTdiag将式TXV右乘1V得TXTV上式是矩阵X的主成分分解的一种表达式，由上式得求解T和V的方法1()TTTVTTTX1()TTXVVV依据矩阵乘法规则即可获得矩阵V和T中每一个矢量的计算公式：/,/TTTTjjjjjjvtXtttXvvt根据上面两个公式可以设计主成分分解的迭代法算法如下：(1)取X中任意一列作为起始的t。(2)由此t计算：/TTTvtXtt(3)将vT归一化：/TTTnewoldoldvvv(4)计算新的t：/TtXvvv(5)比较步骤4所得的t和上一步的t。若二者相等（在给定的误差范围内），则按（Ttt）计算特征值，转第六步继续进行；否则返回第二步继续迭代。(6)从Y中减去Ttv的贡献：TXXtv。返回1，继续运行，直到最后Y趋近于零。从理论上讲，在m空间中，可以获得m个主成分。但是在实际应用中一般只取前几个对方差贡献最大的主成分，这样就使高维空间的数据降到低维，如二维或三维空间，非常有益于数据的观察，同时损失的信息量还不会太大。取前p个主成分的依据为比率(%)11/pmiiii一般推荐，比率(%)≥80%1.4偏最小二乘法算法(1)矩阵X和Y的标准化处理(2)取Y中任意一列赋给作为起始的u对于X矩阵(3)wT=uTX/uTu(4)归一化：/||||TTTnewoldold(5)计算新的t：t=Xw/wTw对于Y矩阵(6)qT=tTY/tTt(7)归一化：/||||TTTnewoldoldqqq(8)u=Yq/qTq收敛判据：(9)将步骤8所得的u与前一次迭代的结果相比较，若等于（在允许误差范围内）,到步骤10，否则返回3。(10)pT=tTX/tTt(11)归一化：/||||TTTnewoldoldppp(12)tnew=told·||pold||(13)/||||TTTnewoldoldwwp计算回归系数b以用于内部关联：(14)b=uTt/tTt对于主成分h计算残差：(15)1ThhhhEEtp(16)1ThhhhhFEbtwq之后回到步骤(2),去进行下一主成分的运算，直到残差趋近于零。未知样品预测(17)如校正部分，将X矩阵标准化(18)h=0，y=0(19)h=h+1,ThhtXW,Thhhyybtwq,Thhxxtp(20)若ha(主成分数)，到步骤(21)。否则返回步骤(19)(21)得到的Y已经标准化，因此需要按标准化步骤的相反操作，将之复原到原坐标注意的是对预测集进行标准化处理的时，使用的是训练集的均值和标准偏差。因此，在进行反标准化操作时，使用的也应该是训练集的均值和标准偏差。1.5程序框图与程序代码程序框图：输入数据矩阵：训练集X，Y和预测集X’数据标准化：Xij’=(Xij-Xj)/Sj从Y中任选一列作为初始u，h=h+1wT=uTX/uTu,/||||TTTnewoldoldt=Xw/wTw,qT=tTY/tTt,/||||TTTnewoldoldqqq,u=Yq/qTqu=Yq/qTqqT=tTY/tTtt=Xw/wTw两次迭代t相等吗pT=tTX/tTt，/||||TTTnewoldoldppp，tnew=told·||pold||，/||||TTTnewoldoldwwp，b=uTt/tTt，1ThhhhEEtp，1ThhhhhFEbtwqh=H?h=h+1,ThhtXW，ThhxxtpThhhyybtwq，Thhxxtp令h=0,y=0h=H?反标准化操作，输出结果用C语言编制的PLS程序源代码：C语言的源程序的各函数和功能如下：data_input();/*数据的输入*/data_standardization();/*数据的标准化*/pctest()/*在主成分分解时判断t的收敛*/principalcomponent();/*主成分分解*/factornum();/*确定必要的主成分分数*/normaleq();/*正规方程的建立*/test(intk)/*迭代求解系数矩阵的收敛检验*/iteration();/*迭代法求解*/predict();/*未知样品预测*/bias();/*计算预测标准偏差*/report();/*输出结果*/main(argc,argv)/*主函数*/PLS程序源代码#includestidio.h#includemath.h#includeio.h#includestdlib.h#defineN25#defineM11#defineL4#defineTN8#defineH5FILE*infp,*outfp;doubletrain_x[N][M],train_y[N][L];doubletest_x[TN][M],test_y[TN][L];doubleavg_x[M],std_x[M]；doubleavg_y[L],std_y[L]；doubleerror_value[TN][L]doubleu[N],new_u[N];doublesave_p[H][N],save_q[H][L],save_w[H][M],save_d[H];doublepredictvalue[TN][L],bias[L];voiddata_input()/*数据的输入*/{inti,j;for(j=0;jM;j++)for(i=0;iN;i++)fscanf(infp,“%lf”，&train_x[i][j]);for(i=0;iN;i++)for(j=0;jL;j++)fscanf(infp,“%lf”，&train_y[i][j]);for(j=0;jM;j++)for(j=0;jL;j++)new_u[i]+=train_y[i][j]*q[j];new_u[i]/=sq;}mask=pctest();}while(mask);for(j=0;jL;j++)save_q[ih][j]=q[j];for(st=0,i=0;iN;i++)st+=t[i]*t[i];for(j=0;jM;j++){p[j]=0.0;for(i=0;iN;i++)p[j]+=t[i]*train_x[i][j]；p[j]/=st;}for(sp=0,j=0;jM;j++)sp+=p[j]*p[j];sp=sqrt(sp);for(j=0;jM;j++){p[j]/=sp;save_p[ih][j]=p[j];}for(i=0;iN;i++)t[i]*=sp;for(i=0;iM;i++){w[i]*=sp;save_w[ih][i]=w[i];for(i=0;iTN;i++)fscanf(infp,“%lf”，&test_x[i][j]);for(i=0;iTN;i++)for(j=0;