线性回归

一、问题描述

线性回归问题非常简单。对于给定的m个样本$\left\{ \left( {{x}^{\left( 1 \right)}},{{y}^{\left( 1 \right)}} \right),\left( {{x}^{\left( 2 \right)}},{{y}^{\left( 2 \right)}} \right),...,\left( {{x}^{\left( i \right)}},{{y}^{\left( i \right)}} \right),...,\left( {{x}^{\left( m \right)}},{{y}^{\left( m \right)}} \right) \right\}$ ，其中${{x}^{\left( i \right)}}\in {{\mathbb{R}}^{n}}$为第 i个样本的自变量，表示成n维列向量，即${{x}^{\left( i \right)}}={{\left[ \begin{matrix}x_{1}^{\left( i \right)}, & x_{2}^{\left( i \right)}, & x_{3}^{\left( i \right)}, & \cdots  & ,x_{n}^{\left( i \right)}  \\ \end{matrix} \right]}^{T}}$ (这里为了方便已经对 ${{x}^{\left( i \right)}}$ 进行了增广，即对n-1维的变量增加取值为1的一维，组成n维，这样可以方便地将公式写成向量形式），${{y}^{\left( i \right)}}\in {{\mathbb{R}}^{1}}$ 为对应的函数值。

线性回归的目标是找到自变量和函数值之间的一个最佳线性拟合${{h}_{\theta }}\left( x \right)={{\theta }^{T}}x$，其中 $\theta \in {{\mathbb{R}}^{n}}$ ，表示为n维列向量，即$\theta ={{\left[ \begin{matrix}{{\theta }_{1}}, & {{\theta }_{2}}, & {{\theta }_{3}}, & \cdots  & ,{{\theta }_{n}}  \\ \end{matrix} \right]}^{T}}$，我们要找到一个$\theta $使得这种拟合是最优的。

二、求解

为了找到最优的$\theta $，需要首先定义一个目标函数来说明什么样的$\theta $是最优的。线性回归采用的是最小平方误差准则。

定义目标函数$\ell \left( \theta  \right)=\frac{1}{2}\sum\limits_{i=1}^{m}{{{\left( {{h}_{\theta }}\left( {{x}^{\left( i \right)}} \right)-{{y}^{\left( i \right)}} \right)}^{2}}}=\frac{1}{2}\sum\limits_{i=1}^{m}{{{\left( {{\theta }^{T}}{{x}^{\left( i \right)}}-{{y}^{\left( i \right)}} \right)}^{2}}}$，显然，求和式内部的每一项为预测值和原始值的平方误差差，所以目标函数定义了所以样本的预测值和原始值平方差的和，要使得线性拟合最优，那么就需要选取使得这个函数取值最小。下面来求解$\theta $。这是一个优化问题，可以用不同的优化方法来做，比如梯度下降或者牛顿法啊什么的，当然能直接得到解析解最好了~对于现在这个问题，我们可以很方便求出它的解析解。

为了方便，我们将目标函数写成向量形式。

令

$\begin{align} X=\left[ \begin{matrix}{{x}^{\left( 1 \right)T}}  \\{{x}^{\left( 2 \right)T}}  \\{{x}^{\left( 3 \right)T}}  \\ \vdots   \\{{x}^{\left( m \right)T}}  \\ \end{matrix} \right] y=\left[ \begin{matrix} {{y}^{\left( 1 \right)}}  \\{{y}^{\left( 2 \right)}}  \\{{y}^{\left( 3 \right)}}  \\\vdots   \\{{y}^{\left( m \right)}}  \\ \end{matrix} \right] \theta =\left[ \begin{matrix}{{\theta }_{1}}  \\{{\theta }_{2}}  \\{{\theta }_{3}}  \\ \vdots   \\{{\theta }_{n}}  \\\end{matrix} \right] \end{align}$

这样目标函数

$\begin{align}\ell (\theta )&=\frac{1}{2}{{(X\theta -y)}^{T}}(X\theta -y) \\& =\frac{1}{2}({{\theta }^{T}}{{X}^{T}}-{{y}^{T}})(X\theta -y) \\& =\frac{1}{2}({{\theta }^{T}}{{X}^{T}}X\theta -2{{y}^{T}}X\theta -{{y}^{T}}y) \end{align}$

为了求目标函数极值，对其求导：\[\frac{d\ell (\theta )}{d\theta }=\frac{d\frac{1}{2}({{\theta }^{T}}{{X}^{T}}X\theta -2{{y}^{T}}X\theta +{{y}^{T}}y)}{d\theta }={{X}^{T}}X\theta -{{X}^{T}}y\]

（上式利用了求导公式$\frac{d({{x}^{T}}Ax)}{dx}=2Ax$和$\frac{d({{b}^{T}}x)}{dx}=b$，令$A={{X}^{T}}X$，${{b}^{T}}={{y}^{T}}X$）

令倒数为0：\[\begin{align}&&{{X}^{T}}X\theta -{{X}^{T}}y&=0 \\ &\Rightarrow &{{X}^{T}}X\theta &={{X}^{T}}y \\ &\Rightarrow &\theta &={{({{X}^{T}}X)}^{-1}}{{X}^{T}}y \\\end{align}\]

这样便求出了最优的参数$\theta $，回归曲线${{h}_{\theta }}(x)$便也求出了。

回过头来再看，上面的推导有两点问题：

（1）对于目标函数$\ell \left( \theta  \right)=\frac{1}{2}({{\theta }^{T}}{{X}^{T}}X\theta -2{{y}^{T}}X\theta +{{y}^{T}}y)$，我们求它的最小值，那么最小值存在吗？感性的认识一下，由于我们找的拟合直线（或超平面）是距离所有样本偏差平方和最小的，那么目标函数最新值应该是存在的，直观上好像总能找到这样一个直线（或超平面），如果要找最大值，那当然不存在了，我们可以将直线（或超平面）放置的离样本无穷远。但是现在我们这里求解的是解析解，所以还是需要一下说明的。我们知道，如果一个函数的Hessien矩阵正定或者半正定，那么这个函数是一个凸函数（开口向上），所以是梯度为零的驻点就是全局函数值最小点。显然，现在的目标函数是关于的二次函数，其Hessien矩阵为${{X}^{T}}X$，那么它半正定吗？其实我也不知道。。。

查书找到一个结论，如果为$X$方阵，那么${{X}^{T}}X$是半正定的，可惜这里$X$不是方阵，网上问了一下又说${{X}^{T}}X$是半正定的，但不知道怎么证明，所以也不敢确定喽~

（2）对于${{X}^{T}}X\theta ={{X}^{T}}y\Rightarrow \theta ={{({{X}^{T}}X)}^{-1}}{{X}^{T}}y$这一步推导 ，需要${{X}^{T}}X$可逆，那么它可逆吗？如果问题（1）中${{X}^{T}}X$正定，那么正定矩阵是可逆的，就不存在这个问题了。但如果如（1）中网友所言${{X}^{T}}X$是半正定，那么能说明它可逆吗？这个问题我也不知道怎么弄。

在实际中，如果${{X}^{T}}X$不可逆，为了得到一个解，我们可以用它的伪逆当作逆来用，这样做在数学上是怎么回事我也就不知道了，人家数学家说可以这样弄个，我很天真，就相信了。。。

三、概率模型解释

下面从概率模型角度来解释一下上面线性回归的目标函数为何是最小二乘。

对于一个样本的自变量${{x}^{\left( i \right)}}$，线性回归对其函数值预测值为${{h}_{\theta }}({{x}^{(i)}})$，而它的真实值为${{y}^{\left( i \right)}}$，两者之间是有一定偏差的。假设这个偏差$\xi$本身是服从均值为0，方差为${{\sigma }^{2}}$的正态分布的，即$\xi \sim N(0,{{\sigma }^{2}})$。

那么由于${{y}^{\left( i \right)}}={{h}_{\theta }}({{x}^{\left( i \right)}})+\xi$，所以${{y}^{(i)}}\sim N({{h}_{\theta }}({{x}^{(i)}}),{{\sigma }^{2}})$，这样，根据正态分布的概率密度函数，可以得出$p({{y}^{\left( i \right)}}|{{x}^{\left( i \right)}};\theta )=\frac{1}{\sqrt{2\pi }\sigma }\exp (-\frac{{{({{h}_{\theta }}({{x}^{(i)}})-{{y}^{(i)}})}^{2}}}{2{{\sigma }^{2}}})=\frac{1}{\sqrt{2\pi }\sigma }\exp (-\frac{{{({{\theta }^{T}}{{x}^{(i)}}-{{y}^{(i)}})}^{2}}}{2{{\sigma }^{2}}})$

现在，可以写出对参数$\theta$估计的对数似然函数了

$\begin{align}L(\theta )&=\log \prod\limits_{i=1}^{m}{p({{y}^{(i)}}|{{x}^{(i)}};\theta )} \\& =\sum\limits_{i=1}^{m}{\log (\frac{1}{\sqrt{2\pi }\sigma }\exp (-\frac{{{({{\theta }^{T}}{{x}^{(i)}}-{{y}^{(i)}})}^{2}}}{2{{\sigma }^{2}}})}) \\& =\sum\limits_{i=1}^{m}{(\log \frac{1}{\sqrt{2\pi }\sigma }-\frac{{{({{\theta }^{T}}{{x}^{(i)}}-{{y}^{(i)}})}^{2}}}{2{{\sigma }^{2}}}}) \end{align}$

为了求解最优参数$\theta$，需要对似然函数最大化，而对上面似然函数的最大化就等价与对$\sum\limits_{i=1}^{m}{{{\left( {{\theta }^{T}}{{x}^{\left( i \right)}}-{{y}^{\left( i \right)}} \right)}^{2}}}$最小化，这和我们开始的最小二次的目标是一致的，这就从概率模型角度解释了线性回归采用最小二次作为目标的合理性。

好吧，我理解的线性回归就这些了~~

四、简单实现

一维的吧(${{x}^{\left( i \right)}}$增广后为二维），多维的结果不好展示~

n <- 100; #样本数目
x <- runif(n, min = 1, max = 10); #随机生成均匀分布数据
x <- matrix(x,nrow=n); #为了计算搞成矩阵类型
y <- 3*x+4+rnorm(n,0,2); #生成y，并加入高斯噪声
y <- matrix(y,nrow=n);
plot(x,y);

#R自带的线性回归函数
cof <- lm(y~x)$coefficients;
print("cof")
print(cof);
abline(cof[1],cof[2],col="blue");

#套公式自己实现试试
x <- cbind(rep(1,n),x); #对x增广
theta <- solve(t(x)%*%x) %*% t(x) %*% y;
print("theta");
print(theta);
abline(theta[1],theta[2],col="red");

从下面结果可以看出，R自带的lm函数和自己计算的结果是一样的，但是和我们设置的解[4,3]还差一丢丢，样本和噪声的问题吧~~

 

就这样吧~

参考资料：

[1]Andrew Ng机器学习视频和课件：http://cs229.stanford.edu/