偏最小二乘法算法

偏最小二乘法 1.1 基本原理

偏最小二乘法（PLS ）是基于因子分析的多变量校正方法，其数学基础为主成分分析。但它相对于主成分回归（PCR ）更进了一步，两者的区别在于PLS 法将浓度矩阵Y 和相应的量测响应矩阵X 同时进行主成分分解：

X=TP+E Y=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=(T T T -1)T T U

因此，偏最小二乘法的校正步骤包括对矩阵Y 和矩阵X 的主成分分解以及对关联矩阵B 的计算。 1.2主成分分析

主成分分析的中心目的是将数据降维，以排除众多化学信息共存中相互重叠的信息。他是将原变量进行转换，即把原变量的线性组合成几个新变量。同时这些新变量要尽可能多的表征原变量的数据结构特征而不丢失信息。新变量是一组正交的，即互不相关的变量。这种新变量又称为主成分。

如何寻找主成分，在数学上讲，求数据矩阵的主成分就是求解该矩阵的特征值和特征矢量问题。下面以多组分混合物的量测光谱来加以说明。假设有n 个样本包含p 个组分，在m 个波长下测定其光谱数据，根据比尔定律和加和定理有：

A n×m =C n×p

B p×m

如果混合物只有一种组分，则该光谱矢量与纯光谱矢量应该是方向一致，而大小不同。换句话说，光谱A 表示在由p 个波长构成的p 维变量空间的一组点（n 个），而这一组点一定在一条通过坐标原点的直线上。这条直线其实就是纯光谱b 。因此由m 个波长描述的原始数据可以用一条直线，即一个新坐标或新变量来表示。如果一个混合物由2个组分组成，各组分的纯光谱用b1，b2表示，则有：

1122

T T T i i i a c b c b =+ 有上式看出，不管混合物如何变化，其光谱总可以用两个新坐标轴b1,b2来表示。因此可以

推出，如果混合物由p 个组分组成，那么混合物的光谱就可由p 个主成分轴的线性组合表示。

因而现在的问题就变成了如何求解这些主成分轴。而寻找这些坐标轴的基本原则是使新坐标轴包含原数据的最大方差。即沿着新坐标轴的方向，使方差达到最大。而其他方向，使方差达到最小。从几何角度看，就是变量空间中所有的点到这个新坐标轴的距离最短。

以二维空间的为例说明如何寻找主成分坐标轴。变量空间的每一个数据点（一个样本）都可以用通过该点与坐标原点的一个矢量x i 表征。

x2

x1

上图中直角三角形的三个边长分别以a,b,c 表示，那么这n 个点到第一个主成分轴v 1距离的平方和可以通过勾股定理与矢量点积得出：

222

1

1

()n

n

i i i i i b c a ===-∑∑

因为2||||i i c x =与 11||||||||cos T i i x v x v θ=??，所以

22

211

1(||||

())n n

T i

i

i i i b x x v ===-∑∑

2

2111||||()n

n

T i i i i x x v ===-∑∑

2

1111

||||()()n

n

T

T i i i i i x v x x v ===-∑∑

2111

||||n

T T i i x v X Xv ==-∑ min 上式等价于

11T

T v X Xv max （最大特征值λ）

上式中v 1表示第一个主成分轴矢量，即第一个特征矢量，所对应的最大值称为特征值，

用λ1表示。从上面推导看出，寻找主成分轴就是求X 矩阵的协方差矩阵X T X 中的最大特征

值（λi ）和特征向量（v i ）。

下面考虑变量数为m 的一般情况。 在m 为空间中新变量可以表示为：

11111221221122221122i i i m im i i i m im im m i m i mm im

u v x v x v x u v x v x v x u v x v x v x =+++=+++=+++

其中系数矩阵V 为

V =

1112121221112

m

m m mm v v v v v v v v v ??

??

????

?

?

??

用u 和x 分别表示新变量和原始矢量，则

12

i i im u u u u ??????=?????? ，12i i im x x x x ??????=??????

u V x =?

上述m 维主成分系数必须满足下面两个条件

（1） 正交条件：任意两个主成分u k 、u r ，其系数的乘积之和为0。

11220k r k r km rm v v v v v v +++=

（2）归一化条件：对于任一主成分系数的平方和等于1。

222121k k km v v v +++=

满足这两个条件的矩阵，称之为正交矩阵。正交矩阵具有如下性质：

1,T T V V I V V -==

1.3 矩阵的主成分分解

根据特征向量和特征值的定义

,1,2,,T T i i i v X Xv i m λ== (*)

同时令X 的协方差矩阵为

T Z X X =

(*)式两边同时左乘v i ，有

,1,2,,i i i Zv v i m λ==

主成分系数矩阵V 也可写为

12(,,,)m V v v v =

因此可得

{}i ZV V diag λ=?

其中{}i diag λ表示一个对角矩阵，即对角线元素为i λ，非对角线元素为0的矩阵。 上式两边同时左乘V T ，得

{}T i V ZV diag λ=

(){}T T T T i V ZV V X XV XV XV diag λ===

令T XV =,则上式变为{}T i T T diag λ=将式T XV =右乘1

V -得

T X TV =

上式是矩阵X 的主成分分解的一种表达式，由上式得求解T 和V 的方法

1()T T T V T T T X -= 1()T T XV V V -=

依据矩阵乘法规则即可获得矩阵V 和T 中每一个矢量的计算公式：

/,/T T T T

j j j j j j v t X t t t Xv v t ==

根据上面两个公式可以设计主成分分解的迭代法算法如下：

(1) 取X 中任意一列作为起始的t 。 (2) 由此t 计算：/T T T

v t X t t =

(3) 将v T 归一化：/T T T

new old old v v v =

(4) 计算新的t ：/T

t Xv v v =

(5) 比较步骤4所得的t 和上一步的t 。若二者相等（在给定的误差范围内），则按（T

t t λ=）计算特征值，转第六步继续进行；否则返回第二步继续迭代。

(6) 从Y 中减去T

t v ?的贡献：T

X X t v =-?。返回1，继续运行，直到最后Y 趋近于零。 从理论上讲，在m 空间中，可以获得m 个主成分。但是在实际应用中一般只取前几个对方差贡献最大的主成分，这样就使高维空间的数据降到低维，如二维或三维空间，非常有益于数据的观察，同时损失的信息量还不会太大。取前p 个主成分的依据为

比率(%)1

1

/p m

i

i

i i λλ===

∑∑

一般推荐，比率(%)≥80% 1.4偏最小二乘法算法

(1) 矩阵X 和Y 的标准化处理

(2) 取Y 中任意一列赋给作为起始的u 对于X 矩阵 (3) w T =u T X/u T u

(4) 归一化：/||||T T T

new old old w w w =

(5) 计算新的t ：t=Xw/w T w 对于Y 矩阵 (6) q T =t T Y/t T t

(7) 归一化：/||||T T T

new old old q q q =

(8) u=Yq/q T q 收敛判据：

(9) 将步骤8所得的u 与前一次迭代的结果相比较，若等于（在允许误差范围内）,到步骤

10，否则返回3。 (10) p T =t T X/t T t

(11) 归一化：/||||T T T

new old old p p p =

(12) t new = t old ·|| p old ||

(13) /||||T T T new old old w w p =

计算回归系数b 以用于内部关联： (14) b=u T t/t T t 对于主成分h 计算残差： (15) 1T

h h h h E E t p -=-

(16) 1T

h h h h h F E b t wq -=-

之后回到步骤(2),去进行下一主成分的运算，直到残差趋近于零。 未知样品预测

(17) 如校正部分，将X 矩阵标准化 (18) h=0，y=0

(19) h=h+1,T h h t XW =,T h h h y y b t wq =+,T

h h x x t p =-

(20) 若h>a(主成分数)，到步骤(21)。否则返回步骤(19)

(21) 得到的Y 已经标准化，因此需要按标准化步骤的相反操作，将之复原到原坐标注意的是对预测集进行标准化处理的时，使用的是训练集的均值和标准偏差。因此，在进行反标准化操作时，使用的也应该是训练集的均值和标准偏差。

1.5 程序框图与程序代码程序框图：

用C语言编制的PLS程序源代码：

C语言的源程序的各函数和功能如下：

data_input(); /*数据的输入*/

data_standardization (); /*数据的标准化*/

pctest() /*在主成分分解时判断t 的收敛*/ principalcomponent(); /*主成分分解*/

factornum(); /*确定必要的主成分分数*/ normaleq(); /*正规方程的建立*/

test(int k) /*迭代求解系数矩阵的收敛检验*/ iteration(); /*迭代法求解*/

predict(); /*未知样品预测*/

bias(); /*计算预测标准偏差*/

report(); /*输出结果*/

main(argc,argv) /*主函数*/

PLS程序源代码

#include

#include

#include

#include

#define N 25

#define M 11

#define L 4

#define TN 8

#define H 5

FILE *infp,*outfp;

double train_x[N][M], train_y[N][L];

double test_x[TN][M], test_y[TN][L];

double avg_x[M],std_x[M]；

double avg_y[L],std_ y[L]；

double error_value[TN][L]

double u[N],new_u[N];

double

save_p[H][N],save_q[H][L],save_w[H][M],save_d[H];

double predictvalue[TN][L],bias[L];

void data_input() /*数据的输入*/

{ int i,j;

for(j=0;j

for(i=0;i

fscanf(infp,“%lf”，&train_x[i][j]);

for(i=0;i

for(j=0;j

fscanf(infp,“%lf”，&train_y[i][j]);

for(j=0;j

for(j=0;j

new_u[i]+=train_y[i][j]*q[j];

new_u[i]/=sq;

}

mask=pctest();

}

while(mask);

for(j=0;j

save_q[ih][j]=q[j];

for(st=0,i=0;i

st+=t[i]*t[i];

for(j=0;j

{

p[j]=0.0;

for(i=0;i

p[j]+=t[i]*train_x[i][j]；

p[j]/=st;

}

for(sp=0,j=0;j

sp+=p[j]* p[j];

sp=sqrt(sp);

for(j=0;j

{ p[j]/=sp;

save_p[ih][j]=p[j];

}

for(i=0;i

t[i]*=sp;

for(i=0;i

{ w[i]*=sp;

save_w[ih][i]=w[i];

for(i=0;i

fscanf(infp,“%lf”，&test_x[i][j]);

for(i=0;i

for(j=0;j

fscanf(infp,“%lf”，&test_y[i][j]);

fclose(infp);

}

void data_standardization () /*数据的标准化*/

{

int i,j;

for(j=0;j

{ avg_x[j]=std_x[j]=0.0;

for(i=0;i

avg_x[j]+=train_x[i][j];

avg_x[j]/=N；

for(i=0;i

std_x[j]+=pow((train_x[i][j]-avg_x[j]),2); std_x[j]=sqrt(std_x[j])/(N-1);

}

for(i=0;i

for(j=0;j

train_x[i][j]= (train_x[i][j]-avg_x[j])/std_x[j]; for(i=0;i

for(j=0;j

train_x[i][j]= (train_x[i][j]-avg_x[j])/std_x[j]; for(j=0;j

{ avg_y[j]=std_y[j]=0.0;

for(i=0;i

avg_y[j]+=train_y[i][j];

avg_y[j]/=N;

for(i=0;i

std_y[j]+=pow((train_y[i][j]-avg_y[j]),2); std_y[j]=sqrt(std_y[j])/(N-1);

}

for(i=0;i

for(j=0;j

train_y[i][j]= (train_y[i][j]-avg_y[j])/std_y[j]; }

pctest() /*检验收敛*/

{ int i,j;

double count=0.0;

for(i=0;i

{ count+=pow((new_u[i]-u[i]),2);

}

for(st=0,i=0;i

st+=t[i]*t[i];

for(b=0,i=0;i

b+=new_u[i]*t[i]；

b/=st;

save_d[ih]=b;

for(i=0;i

for(j=0;j

train_x[i][j]-=t[i]*p[j];

for(i=0;i

for(j=0;j

train_y[i][j]-=b*t[i]*q[j];

ih++;

}

while(ih

}

predict() /*计算校正集的浓度*/ {

int i,j,k;

double pt[TN]

for(i=0;i

for(j=0;j

predictvalue[i][j]=0.0;

for(k=0;k

{ for(i=0;i

{ pt[i]=0.0;

for(j=0;j

pt[i]+=test_x[i][j]*save_w[k][j];

}

for(i=0;i

for(j=0;j

predictvalue[i][j]+=save_d[k]*pt[i]*save_q[k][j]; for(i=0;i

for(j=0;j

test_x[i][j]-=pt[i]*save_p[k][j];

}

for(i=0;i

for(j=0;j

test_x[i][j]-=pt[i]*save_p[k][j];

}

for(i=0;i

for(j=0;j

{ predictvalue[i][j]*=std_y[i];

predictvalue[i][j]+=avg_y[j];

}

}

statistics() /*计算结果评估*/

u[i]=new_u[i];

}

count=sqrt(count);

if(count<1e-12)

return 0;

else

return 1;

}

calibration() /*模型的建立*/

{ int i,j,ih,mask;

double su,sw,st,sq,sp,b;

double t[N],w[M],p[M],q[L];

ih=0;

do{ for(i=0;i

u[i]=train_y[i][l];

do{ for(su=0,i=0;i

su+=u[i]*u[i];

for(j=0;j

{ w[j]=0.0;

for(i=0;i

w[j]+=u[i]*train_x[i][j];

w[j]/=su;

}

for(sw=0,j=0;j

sw+=w[j]*w[j];

sw=sqrt(sw);

for(j=0;j

w[j]/=sw;

for(sw=0,j=0;j

sw+=w[j]*w[j];

for(i=0;i

{ t[i]=0.0;

for(j=0;j

t[i]+=train_x[i][j]*w[j];

t[i]/=sw;

}

for(st=0,i=0;i

st+=t[i]*t[i];

for(j=0;j

{ q[j]=0.0;

for(i=0;i

q[j]+=t[i]*train_y[i][j];

{ int i,j;

for(i=0;i

for(j=0;j

error_value[i][j]=test_y[i][j]-predictvalue[i][j]; for(j=0;j

{ bias[j]=0.0;

for(i=0;i

bias[j]/=TN;

bias[j]=sqrt(bias[j]);

}

}

report() /*输出结果*/

{ int i,j;

fprintf(outfp,“预测的浓度值为：\n”);

for(i=0;i

for(j=0;j

fprintf(outfp, “%9.4lf”, predictvalue[i][j]);

fprintf(outfp,“\n”);

}

fprintf(outfp,“预测结果的误差为：\n”);

for(i=0;i

{

for(j=0;j

fprintf(outfp, “%9.4lf”, error_value[i][j]);

}

fprintf(outfp,“\n”);

fprintf(outfp,“相对于各组分的校正标准偏差为\n”); for(i=0;i

fprintf(outfp, “%10.4lf”, bias[i]);

}

void main(argc,argv)

int argc；

char*argv[];

{ if(argc!=3)

{ printf(“you’v forgot enter a filename.\n”);

exit(1);

if((infp=fopen(argv[1],“r”))==0

{ printf(“can’t open data file\n”);

q[i]=/st;

}

for (sq=0, j=0;j

sq+=q[j]*q[j];

sq=sqrt(sq);

for(j=0;j

q[j]=/sq;

for (sq=0, j=0;j

sq+=q[j]*q[j];

for(i=0;i

{ new_u[j]=0.0;

exit(1);}

if((outfp=fopen(argv[2],“w”)) = = 0)

{ printf(“can’t open output file\n” ); exit(1);

data_input();

data_standardization(); calibration();

predict();

statistics();

report();

}