一笑出门去，千里落花风

数据挖掘|基于乔利斯基分解运算的逐步回归优化

目录

本文介绍优化的逐步回归方法：基于Cholesky（乔利斯基）分解实用算法和一定矩阵计算简化方法的逐步回归。原文来自马景义老师的数据挖掘教学。

前排警告：本文变量众多，请耐心按顺序阅读研究。

理论部分：Cholesky分解法

Cholesky(乔利斯基)分解法,又称为平方根法，是求解对称正定线性方程组最常用的方法之一。Cholesky分解定理如下：

若$A \in R^ {n \times n}$对称正定，则存在一个对角元均为正数的下三角阵$L \in R^ {n \times n}$，使得
$\begin{equation}A=LL^T. \end{equation}$

因此，我们可按如下步骤求解方程组$Ax=b​$：
(1)计算$A​$的Cholesky分解：$A=LL^T​$;
(2)求解$Ly=b​$得$y​$;
(3)求解$L^Tx=y​$得$x​$.

Cholesky分解 实用算法part

$\begin{equation}A=LL^T.\end{equation}$
$A$中的元素记为$a_{ij}$
通常简单而实用的Cholesky分解方法是通过直接比较$A=LL^T$两边的对应元素来计算$L$。设

则$LL^T$的元素为

比较两边元素，有关系式：

可直接计算$L$的第一列元素：

……

进一步地，假定已算出$L$的前$k-1$列元素。由

则可得

再由

得

由此便可求出矩阵$L$。

以上介绍可能过于数学化，下面以一个三阶方阵为例进行说明。设

则$M_1=LL^T$的元素为

可以以列为单位求解$L​$的元素。根据之前的叙述，$L​$的第一列元素可以很容易计算得到。另外，注意到$M_1​$删去第一行和第一列后余下的$2 \times 2​$矩阵（记为$M_{12}​$）的每一个分量的第一项为

因此，在计算出$L​$的第一列元素后，从$M_{12}​$中减去上述矩阵，得到$M_2​$：

注意到这与$M_1​$右下部分的形式完全一致，可以继续计算$L​$的第二列，进而求出$L​$.

前/回代法（Forward/Backward Substitution）

对于方程组$Lx=b$和$Ux=b$（其中$L$和$U$分别是下三角矩阵和上三角矩阵），可以分别采用前代法和回代法求解。此处以下三角矩阵的前代法为例。

设方程组$Lx=b$的$L$是已知的非奇异下三角阵（则其主对角线元素非零），则方程组的矩阵形式为：

则由方程组的第一个方程可得：

进一步地，如果已求出$x_1,…,x_{i-1}$，就可以根据方程组的第$i$个方程

求出

从而可以求出方程组的解。
求解上三角方程组的回代法的做法类似。

氦核注：利用回代法求解矩阵的目的是优化r语言内部函数的矩阵计算，提升程序运行速度。
其中某些步骤使用反复回代，修改其中部分代码实现，逻辑比较复杂，需要数学验证。

实现代码

Cholesky分解

mchol函数

mchol函数能够实现对对称正定矩阵的Cholesky分解。它以一个对称正定方阵为参数，返回其Cholesky因子。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
mchol <- function(x)
{
  #目标是构造L矩阵，X = L*L'
  mn <- dim(x)
  m <- mn[1] #行数
  n <- mn[2] #列数
    
  #检验传入参数是否符合要求
  if(m != n) 
  {
    return ("Wrong dimensions of matrix!") #检查是否为方阵
  }
  if(sum(t(x) != x) > 0) 
  {
    return ("Input matrix is not symmetrical!") #检查是否对称
  }

  L <- matrix(0, m, m) # 构造一个L矩阵，容纳L_ij
  
  for(i in 1:m)  #以列为单位求解L
  {
    L[i,i] <- sqrt(x[i,i])
    if(i < m)
    {
      L[(i+1):m,i] <- x[(i+1):m,i]/L[i,i]

      #更新本次循环计算后的矩阵，便于下一次同样方法计算
      TLV <- L[(i+1):m,i] 
      #TLM <- matrix(TLV, m-i, m-i)
      #TLM <- sweep(TLM, 2, TLV, "*")
      TLM <- outer(TLV,TLV) #L的第一列正交
      x[(i+1):m,(i+1):m] <- x[(i+1):m,(i+1):m] - TLM #作差消去外层，得到新矩阵，再计算下一列
    }
  }
  L  
}

氦核注：这个函数妙处在于，虽然由两层循环构成，但是只需要变动一个参数i就可以完成整个下三角矩阵的计算。
下面是该函数的一个使用例子，直接出结果，没什么看的。

1
2
3
4
y=matrix(rnorm(20),5)
x=t(y)%*%y  #构造对称正定矩阵
mchol(x)
t(chol(x)) #r语言自带函数，做出来是上三角阵，我们的函数为了方便思考，得到的结果正好对称

2.3637979 0.0000000 0.0000000 0.000000
1.2068850 1.7531628 0.0000000 0.000000
-0.3977667 -0.2948317 0.9732456 0.000000
0.8956632 0.6154149 0.6515290 1.175918

2.3637979 0.0000000 0.0000000 0.000000
1.2068850 1.7531628 0.0000000 0.000000
-0.3977667 -0.2948317 0.9732456 0.000000
0.8956632 0.6154149 0.6515290 1.175918

前/回代法求解

前代法 - mforwardsolve函数

mforward()以方程组的系数矩阵L（下三角矩阵）和右端常数项b为参数，使用前代法求解方程组$Lx=b$，返回方程组的解。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
mforwardsolve=function(L,b)
{
  mn=dim(L); m=mn[1]; n=mn[2]
  if(m!=n) return ("Wrong dimensions of matrix L!")
  if(m!=length(b)) return ("Wrong dimensions of matrix L or vector b!") #检查输入参数是否符合要求
  x=rep(0,m)
  for(i in 1:m)
  {
    x[i]=b[i]/L[i,i]
    if(i<m) b[(i+1):m]=b[(i+1):m]-x[i]*L[(i+1):m,i]   
    #更新右端项，逐步减去已求出的未知量对右端项的贡献   
  }
  x  
}

下面是使用mforward()的一个例子。

1
2
3
4
5
y=matrix(rnorm(20),5)
x=t(y)%*%y
L=mchol(x); b=1:4
mforwardsolve(L,b)
forwardsolve(L,b) #R语言中自带函数的求解结果

0.510407168928626 0.954253346647569 1.24831106733918 3.88142767768075
0.510407168928626 0.954253346647569 1.24831106733918 3.88142767768075

回代法 - mbacksolve函数（与前代法结果相同，形式不同）

mbacksolve()以方程组的系数矩阵L（上三角矩阵）和右端常数项b为参数，使用回代法求解方程组，返回方程组的解。

1
2
3
4
5
6
7
8
9
10
11
12
mbacksolve=function(L,b){
  mn=dim(L); m=mn[1]; n=mn[2]
  if(m!=n) return ("Wrong dimensions of matrix L!")
  if(m!=length(b)) return ("Wrong dimensions of matrix L or vector b!") #检查输入参数是否符合要求
  x=rep(0,m)
  for(i in m:1){ 
    #这里的循环条件是i从大到小，从下向上运算，因此得到的是上三角矩阵，原理相同
    x[i]=b[i]/L[i,i]
    if(i>1) b[(i-1):1]=b[(i-1):1]-x[i]*L[(i-1):1,i]   #更新右端项   
  }
  x  
}

下面是使用mbacksolved()的一个例子

1
2
3
4
5
y=matrix(rnorm(20),5)
x=t(y)%*%y
L=mchol(x); b=1:4
mforwardsolve(L,b)
forwardsolve(L,b)  #R语言自带函数的求解结果

0.471045013164262 1.75905390934574 1.69960995026223 16.4369180950282
0.471045013164262 1.75905390934574 1.69960995026223 16.4369180950282

逐步回归

减少一个变量: $Lk$ 的计算

逐步回归时,在当前模型的基础之上减去一个变量,求解回归系数时,可以不需要重新计算新的$X_{k}^T X_{k}$的分解,而是可以借助上一步的分解$L$,使用Givens变换简化计算.具体的操作方法就是:如果现在要删掉第k个变量，那么对应$X_{k}^T X_{k}$的Cholesky因子$L_{k}$就是上一步的$L$删掉第k行,然后右乘Givens矩阵进行列变换,化为下三角矩阵$L_{k}$.下面举个例子.

Givens变换 - gives函数（先决知识）

变换的目标是
$ \left[\begin{matrix} a & b \end{matrix}\right]
\left[\begin{matrix} g_{1} & g_{2}\\ g_{3} & g_{4} \end{matrix}\right]
=\left[\begin{matrix} c & 0 \end{matrix}\right] ​$，
即通过变换，将最后一列变成0。
这个$ G ​$矩阵就是Givens变换要构造的矩阵。

通常使用三角函数矩阵$\left[\begin{matrix} cos(\theta) & -sin(\theta)\\ sin(\theta) & cos(\theta) \end{matrix}\right] ​$。
这个矩阵与自己是正交的，不影响向量的长度$ \sqrt{a^2 + b^2}​$，
如$ \left[\begin{matrix} a & b \end{matrix}\right] ​$中的a，b分别为3，4时，通过Givens变换的到的矩阵为$ \left[\begin{matrix} 5 & 0 \end{matrix}\right] ​$，

故可以构造Givens矩阵为$ \left[\begin{matrix} a/(\sqrt{a^2 + b^2}) & -b/(\sqrt{a^2 + b^2})\\ b/(\sqrt{a^2 + b^2}) & a/(\sqrt{a^2 + b^2}) \end{matrix}\right] $。

首先给出构造Givens变换矩阵的函数.

1
2
3
4
5
6
7
8
9
10
##构造Givens矩阵的函数
gives <- function(mx, lmx){
    mc <- mx[1]/lmx #lmx是mx的长度，mc是cos，ms是sin
    ms <- mx[2]/lmx
    matrix(c(mc,ms,-ms,mc),ncol=2)
}
x1 <- c(3,4)
lmx1 <- sqrt(sum(x1*x1))
lmx1
x1%*%gives(x1,lmx1)

5
5 -4.440892e-16

1
2
3
4
5
##原始的X矩阵和其Cholesky分解L
x=matrix(rnorm(60),10)
xtx=t(x)%*%x
L=mchol(xtx)
L

1.8958211 0.00000000 0.00000000 0.00000000 0.000000 0.000000
-1.2409637 3.60216335 0.00000000 0.00000000 0.000000 0.000000
-0.5486718 -0.44502336 2.87876658 0.00000000 0.000000 0.000000
-0.2459321 0.27392813 -0.52737167 2.81602139 0.000000 0.000000
-0.9724829 0.52179422 0.03686773 -0.05877291 1.500547 0.000000
0.3398855 -0.00289822 -0.23142886 0.80326680 -1.803605 1.795161

下面的例子显示了$L$在删掉第二行后不再是一个下三角矩阵,而是在对角线上多出来一些非零元素.要得到$X_{k}^T X_{k}$的Cholesky因子$L_{k}$,就是要通过Givens变换将这些对角线上的非零元一步步化为0.下面是以第一步化简为例.

1
2
3
4
5
##删掉了第二个变量后，对应的新的X的分解L只需在原来L的基础之上删去第二行，并使用Givens变换进行列变换化为下三角阵即可
##下面是以第一步化简为例
L2 =L[-2,] 
LL2 <- L2
L2 #注意！此处L2是一个5行6列的矩阵

1.8958211 0.00000000 0.00000000 0.00000000 0.000000 0.000000
-0.5486718 -0.44502336 2.87876658 0.00000000 0.000000 0.000000
-0.2459321 0.27392813 -0.52737167 2.81602139 0.000000 0.000000
-0.9724829 0.52179422 0.03686773 -0.05877291 1.500547 0.000000
0.3398855 -0.00289822 -0.23142886 0.80326680 -1.803605 1.795161

1
2
3
4
5
6
7
8
mx <- L2[2,2:3]  #需要将L2[2,3]这个对角线上的非零元化为0
lmx <- sqrt(sum(mx*mx))
L2[2,2:3] <- c(lmx,0) 
#givens变换将上面加粗的位置转化成了0
L2[3:5,2:3] <- L2[3:5,2:3] %*% gives(mx, lmx) #该向量下方的两列元素也被Given矩阵变换了，需要更新

L2%*%t(L2)
LL2%*%t(LL2)

3.8059636 3.7800191 0.8021618 0.3828310 0.9020923
3.7800191 9.6332301 -2.7118229 0.5274355 -3.6169142
0.8021618 -2.7118229 7.5498592 -2.0486273 1.2750979
0.3828310 0.5274355 -2.0486273 3.3299015 0.6356107
0.9020923 -3.6169142 1.2750979 0.6356107 21.5283570

3.8059636 3.7800191 0.8021618 0.3828310 0.9020923
3.7800191 9.6332301 -2.7118229 0.5274355 -3.6169142
0.8021618 -2.7118229 7.5498592 -2.0486273 1.2750979
0.3828310 0.5274355 -2.0486273 3.3299015 0.6356107
0.9020923 -3.6169142 1.2750979 0.6356107 21.5283570

氦核注：虽然矩阵内部数字改变，但是正交结果不会变化，我们成功将该位置元素转化为0，而且没有损失信息。

1
L[2:2,]

-0.00705475494350042 2.05665301740704 0 0 0 0

整个化为下三角的过程可以用一个while控制.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
p <- dim(x)[2] #x的列数
k <- 2  #被删掉的变量的下标
Lk <- L[-k,]
Lk #删掉后第k个变量后的系数矩阵L
mk <- k

##从第k列起，逐列将对角线以上的那个非零元约化为0
#关键句if：使进行givens变换后那两项下面的部分，进行循环更新
while( mk < p ){
    mx <- Lk[mk,mk:(mk+1)] #进行givens变换前的两项
    lmx <- sqrt(sum(mx*mx))
    Lk[mk,mk:(mk+1)] <- c(lmx,0) #进行givens变换后的两项
    if( mk < p-1 )Lk[(mk+1):(p-1), mk:(mk+1)] <- Lk[(mk+1):(p-1), mk:(mk+1)] %*% gives(mx, lmx) 
    #注：p-1是为了指示系数发生变化的位置有没有到达最后一列（因为涉及到mk+1列上元素）
    mk <- mk + 1
}
Lk
Lk <- Lk[,-p]
Lk

#与重新进行Cholesky分解的结果进行对比
xtxk <- xtx[-k, -k]
mchol(xtxk)

mgives函数

上面以一个例子展示了原理.下面的mgives函数就是实现这个功能的函数.这个函数在原有的$X_{k}^T X_{k}$的Cholesky分解的$L$基础之上,利用Givens变换计算删去第k个变量后新的$X_{k}^T X_{k}$的Cholesky分解$L_k$.它以上一步的$L_k$和这一次被删掉的变量的下标k为参数,返回$L_k$.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
mgives <- function(L,k){
    p <- dim(L)[1]
    if( k>p ) return ("Wrong input of k!")
    Lk <- L[-k,]
    mk <- k 
    while( mk < p ){
        mx <- Lk[mk,mk:(mk+1)]
        lmx <- sqrt(sum(mx*mx))
        Lk[mk,mk:(mk+1)] <- c(lmx,0)
        if( mk < p-1 ){
            Lk[(mk+1):(p-1), mk:(mk+1)] <- Lk[(mk+1):(p-1), mk:(mk+1)] %*% gives(mx, lmx)
        }
        mk <- mk + 1
    }
    return(Lk[,-p])
}

##查看每次删去一个变量后的Lk,并与重新进行Cholesky分解的结果比较
for (k in 1:dim(L)[1]){
    print("k=")
    print(k)
    print(mgives(L,k))
    xtxk <- xtx[-k, -k]
    print(mchol(xtxk)) #可以用paste写的更漂亮些。。。
}

[1] “k=”
[1] 1
[,1] [,2] [,3] [,4] [,5]
[1,] 2.05666512 0.0000000 0.0000000 0.0000000 0.00000
[2,] 0.15840137 3.4457975 0.0000000 0.0000000 0.00000
[3,] -0.74932067 -0.6720807 3.1096792 0.0000000 0.00000
[4,] -0.02895095 1.4368864 0.3519308 2.7493046 0.00000
[5,] 0.52731588 1.3429521 1.0401498 0.6622768 2.21406
[,1] [,2] [,3] [,4] [,5]
[1,] 2.05666512 0.0000000 0.0000000 0.0000000 0.00000
[2,] 0.15840137 3.4457975 0.0000000 0.0000000 0.00000
[3,] -0.74932067 -0.6720807 3.1096792 0.0000000 0.00000
[4,] -0.02895095 1.4368864 0.3519308 2.7493046 0.00000
[5,] 0.52731588 1.3429521 1.0401498 0.6622768 2.21406
[1] “k=”
[1] 2
[,1] [,2] [,3] [,4] [,5]
[1,] 2.4663786 0.0000000 0.0000000 0.0000000 0.000000
[2,] 1.0463703 3.2869014 0.0000000 0.0000000 0.000000
[3,] -1.4298928 -0.2854816 2.9252667 0.0000000 0.000000
[4,] 0.7702886 1.2597355 0.5508706 2.6952319 0.000000
[5,] 0.2755799 1.3455560 0.9281248 0.6243126 2.316575
[,1] [,2] [,3] [,4] [,5]
[1,] 2.4663786 0.0000000 0.0000000 0.0000000 0.000000
[2,] 1.0463703 3.2869014 0.0000000 0.0000000 0.000000
[3,] -1.4298928 -0.2854816 2.9252667 0.0000000 0.000000
[4,] 0.7702886 1.2597355 0.5508706 2.6952319 0.000000
[5,] 0.2755799 1.3455560 0.9281248 0.6243126 2.316575
[1] “k=”
[1] 3
[,1] [,2] [,3] [,4] [,5]
[1,] 2.466378596 0.00000000 0.0000000 0.000000 0.000000
[2,] -0.007054755 2.05665302 0.0000000 0.000000 0.000000
[3,] -1.429892844 -0.75422991 2.8407433 0.000000 0.000000
[4,] 0.770288640 -0.02630887 0.4336785 2.994311 0.000000
[5,] 0.275579906 0.52826428 0.9607745 1.164281 2.422916
[,1] [,2] [,3] [,4] [,5]
[1,] 2.466378596 0.00000000 0.0000000 0.000000 0.000000
[2,] -0.007054755 2.05665302 0.0000000 0.000000 0.000000
[3,] -1.429892844 -0.75422991 2.8407433 0.000000 0.000000
[4,] 0.770288640 -0.02630887 0.4336785 2.994311 0.000000
[5,] 0.275579906 0.52826428 0.9607745 1.164281 2.422916
[1] “k=”
[1] 4
[,1] [,2] [,3] [,4] [,5]
[1,] 2.466378596 0.00000000 0.000000 0.0000000 0.000000
[2,] -0.007054755 2.05665302 0.000000 0.0000000 0.000000
[3,] 1.046370303 0.16199157 3.282907 0.0000000 0.000000
[4,] 0.770288640 -0.02630887 1.262566 2.7495274 0.000000
[5,] 0.275579906 0.52826428 1.321126 0.8128221 2.396478
[,1] [,2] [,3] [,4] [,5]
[1,] 2.466378596 0.00000000 0.000000 0.0000000 0.000000
[2,] -0.007054755 2.05665302 0.000000 0.0000000 0.000000
[3,] 1.046370303 0.16199157 3.282907 0.0000000 0.000000
[4,] 0.770288640 -0.02630887 1.262566 2.7495274 0.000000
[5,] 0.275579906 0.52826428 1.321126 0.8128221 2.396478
[1] “k=”
[1] 5
[,1] [,2] [,3] [,4] [,5]
[1,] 2.466378596 0.0000000 0.0000000 0.000000 0.000000
[2,] -0.007054755 2.0566530 0.0000000 0.000000 0.000000
[3,] 1.046370303 0.1619916 3.2829071 0.000000 0.000000
[4,] -1.429892844 -0.7542299 -0.2486123 2.829844 0.000000
[5,] 0.275579906 0.5282643 1.3211265 1.080541 2.288278
[,1] [,2] [,3] [,4] [,5]
[1,] 2.466378596 0.0000000 0.0000000 0.000000 0.000000
[2,] -0.007054755 2.0566530 0.0000000 0.000000 0.000000
[3,] 1.046370303 0.1619916 3.2829071 0.000000 0.000000
[4,] -1.429892844 -0.7542299 -0.2486123 2.829844 0.000000
[5,] 0.275579906 0.5282643 1.3211265 1.080541 2.288278
[1] “k=”
[1] 6
[,1] [,2] [,3] [,4] [,5]
[1,] 2.466378596 0.00000000 0.0000000 0.000000 0.000000
[2,] -0.007054755 2.05665302 0.0000000 0.000000 0.000000
[3,] 1.046370303 0.16199157 3.2829071 0.000000 0.000000
[4,] -1.429892844 -0.75422991 -0.2486123 2.829844 0.000000
[5,] 0.770288640 -0.02630887 1.2625664 0.546270 2.694715
[,1] [,2] [,3] [,4] [,5]
[1,] 2.466378596 0.00000000 0.0000000 0.000000 0.000000
[2,] -0.007054755 2.05665302 0.0000000 0.000000 0.000000
[3,] 1.046370303 0.16199157 3.2829071 0.000000 0.000000
[4,] -1.429892844 -0.75422991 -0.2486123 2.829844 0.000000
[5,] 0.770288640 -0.02630887 1.2625664 0.546270 2.694715

增加一个变量:$L_k$的计算

forupdate函数

上面一部分讲了减少一个变量后求解$L_k$的简化计算方法。那在增加一个变量的时候，同样可以通过矩阵运算简化计算。当在$X$的末尾再加一列$x_k$时,新的$X_{k}^T X_{k}$的分解$L_k$是在原来$L$的基础之上在底下多加了一行,多加的元素的计算可以通过矩阵运算很快的得到.下面的这个forupdate函数就是完成在$L$的基础之上计算$L_k$的功能的.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
forupdate <- function(L, xxk, xkxk){
    lk <- mforwardsolve(L, xxk) #xxk是新的变量在变量列表里的位置（序号）
    lkk <- sqrt(xkxk - sum(lk*lk)) #lkk是化简用的矩阵，运算中间步骤
    return( as.matrix( rbind( cbind(L,0),c(lk,lkk) ) ) )
}

x <- matrix(rnorm(60),10)
xtx <- t(x)%*%x


#在3,2,4个变量的基础之上，加入第5个变量
A <- c(3,2,4)  #注意：A里面带有变量的顺序，因为有些情形需要对变量的重要性等进行目标规划
L <- mchol(xtx[A,A])  #三个变量的xtx的分解
k <- 5 #加入第五个变量
xxk <- xtx[A,k,drop=T] 
xkxk <- xtx[k,k] #计算扩充L时所需要的数
forupdate (L, xxk, xkxk)

A <- c(A, k)
mchol(xtx[A,A]) #和重新分解对比

2.18726152 0.0000000 0.000000 0.000000
-0.31876255 2.5768506 0.000000 0.000000
-0.09862756 -1.7100174 2.244832 0.000000
0.17257491 0.2771965 1.469177 2.626423

2.18726152 0.0000000 0.000000 0.000000
-0.31876255 2.5768506 0.000000 0.000000
-0.09862756 -1.7100174 2.244832 0.000000
0.17257491 0.2771965 1.469177 2.626423

逐步回归的实现(实际操作)

很长的铺垫过后，我们已经得到逐步回归增减变量时简化计算的工具,下面就可以开始实现逐步回归算法了.首先导入数据并查看使用R自带的函数的回归结果.

1
2
3
4
library(ElemStatLearn)
data(prostate)
data <- prostate[,-10]
head(data)

>
lcavol lweight age lbph svi lcp gleason pgg45 lpsa
-0.5798185 2.769459 50 -1.386294 0 -1.386294 6 0 -0.4307829
-0.9942523 3.319626 58 -1.386294 0 -1.386294 6 0 -0.1625189
-0.5108256 2.691243 74 -1.386294 0 -1.386294 7 20 -0.1625189
-1.2039728 3.282789 58 -1.386294 0 -1.386294 6 0 -0.1625189
0.7514161 3.432373 62 -1.386294 0 -1.386294 6 0 0.3715636
-1.0498221 3.228826 50 -1.386294 0 -1.386294 6 0 0.7654678

1
2
3
4
lm1= lm(lpsa~.,data)
summary(lm1)
extractAIC(lm1) #  nln(RSS/n)+ 2*9, not AIC
step(lm1, direction="both")
1
2
3
summary(lm1)
extractAIC(lm1)
AIC(lm1)

结果没贴，大家粘贴一运行便知。

下面实现我们写的逐步回归。

首先是把所有变量都加入模型进行回归,求解系数,并计算该模型的RSS和AIC。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
np <- dim(data)
n <- np[1]
p <- np[2]-1
xn <- names(data)[1:p]
x <- as.matrix( cbind(1,data[,1:p]) ) #将1与矩阵并联起来，将生成一列1，这是代码高手的简化写法（划重点）
y <- data[,p+1]


xtx <- t(x)%*%x #x平方阵
xty <- drop(t(x)%*%y) 
yty <- sum(y*y) #y平方阵


L <- mchol(xtx) #进行Cholesky分解，下三角矩阵L
tb <- mforwardsolve(L, xty) 
b <- mbacksolve(t(L), tb)
b

RSS <- yty - sum(tb*tb)  ###神奇操作，如何理解？
AICF <- n*log(RSS/n) + 2*(p+1)
AICF

#A1 <- rep(TRUE, p)
#A <- c(TRUE, A1)

A <- 1:(p+1)
LA <- L
MAIC <- AICF
mAIC <- AICF
flag <- rep(FALSE,p)

对于下段代码中的

RSS <- yty - sum(tb*tb)

有如下解释：

mforwardsolve(…)的结果为$L^T \cdot \beta$，其中$L$是$X^T X$经过Cholesky分解后的下三角阵，$L L^T=X^T X$。
$ \begin{equation}
RSS=(Y-X\beta)^T(Y-X\beta)\\
=Y^T Y-\beta^T X^T Y-Y^T X\beta+\beta^T X^T X \beta
\end{equation} $
因为$\beta^T X^T Y$是一个数，相乘的项一样，因此所得结果也必然一样，所以$\beta^T X^T Y=Y^T X\beta$，所以
$ \begin{equation}
RSS=Y^T Y-\beta^T X^T Y\\
=Y^T Y-\beta^T X^T X \beta\\
=Y^T Y-\beta^T L L^T \beta\\
=Y^T Y-(L^T \beta)^T(L^T \beta)
\end{equation} $

整个逐步回归的过程由一个repeat完成.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
repeat{
    
    #print(flag)
    if( length(A) < p+1 ){
        #B <- setdiff( 1:(p+1), A)
        B <- (2:(p+1))[flag]
        #flag[B-1] <- TRUE
        #print(flag)
    }else{
        B <- NULL
    }
    AB <- c(A,B)
    bm <- matrix(0,p,p+1)
    AICm <- rep(0,p)

    #A <- c(TRUE, A1)
    #a1 <- (1:(p+1))[A]
    #c1 <- (1:(p+1))[!A]
    #tp <- length(a1)
    ff <- 2
    #print(AB)
    for(k in AB[2:(p+1)]){
        if( !flag[k-1]){
            
            Lk <- mgives(LA,ff)
            tA <- A[-ff]
            xtyk <- xty[ tA ]
            #print(Lk)
            #print(xtyk)
            tbk <- mforwardsolve(Lk, xtyk)
            bk <- mbacksolve(t(Lk), tbk)
            #tA <- c(1, tA)
            #print(bk)
            #print(tA)
            bm[k-1,tA] <- bk
            #print(bk)
            RSSk <- yty - sum(tbk*tbk)  
            AICk <- n*log(RSSk/n) + 2*( length(A)-1)
            AICm[k-1] <- AICk
            ff <- ff+1
            if( AICk < mAIC ){
                mink <- k
                mtA <- tA
                mAIC <- AICk
                mLA <- Lk
                hb <- bm[k-1,]
            }
            
        } else {
            #print (k)
            xxk <- xtx[A,k,drop=T] 
            xkxk <- xtx[k,k]
            Lk <- forupdate (LA, xxk, xkxk)
            tA <- c(A,k)
            xtyk <- xty[tA]
            tbk <- mforwardsolve(Lk, xtyk)
            bk <- mbacksolve(t(Lk), tbk)
            
            
            bm[k-1,tA] <- bk
            RSSk <- yty - sum(tbk*tbk)   ###神奇操作
            AICk <- n*log(RSSk/n) + 2*(length(A)+1)
            AICm[k-1] <- AICk
            if( AICk < mAIC ){
                mink <- k
                mtA <- tA
                mAIC <- AICk
                mLA <- Lk
                hb <- bm[k-1,]
            }            
        }
    }
    #print(bm)
    #print(AICm)
    #print(flag)
    
    if(mAIC >= MAIC) break
    if (mAIC<MAIC ){
        #print(flag)
        flag[mink-1] = !flag[mink-1]
        A <- mtA
        MAIC <- mAIC
        print(MAIC)
        hhb <- hb
        LA <- mLA
    }

}
re <- data.frame(matrix(hhb[c(TRUE,!flag)],nrow=1))
names(re) <- c("inter", xn[!flag])
re

#}

将上述逐步回归的过程函数化,用一个steplm()来完成.它以样本数据(最后一列为因变量)为参数,返回最优的模型的结果.

下面以上面加载的数据data为例,使用了这个函数.并且例子显示,不论变量顺序如何都不影响最终结果.

最终函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
steplm <- function(data){
    np <- dim(data)
    n <- np[1]
    p <- np[2]-1
    xn <- c("inter", names(data)[1:p])
    x <- as.matrix( cbind(1,data[,1:p]) ) 
    y <- data[,p+1]


    xtx <- t(x)%*%x #x的平方阵
    xty <- drop(t(x)%*%y)
    yty <- sum(y*y) #y的平方阵


    L <- mchol(xtx)
    tb <- mforwardsolve(L, xty)
    b <- mbacksolve(t(L), tb)
    #b

    RSS <- yty - sum(tb*tb) 
    AICF <- n*log(RSS/n) + 2*(p+1)
    #AICF

#A1 <- rep(TRUE, p)
#A <- c(TRUE, A1)

    A <- 1:(p+1)
    LA <- L
    MAIC <- AICF
    mAIC <- AICF
    MFLAG <- c(TRUE, rep(FALSE, p))
    flag <- rep(TRUE,p+1)
    hbb <- b
    repeat{
        if( length(A) < p+1 ){
             B <- (1:(p+1))[!flag]
        }
        else{
            B <- NULL
        }
        AB <- c(A,B)
        bm <- matrix(0,p,p+1)
        AICm <- rep(0,p)

        ff <- 1
        #print(AB)
        for(k in AB){
            if(MFLAG[k]){
                ff <- ff+1
            } else {
                if(flag[k]){
                    Lk <- mgives(LA,ff)
                    tA <- A[-ff]
                    xtyk <- xty[ tA ]
                    ff <- ff+1
                } else {
                    xxk <- xtx[A,k,drop=T] 
                    xkxk <- xtx[k,k]
                    Lk <- forupdate (LA, xxk, xkxk)
                    tA <- c(A,k)
                    xtyk <- xty[tA]
                }
                #print(A)
                #print(B)
                #print(k)
                tbk <- mforwardsolve(Lk, xtyk)
                bk <- mbacksolve(t(Lk), tbk)
            #tA <- c(1, tA)
                #print(tbk)
            #print(tA)
                bm[k-1,tA] <- bk
            #print(bk)
                RSSk <- yty - sum(tbk*tbk)  
                AICk <- n*log(RSSk/n) + 2*length(tA)
                AICm[k-1] <- AICk
       
                if( AICk < mAIC ){
                    mink <- k
                    mtA <- tA
                    mAIC <- AICk
                    mLA <- Lk
                    hb <- bm[k-1,]
                }
            }
        }

        if(mAIC >= MAIC) break
        if ( mAIC<MAIC ){
            flag[mink] = !flag[mink]
            A <- mtA
            MAIC <- mAIC
            hhb <- hb
            LA <- mLA
        }

    }
    re <- data.frame(matrix(hhb[c(flag)],nrow=1))
    names(re) <- xn[flag]
    return(re)
}
steplm(data)
data <- data[,c(sample(1:8,8),9)]
steplm(data)

下面是R自带的逐步回归的函数.

1
coef(step(lm1, direction="both"))

注：理论部分参考：徐树方, 高立, 张平文. 数值线性代数.第2版[M]. 北京大学出版社, 2013.


本文链接： https://konelane.github.io/2019/03/22/190322逐步回归-Cholesky分解法/

-- EOF --

转载请注明出处 署名-非商业性使用-禁止演绎 3.0 国际（CC BY-NC-ND 3.0）

￥^￥请氦核牛饮一盒奶~suki