线性回归模型的降维计算：交互投影算法，本文主要内容关键词为：线性论文,算法论文,模型论文，此文献不代表本站观点，内容供学术参考，文章仅供参考阅读下载。
线性回归模型的自变量如果太多，不仅影响模型的统计性质，而且给计算带来困难。我们知道矩阵X＇X的阶数，就是线性回归模型自变量的个数，X＇X的阶数如果太高，求逆运算可能无法进行。研究线性模型的降维计算方法一直是统计计算的一个重要课题。这方面的研究成果比较多。笔者提出的交互投影算法与投影寻踪回归、偏最小二乘方法、单指标回归模型等，都是线性回归模型降维计算的行之有效的方法。
设有线性回归模型
Y[,n×1]＝X[,n×p]β[,p×1]+ε[,n×1]
（1）
它的回归系数的最小二乘估计是
β＝（X＇X）[-1]X＇Y
（2）
现在考虑求（X＇X）[-1]的某些计算困难和解决办法。我们知道，求p阶逆矩阵时乘除运算的规模与p[3]成正比， 如果能设法把求逆的阶数降下来，可以节省计算时间。在微机上使用Fortran语言时对p是有限制的，某种Fortran编译程序对200×100的矩阵就不能编译， 并提出维数太多。此时如果不实现降维计算就无法求解。分块求逆可以实现降维计算，可是它需要较多的工作数组，往往也不能通过编译。那么除了投影寻踪回归、偏最小二乘等方法外，还有别的降维计算办法吗？
还有一个问题。X的列向量可能是复共线，此时行列式│X＇X │近似于0，不仅β的最小二乘估计性质变坏， 而且计算中的扫除变换也可能停止。我们能不能把复共线的两列向量分开呢？
采用凸集间交互投影的迭代算法可以实现求线性模型回归系数的降维计算与病态分离。
下面我们介绍解决这个问题的交互投影的迭代算法。
回归模型（1）可以写为
β
‖Y－Xβ‖→min （3）
它可以被看作是求一定点Y到一个p维线性子空间Xβ，β∈R[p]的投影。现在我们重新剖分X、β。令X＝（X[,1]│X[,2]），β＇＝（β＇[,1]│β＇[,2]），有
β[,1],β[,2]
‖Y－X[,1]β[,1]－X[,2]β[,2]‖ → min（4）X[,2]β[,2]（β[,2]∈R[p2]，p[,2]+p[,1]＝p，当p为偶数时，最好取p[,1]＝p[,2]），是一线性子空间，容易证明Y－X[,1]β[,1]（β[,1]∈R[p1]），是一闭凸集（在点Y－X[,1]β[,11]与点Y－X[,1]β[,12]连线上，任取一点Y－X[,1]β[,11]+t（X[,1]（β[,11]－β[,12]）），0≤t≤1，显然仍属集合Y－X[,1]β[,1]（β[,1]∈R[p1]））。于是上式可以看作是求一闭凸集与一线性子空间的距离，我们可以使用如下的交互投影的迭代方法。
任取初值β[(1)][,1]，计算Y[(1)][,1]＝Y－X[,1]β[(1)][,1]，模型变为
β[,2]
‖Y[(1)][,1]－X[,2]β[,2]‖─→ min （5）
我们采用普通方法来计算它的参数的最小二乘估计。因为复共线的两列向量已被剖分到两个子块X[,1]和X[,2]中，即实现了病态分离，所以计算没有什么困难，求逆矩阵的维数也下降了大约二分之一。此时即可求出从闭凸集Y－X[,1]β[,1]上一点Y－X[,1]β[(1)][,1]到子空间X[,2]β[,2]上的投影X[,2]β[(1)][,2]。
从（5）中计算出来β[(1)][,2]，代入（4）中，可计算出Y[(1)][,2]＝Y－X[,2]β[(1)][,2]。模型（4）变为
β[,1]
‖Y[(1)][,2]－X[,1]β[,1]‖─→ min （6）
同样用普通方法可求得β[(2)][,1]。 于是可以得到从线性子空间X[,2]β[,2]上一点X[,2]β[(1)][,2]到闭凸集Y－X[,1]β[,1] 的投影Y－X[,1]β[(2)][,1]。
用这个方法得到序列
β[(1)][,1]，β[(1)][,1]，β[(2)][,1]，β[(2)][,2]，…，β[(k)][,1]，β[(k)][,2] （7）对于给定的控制参数ε，当
‖β[(k+1)][,1]－β[(k)][,1]‖≤ε，‖β[(k+1)][,2]－β[(k)][,2]‖≤ε时，停止迭代过程，可得到原来模型的一个渐近解
β[(k)＇]＝（β[(k)＇][,1]│β[(k)＇][,2]）（8）
现在我们需要回答两个问题。第一，序列｛β[(k)]｝是否收敛？
第二，序列｛β[(k)]｝是否收敛到原来模型的最小二解β？
下面我们给出交互投影收敛性的证明。
在一般情况下我们来让明两个闭凸集间的交互投影是收敛的。
设Ω是一个内积空间，d是在Ω中由内积定义的距离。 我们首先给出在Ω中由一点到一个闭集的投影的定义。
定义1 设点b∈Ω，闭集A
Ω，点a[,0]∈A，a∈A，如果
d（b，a[,0]）＝inf d（b，a）＝d（b，A） （9）
a∈A则称a[,0]是点b在闭集A上的投影。
很明显，如果A是一个子空间， 上述定义就和一般内积空间投影定义一致，下面再给出两个闭集间局部距离的定义。
定义2 设闭集 AΩ、BΩ；点a∈A，b∈B，如果
d（a，b）＝d（a，B）＝d（b，A）（10）则d（a，b）是闭集A和B之间的一个局部距离。
考虑两个闭凸集间的局部距离与整体距离之间的关系，有如下定理：
定理1 设A、B是内积空间中两个闭凸集，d（a，b）是A、B间的一个局部距离，则
d（A，B）＝d（a，b） （11）
证明 如果d（a，b）＝0，则定理显然正确。假定d（a，b）≠0，
─可分别通过a，b作两个超平面π[,A]和π[,B]，使之都垂直于线段ab。两个闭凸集A和B分别位于π[,A]和π[,B]的外侧。d（a，b）是π[,A]和π[,B]的距离，同时也是闭凸集A和B的距离。证毕。
现在我们叙述两个闭凸集间交互投影的过程，它可以用来寻找两个闭凸集间的距离。如果这两个闭凸集相交，它也是寻找其一个交点的方法。
设A和B是一内积空间中的两个闭凸集，任给a[,0]∈A，求得a[,0]在B点的投影b[,0]。再进一步寻求b[,0]在A上的投影a[,1]。如果继续下去，当得到a[,i]在B上的投影b[,i]后，再寻求b[,i]在A上的投影a[,i]+1。这样便可得到成对点列
a[,1]，b[,1]，a[,2]，b[,2]，…；a[,i]，b[,i] （12）
上述迭代过程收敛的含义是：存在点对a[*]∈A，b[*]∈B，i →∞，有a[,ii]→a[*]，b[,i]→b[*]，且d（a[*]，b[*]）＝d（A，B）。
定理2 两个闭凸集间的交互投影过程收敛。
证明 显然有
d（a[,0]，b[,0]）≥d（a[,1]，b[,0]）≥d（a[,1]，b[,1] ）≥…≥d（a[,i]，b[,i]）≥d（a[,i+1]，b[,i]）≥…单调递减下有界数列必有下确界，设其收敛于d（a[*]，b[*]）。因A、B是闭集，距离是可达的，a[*]∈A，b[*]∈B。因A、B是凸集，d（a[*]，b[*]）也是A、B间的整体距离，故定理成立。证毕。
如果A、B相交，则其距离为0，此时a[*]＝b[*]是A、B 的一个交点。
现在我们可以回答前面提出的两个问题。原模型（3 ）是寻求一点Y到一子空间Xβ的距离，它等价于模型（4 ）寻求一个闭凸集到一个子空间之间的距离。而这个子空间也可以看作闭凸集，它们之间的交互投影过程收敛。由于X、X[,1]、X[,2]都列满秩， 所作的最小二乘估计都是唯一的，所以前面两个问题答案都是肯定的。
一、算例 交互投影迭代算法
（一）病态分离
在模型Y＝Xβ+ε中，取n＝8，p＝6，令Y＝（4，3，3，2，2，4，3，3）’，X＝（X[,1]│X[,2]）
┌1 0 0 1 0 1 0 1┐
┌1 0 0 1 0 1 0 1┐
X'[,1]＝│0 1 1 0 1 0 1 0│， X'[,2]＝│1 1 0 0 1 1 0 0│
└0 0 1 0 0 1 1 0┘
└1 1 1 0 0 0 1 1┘
注意X[,1]与X[,2]的第一列向量相等，行列式│X'X│＝0，我们可以采用广义逆矩阵去计算参数的LSE （在这里采用的是病态分离后的交
互投影方法）。取初值β[,0]＝（1.5，1.5，1.5）＇， 算得β＝（1.625，1.0，1.0，0.375，1.0，1.0）＇。再取初值β[,0]＝（0.8，0.8，0.8）＇，算得β＝（0.75，1.0，1.0，1.25，1.0，1.0）＇。容易验证它们都是合适的，一般地，因秩rk（X）＝5，有β＝β[,1]＝tβ[,2]，组合简化得β[,1]＝（1，1，1，1，1，1），β[,2]＝（1，0，0，－1，0，0），t∈R。
（二）降维计算
在模型Y＝Xβ+ε中，取n＝200，p＝80。 注意在微机上某种Fortran编译程序不能编译一个200×100的数组，因数组太大， 单独一个200×80的数组可通过编译，但加上一个80×80 的数组就不能通过编译。用分块求逆的办法需要有两个200×40和两个40×40的数组， 合在一起也不能通过编译。我们只好利用交互投影的方法，不过需要一点编程技巧。
首先来构造一个大规模回归模型， 用发生随机数的办法作出两个200×40的矩阵X[,1]和X[,2]，1个200×1的数组Y，回归模型就算作好了。把X[,1]和X[,2]分别存入两个数据文件。主程序里只设一个200 ×40和两个40×40的工作数组，可以通过编译。每次迭代用工作数组读入数据文件，就可节省一个200×40的数组。用后除变换求X'[,1]X[,1]的逆，不需要先算出X'[,1]X[,1]，只需要一个工作数组X[,1]，这样就可以实现交互投影计算了。
为了验证计算结果，可以先令β＝（1，1，…，1）＇，发生随机数ε～N（0，σ[2]），计算出Y＝Xβ+ε。 再用交互投影方法计算得
β≈（1，1，…1）＇。计算是成功的。