质量守恒约束下的多变量线性拟合算法
算法背景
之前介绍了有边界的线性拟合算法，《有边界的线性拟合算法》，它的约束条件是所有系数都须在[0,1]范围内，但实际生产环境不仅仅是这个基本的限制。
已知条件
自变量矩阵X，在生产环境中称为进料矩阵：
X是k*m矩阵，k是样本数，m是自变量数，xij是第i个样本的第j个自变量。
因变量矩阵Y，在生产环境中称为出料矩阵：
Y是k*n矩阵，k是样本数，n是因变量数，yij是第i个样本的第j个因变量。
基础系数矩阵B，在生产环境中称为基础收率：
B是n*m矩阵。bij是第i个因变量（进料）第j个自变量（出料）的基础收率，基础收率的列和等于1。
边界矩阵BD
BD是n*m矩阵。wdij是第i个因变量第j个自变量的下边界，wuij是第i个因变量
算法目标
利用X，Y，B，BD在约束条件下拟合出系数矩阵W，生产环境中称为收率矩阵：
其中W是n*m矩阵。wij是第i个因变量（进料）第j个自变量（出料）的收率。
约束条件
1. 所有系数属于[0,1]区间，生产环境中收率不可以小于0或者大于1：
wij∈[0,1]
2. 系数矩阵W的各列和为1，生产环境中物料守恒，进料等于出料：
3. 系数矩阵W不可以偏离基础系数矩阵太远，只允许在一个相对较小的范围内取值，生产环境中实际收率只允许在基础收率附近波动。
wbij+β1≤wij≤wbij+β2
其中|β1|，|β2|相对于wbij是一个比较小的数。
算法原理
把自变量矩阵与基础系数矩阵的乘积称为基础因变量矩阵，记为Y’，把因变量矩阵与基础因变量矩阵的差称为误差矩阵，记为E。
Y’=X*BT
E=Y-Y’
利用有边界的线性拟合算法拟合误差矩阵，得到偏差系数矩阵，记为WE，再通过线性变换将偏差系数矩阵的列和调整为0，将变换后的偏差系数矩阵记为WE’，最后加上基础系数矩阵即可得到系数矩阵W。
WdT =linefit_bd(X,E,BD)
Wd’= linetrans(Wd,B)
W=B+Wd’
其中linefit_bd()是由边界的线性拟合函数，BD是允许的边界矩阵；linetrans()是线性变换函数。
如何满足约束条件：
有边界的线性拟合算法将Wd控制在一个较小的范围内，使得最后的系数W不会偏离B太远；基础系数B中的所有系数都在[0,1]之间且列和为1，线性变换时会根据B来变换，将其控制在[0,1]之间，同时使Wd’的列和为0，保证W的所有系数都在[0,1]之间且列和为1；
算法过程
E=Y-X*BT
2. 有边界线性拟合算法拟合偏差系数Wd：
WdT=linefit_bd(X,E,BD)
linefit_bd()算法过程参见《有边界的线性拟合算法》。
3. 对Wd线性变换得到Wd’
Wd各列和可能不为0，对其进行线性变化，设Wd是如下矩阵：
（1） 计算Wd的各列和Wcl：
（2） 把Wd每一行看作一个向量，第i行向量为Wdi
Wdi=[wdi1 wdi2 … wdim]
偏差系数矩阵可表示成：
（3） 计算Wdi以外的其他向量的和组成的矩阵Wds
Wds是n*m矩阵，其第i个向量Wdsi=Wd1+Wd2+…+Wdi-1+Wdi+1+...+Wdm，是Wd中去掉第i个向量以外的其他向量之和。
（4） 线性变换得到Wd’
Wdi，Wdsi，Wcl三者有如下关系：
Wdi+Wdsi=Wcl
已知Wdi是准确的，则Wdsi=Wcl-Wdi也是准确的，进行如下线性变换：
其中Wdci=[wdci1 wdci2 ... wdcim] 是第i个因变量，m个自变量线性变换后的向量矩阵。
线性变换后，相当于各列的偏差系数减去该列的系数均值。
W=Wd’+B
即wij=wd’ij+bij。
如果W的各列和符合约束条件wij∈[0,1]，则拟合结束，线性变换后的偏差矩阵就是Wd’。
但当存在wd’ij使得wij∉[0,1]时，需要调整变换规则。
假设Wd的第j列Wdj有z个元素不符合约束条件或者等于0。
其中wdi(p)j是第i(p)个因变量，第j个自变量的偏差系数，它不符合约束条件或者等于0，即：
wdi(p)j+ bi(p)j∉[0,1] or wdi(p)j =0
这些元素组成的矩阵是Wdpj：
计算Wdsj：
其中wdsij是第j列第i个元素以外的其他元素的和，即：
wdsij=wd1j+wd2j+…+wd(i-1)j+wd(i+1)j+wdnj
变换方程如下：
其中Wdcj是Wdc的第j列元素，a是变换矩阵，是n*1矩阵。
为了使Wdcj满足约束条件，当wdij不满足约束条件时，ai=0，解上述方程如下：
其中z是不符合约束条件或者等于0的偏差系数个数，
是z个不符合约束条件或者等于0的偏差系数之和，ai是a的第i个元素。
求得变换矩阵a后，带入Wdcj=Wdj-Wdsj**a，求得Wdcj。
用上述方法求得所有列的Wdcj，也就得到Wd’。
如果Wd’的所有元素符合约束条件wij=(wdcij+bij)∈[0,1]，则拟合结束，变换后的偏差矩阵就是Wd’。如果仍然有wij∉[0,1]，则把对应位置的wdcij元素视为不符合约束条件或者等于0的元素，加到Wdpj中，递归上述方法，直到wij=(wdcij+wbij)∈[0,1]为止，得到Wd’。
变换后的偏差系数Wd’加基础系数矩阵B，得到W
W=B+Wd’
写出代码
A |
B |
C |
D |
E |
|
1 |
=X |
/参数：X 矩阵 (k*m) |
|||
2 |
=Y |
/参数：Y 矩阵 (k*n) |
|||
3 |
=B |
/参数：基础系数矩阵 (n*m) |
|||
4 |
=BD |
/参数：边界矩阵 BD(n*m) |
|||
5 |
=linefit_bd_file |
/边界算法文件绝对路径，如：”D:/linefit_bd_file.dfx” |
|||
6 |
=func(A8,A1,A2,A3,A4) |
||||
7 |
/约束条件下的线性变换，参数：X,Y,B,BD |
||||
8 |
func |
||||
9 |
=transpose(C8) |
/BT,B转置 |
|||
10 |
=mul(A8,B9) |
/X*BT |
|||
11 |
=B8.(~--B10(#)) |
/E=Y-X*BT |
|||
12 |
=call(linefit_bd_file,A8,B11,D8) |
/WdT=linefit_bd(X,E,BD) |
|||
13 |
=transpose(B12) |
/Wd |
|||
14 |
=func(A17,B13,C8) |
/WT=linetrans(Wd,B) |
|||
15 |
=transpose(B14) |
/W |
|||
16 |
/线性变换，参数：Wd，B |
||||
17 |
func |
||||
18 |
=transpose(A17) |
/Wd转置 |
|||
19 |
=transpose(B17) |
/B转置 |
|||
20 |
=A17.pselect@a(~==~.len()*[0]) |
/Wd全为 0 的行 |
|||
21 |
=B18.(~.count(~==0)) |
/Wd每列 0 的个数 |
|||
22 |
=A17.len() |
||||
23 |
=B18.((s=~.sum(),ind=#,~.(if(~==0,0,s/(B22-B21(ind)))))) |
/0不参与平均，各列求均值 |
|||
24 |
=B18.(~--B23(#)) |
/各列值减平均，Wd' |
|||
25 |
=B19.(~++B24(#)) |
/W=B+Wd' |
|||
26 |
=transpose(B25) |
||||
27 |
=B26.pselect@a(~.pselect(~<0||~>1)>0) |
/查看是否有超范围的系数 |
|||
28 |
if B27.len()==0 |
return B25 |
/没有超范围的返回 W |
||
29 |
else |
=func(A32,B18,B20,B27,B19) |
/有超范围的变换规则 |
||
30 |
return C29 |
||||
31 |
/递归消除不符合要求的系数，参数：WdT, 全为 0 的行号, 超范围的 Y 系数索引,BT |
||||
32 |
func |
||||
33 |
=D32.~.len() |
||||
34 |
=A32.(~.sum()) |
/Wcl |
|||
35 |
=A32.((w=~,w.((ind=#,if(~==0,0,w(to(B33)\ind).sum()))))) |
/Wds转置 |
|||
36 |
=A32.(~.pselect@a(~==0)) |
||||
37 |
for C32 |
=B32&B37 |
/全为 0 的行与超限行合并 |
||
38 |
=B36.(~&C37) |
/Wdpj的索引 |
|||
39 |
=B34.(~/((B33-C38(#).len()-1)*~+A32(#)(C38(#)).sum())) |
/Wclj/[(n-z-1)*Wclj+wdi(p)j.sum()] |
|||
40 |
=A32.((ind=#,~.(if(C38(ind).pos(#)>0,0,C39(ind))))) |
/a |
|||
41 |
=B35.(~**C40(#)) |
/a**Wdsj |
|||
42 |
=A32.((idx=#,if(C40(#)==C40(#).len()*[0],C41(#),~--C41(#)))) |
/Wd‘转置 |
|||
43 |
=D32.(~++C42(#)) |
/BT+Wd'T |
|||
44 |
=transpose(C43) |
/W |
|||
45 |
=C44.pselect@a(~.pselect(~<0||~>1)>0) |
/再次确认是否有超限系数 |
|||
46 |
if C45.len()==0 |
return C43 |
/没有超范围的返回 W |
||
47 |
else |
=B32=C37 |
/有超范围的记录下超限索引 |
||
48 |
=func(A32,A32,B32,C45,D32) |
/递归变换 Wd |
|||
输入参数说明：
1. A1中的X是输入自变量矩阵，形状是k*m。
2. A2中的Y是输入自变量矩阵，形状是k*n。
3. A3中的B是输入基础系数矩阵，形状是n*m。
4. A4中的BD是输入边界矩阵，形状是n*m（这里要注意，每一行是1个Y的所有X的边界）。
应用举例
自变量矩阵X：
因变量矩阵Y：
基础系数矩阵B：
边界矩阵BD：
算法运行结果，最佳系数矩阵W：
最佳系数矩阵W的rmse：
rmse1=1.64
rmse2= 2.72
rmse3= 1.93
rmsei是第i列Y的rmse。
再来看下最小二乘法直接拟合的结果W：
完全不满足约束条件。
