多元线性回归

多元线性回归

假设我们的输入变量为$X=(X_1,X_2,...,X_p)$，输出变量为$Y$。$X\in R^{n\times p}$，代表$n$个样本 $p$个变量。这里线性假设可以写为：

$$f(X_i)=\beta_0+\beta_1X_{1i}+\beta_2X_{2i}\ +\ ...+\ \beta_pX_{pi}$$
写成矩阵形式$f(X)=X\beta$，其中$\beta=(\beta_0,\beta_1,...,\beta_p)^T$，多元线性回归就是输入变量个数多于一个的回归问题。

为了解决这个多元线性回归问题，首先定义下相应的损失函数：

$$J(\beta_0,\beta_1,...,\beta_p) =\frac{1}{2n}\sum_{i=1}^{n}(f(X_i)-Y_i)^2$$
这里$X_i$表示一个输入样本，$Y_i$表示相应的输出值。下面我们用$X_{ji}$表示第$i$个观测样本的第$j$个特征值
多元线性回归梯度下降算法伪代码可以写为：

New Algorithm $(\ p\ge 1\ )$：
Repeat until convergence{
$\ \ \ \ \beta_j:=\beta_j-\alpha \frac{1}{n}\sum_{i=1}^{n}(f(X_i)-Y_i)X_{ji}$
$\ \ \ \ $(simultaneous update $\beta_j$ for $j=0,1,..,p$)
}

在求解多元线性回归时往往还需要对数据进行标准化处理，数据标准化能够消除数据不同量纲带来的影响，也有助于优化程序快速找到最优值。
常用的数据标准化方法有：

• Min-max 标准化
  新数据=（原数据-极小值）/（极大值-极小值），Min-max 标准化将原数据映射到$(0,1)$之间。

• z-score 标准化
  新数据=（原数据-均值）/标准差，z-score 标准化得到的新数据服从标准正态分布。

除了对多个变量做回归，我们还可以把线性回归模型推广到多项式回归模型，比如

$$f(x)=\beta_0+\beta_1x+\beta_2x^2$$
下面我们就用梯度下降算法来解决一个多元回归例子。
（数据来自R语言 MASS程序包，Boston数据集。）

#linear regression by gradient decent algorithm
#auther: JiZhG
#
#gradient decent algorithm for linear regression
#x: input
#y: output
#alpha: learning rate default alpha=0.05
#iter: max iterations default alpha=500
#tol default 1e-10
LinearRegression=function(x,y,alpha=0.05,iter=1e+5, tol=1e-5){
  n=length(y)
  x=cbind(rep(1,n),x)
  p=dim(x)[2]
  beta=rep(0,p)
  temp=rep(0,p)
  i=1
  while(i<=iter){
    beta_old=beta
    diff=x%*%beta-y
    for(j in 1:p){
      temp[j]=1/n*sum(diff*x[,j])
    }
    beta=beta-alpha*temp
    #print(beta)
    if(sum(abs(beta-beta_old))<tol) break
    i=i+1
  }
  beta
}
#
library(MASS)
attach(Boston)
#second-order polynomial
input=cbind(lstat,lstat^2)
#normalization
input=scale(input)
plot(input[,1],medv,pch=3,lwd=2,xlab="lstat")
beta=LinearRegression(input,medv,alpha=0.2)
#beta
#22.53281 -16.65819  10.27913
#prdecition
pred=cbind(rep(1,dim(input)[1]),input)%*%beta
#lines
lines(input[,1][order(input[,1])],pred[order(input[,1])],lwd=2,col="blue")
#third-order polynomial
input=cbind(lstat,lstat^2,lstat^3)
#normalization
input=scale(input)
beta=LinearRegression(input,medv,alpha=0.5)
pred=cbind(rep(1,dim(input)[1]),input)%*%beta
#lines
lines(input[,1][order(input[,1])],pred[order(input[,1])],lwd=2,col="red")
legend("topright",legend=c("second-order","third-order"),col=c("blue","red"),lty=c(1,1),lwd=2,bty="n")

以上代码根据梯度下降算法分别解决了2阶和3阶多项式拟合问题，拟合曲线如图，蓝色代表2阶多项式拟合，红色代表3阶多项式拟合。值得一提的是学习率$\alpha$的选取很可能会决定程序的收敛与否，而且在用梯度下降算法求解优化问题时对数据进行标准化处理非常重要！

以上我们是通过梯度下降算法解决多元回归问题，梯度下降算法通过不停迭代更新参数来找到最优解，下面介绍一种解决多元线性回归更为高效的方法——正规方程组法，正规方程组法是一种解析解法。
为了最小化

$$J(\beta_0,\beta_1,...,\beta_p) =\frac{1}{2n}\sum_{i=1}^{n}(f(X_i)-Y_i)^2$$
我们可以令各个偏导
$$\frac{\partial}{\partial{\beta_j}}=\dots=0,\ \ (for \ every \ j)$$
然后解这个线性方程组，得到最优解$\beta_0,\beta_1,...,\beta_p$,写成矩阵形式有
$$\hat{\beta}=(X^TX)^{-1}X^TY$$
下面根据这个公式来解上面的2阶多项式拟合。

#normalization
input=scale(input)
input=cbind(rep(1,dim(input)[1]),input)
beta=solve(t(input)%*%input)%*%t(input)%*%medv
#beta
#22.53281 -16.65882 10.27976

这里得到的参数值和之前梯度下降算法得到的参数值基本一致。

至于应该使用哪个算法，我们先看一下这两个算法的优缺点。

梯度下降算法	正规方程组法
需要选择学习率$\alpha$	不需要选择学习率
需要多次迭代	不需要迭代
当$n$很大时也能很好的工作	需要计算$(X^TX)^{-1}$,当$n$很大时速度比较慢

总的来说，当样本量不大时，我们可以毫不犹豫的选择正规方程组法，正规方程组法需要计算矩阵逆，矩阵逆的计算复杂度约为$O(n^3)$，当很大甚至上百万时，梯度下降算法会是一个不错的选择。