线性回归

线性回归

用\(X=(x_1, x_2, ...,x_n)^T \in \mathbb{R}^{n\times p}\)表示数据，用\(y=(y_1, y_2, ...,y_n)^T \in \mathbb{R}^{n}\)表示标签。对于一个样本\(x_i\)，它的输出值是其特征的线性组合：

\[\begin{equation} f(x_i) = \sum_{m=1}^{p}w_m x_{im}+w_0={w}^T{x_i} \end{equation} \]

线性回归的目标是用预测结果尽可能地拟合目标label，用最常见的Least square作为loss function：

\[\begin{equation} J(w)=\frac{1}{n}\sum_{i=1}^{n}(y_i - f(x_i))^2=\frac{1}{n}\|y-Xw\|^2 \end{equation} \]

性能指标

回归模型最简单常用的两种指标是平方根误差(Root Mean Square Error,RMSE)和R平方值(\(R^2\)).

  • RMSE 是 MSE的平方根: \(\sqrt{\frac{1}{n}\sum_{i=1}^{n}(y_i - f(x_i))^2}\)
    缺点: 数值大小与数据集相关, 不容易在不同数据集上进行比较. 如果预测值和真实值都比较大,那么RMSE也较大.
  • R平方: 是一个0~1之间的相对值,克服了RMSE的缺点. R平方越接近于1,预测效果越好.

\[R^2 = 1 - \frac{\sum_{i=1}^{n}(y_i - f(x_i))^2}{\sum_{i=1}^{n}(y_i - \bar y))^2} \]

局部加权线性回归

Locally Weighted Lieanear Regression (LWR) 是一种非参数学习方法,不适合于数据量比较大的情况,每次数据改变时均需要重新计算.

LWR: Fit \(\theta\) to minimize

\[\sum_i \omega^{(i)}(y^{(i)}-\theta^T x^{(i)})^2 \]

where \(\omega^{(i)}=\exp(-\frac{(x^{(i)}-x)^2}{2})\)

当\(x^{(i)}\)和x非常接近时,\(\omega^{(i)}\)接近于1;相差很大时,接近于0.
这样离 x 越近的样本权重越大,越远的影响越小。这个公式与高斯分布类似,但有所不同,因为 w (i) 不是随机变量。
此方法称为非参数学习算法,因为误差函数loss随着预测值的不同而不同,这样 θ 无法事先确定,预测一次需要临时计算.
一个更加一般化的公式是在\(\omega^{(i)}\)的公式中分母2的后边加上\(\tau^2\),称为bandwith parameter(波长函数),它控制了权值随距离下降的速率.
\(\omega^{(i)}=\exp(-\frac{(x^{(i)}-x)^2}{2\tau^2})\)

单元多项式回归

上边提到的局部加权线性回归属于多元线性回归,除此之外还有类线性的单元多项式回归. 如多项式假设:\(h_\theta=\theta_0+\theta_1x_1+\theta_2x_1^2+\theta_3x_1^3+\cdots+\theta_mx_1^m\)
单元多项式回归可以轻易地转化为线性回归:将高次项用新变量替换即可,转化为线性假设函数\(h_\theta=\theta_0+\theta_1x_1+\theta_2x_2+\theta_3x_3+\cdots+\theta_mx_m\).

线性回归求解思路

在求一个点的预测值时,在该点附近一段区间内,拟合一条直线,从而求得预测值.

多元线性回归的求解方法可以是梯度下降和Normal Equation,以及二阶的牛顿法等。最小二乘法可以求出全局最优，但是大规模矩阵的求逆运算复杂度较高，因此在数据量较大时常使用梯度下降算法。牛顿法的收敛速度较快,但同样求逆运算复杂度较高.

一元线性回归-最小二乘

当我们希望预测值和真实值的平方和尽量小时，损失函数为\(\sum_{i=1}^n(y_i-\hat y_i)^2\), 这个就是最小二乘线性回归(Ordinary Least Squares, OLS)。普通最小二乘法是线性回归预测问题中一个很重要的概念， 并且十分有用，例如可以用来做推荐系统、资金流动预测等。

当我们希望预测值和真实值的绝对值的和尽量小时，损失函数为\(\sum_{i=1}^n|y_i-\hat y_i|\), 这个就是最小绝对偏差(Least Absolute Deviation, LAD)线性回归。还有其他的不同的更复杂的损失函数，造就了不同的线性回归模型，比如分位数回归，Huber回归，等等。

由于最常用的、最简单的线性回归模型就是最小二乘法，所以通常说到线性回归的时候其实就是指的最小二乘线性回归。下面推导下我们中学数学课本上所学的最小二乘公式。

公式推导

对观察点 \((x_i,y_i),i=1,2,...,n\) 拟合一条直线: \(y=a+bx\),

均方误差为:\(D=\sum_{i=1}^n(y_i-\hat y_i)^2=\sum_{i=1}^n(y_i-a-b x_i)^2\),对a,b分别求偏导,并令导数为0,得:

\[\begin{align} {\partial D\over \partial a} &=-2\big (\sum_i y_i-na-b\sum_i x_i\big )=0 \\ {\partial D\over \partial b} &=\sum_i 2(y_i-a-bx_i)\cdot(-x_i)=0 \end{align} \]

用\(n\bar x=\sum_i x_i, n\bar y=\sum_i y_i\)代入,可得:

\[\begin{cases} a=\bar y-b\bar x \\ b={\sum_i x_iy_i-\bar y\sum_i x_i\over \sum_i x_i^2-\bar x\sum_i x_i}={\sum_i (x_i-\bar x)(y_i-\bar y)\over \sum_i (x_i-\bar x)^2} \end{cases} \]

公式计算参考

多元线性回归-Normal Equation

Normal Equation是一种基础的最小二乘方法，是以矩阵计算的方式求解多元线性回归。相比于梯度下降法是一种更加数学化的方法，可以用来求权重向量θ。

记：

\[\begin{align} \underset{1\le i\le m}{x^{(i)}} = \left[ \begin{matrix} x_0^{(i)}\\ x_1^{(i)}\\ x_2^{(i)}\\ ...\\ x_n^{(i)} \end{matrix} \right] \in R^{(n+1)};\quad y=\left[ \begin{matrix} y^{(1)}\\ y^{(2)}\\ y^{(3)}\\ ...\\ y^{(m)} \end{matrix} \right]; \quad \theta= \begin{bmatrix} \theta_0\\ \theta_1\\ \theta_2\\ ...\\ \theta_n \end{bmatrix} \\ X=\left[ \begin{matrix} (x^{(1)})^T\\ (x^{(2)})^T\\ (x^{(3)})^T\\ ...\\ (x^{(m)})^T \end{matrix} \right]= \left[ \begin{matrix} x_0^{(1)}&x_1^{(1)}&x_2^{(1)}&...&x_n^{(1)}\\ x_0^{(2)}&x_1^{(2)}&x_2^{(2)}&...&x_n^{(2)}\\ x_0^{(3)}&x_1^{(3)}&x_2^{(3)}&...&x_n^{(3)}\\ ...&...&...&...&...\\ x_0^{(m)}&x_1^{(m)}&x_2^{(m)}&...&x_n^{(m)}\\ \end{matrix} \right] \end{align} \]

用矩阵相乘的形式表示为\(h_{\theta}(x)=\theta^TX=X\cdot\theta\),均方误差损失函数\(J(\theta_0,\theta_1,...,\theta_n)=\frac{1}{2m} \sum_{i=1}^m (h_\theta(x^{(i)})-y^{(i)})^2\)可以表示为

\[\begin{align}J(\theta)&=((X\theta)^T-y^T)(X\theta-y) \\ &=(X\theta)^TX\theta-(X\theta)^Ty-y^T(X\theta)+y^Ty \\ &=\theta^TX^TX\theta-2(X\theta)^Ty+y^Ty \end{align} \]

接着对θ求导：

\[\begin{align} \frac{\partial J}{\partial \theta}=2X^TX\theta-2X^{T}y=0 \\ X^TX\theta=X^{T}y \\ \theta=(X^TX)^{-1}X^Ty \end{align} \]

最终得到的计算表达式即为Normal equation：\(\theta=(X^TX)^{-1}X^Ty\).

此时表达式中存在一个问题，如果矩阵\((X^TX)\) 不可逆怎么办呢？原因有二：

1.存在冗余特征。
比如有两个特征

特征 \(x_1\) \(x_2\)
意义 面积(\(feet^2\)) 面积(\(m^2\))

因为 1 m=3.28 feet，所以 \(x_1=(3.28)^2\cdot x_2\)。我们知道线性代数中出现这种线性相关的情况，其行列式值为0，所以不可逆，我们只需确保不会出现冗余特征即可。(根据相关性去除冗余特征)

2.特征数量 n 过多，而训练样本数 m 过少。 解决方法为删除一部分特征，或者增加样本数量。(特征选择)

另外的解决办法是加入扰动项防止\(X^TX\)不可逆，参考 再论最小二乘。

增加λ扰动：\(\theta = (X^{T}X + \lambda I)^{-1}X^{T}y\)同时能够起到防止过拟合的作用. 这个正则项与添加L2 norm的线性回归求导结果相同:
\(J'(\theta) = [\frac{1}{2}(h_{\theta}(x^{(i)})-y^{(i)})^{2}]' + \lambda\sum_{j=1}^{n}\theta_{j}^{2}= \frac{1}{2}(X\theta-y)^{T}(X\theta-y) + \lambda \theta^{T}\theta\)

然后对这个新的目标函数求驻点，得到：\(\theta = (X^{T}X + \lambda I)^{-1}X^{T}y\)

加权最小二乘（WLS）

注：这里说的加权最小二乘与局部加权线性回归 (LWR) 不同的是，每个变量\(x_i\)具有固定的权重\(w_i\).

引入对角矩阵W：

\[W=\left[\begin{array}{cccc} w_{11} & \ldots & 0 & 0 \\ 0 & w_{22} & \ldots & 0 \\ \ldots & \ldots & \ldots & \ldots \\ 0 & 0 & \ldots & w_{m m} \end{array}\right] \]

损失函数变为：\(f(\mathbf{\theta})=\frac{1}{2}\left \| W(X\mathbf{\theta}-y) \right \|^{2}\)

容易看出只需在Normal Equation求解公式中将\(X\)替换为\(WX\)、\(y\) 替换为\(Wy\)，即可得到加权最小二乘Normal equation：

\[\theta=(X^TW^TWX)^{-1}X^TW^TWy=(X^TW^2X)^{-1}X^TW^2y \]

迭代再加权最小二乘(IRLS)

Iteratively reweighted least squares - Wikipedia

迭代再加权最小二乘（IRLS）或称迭代变权最小平方差（IWLS）算法用于解决广义线性回归的最优化问题，这个最优化问题的目标函数(p-norm 形式)如下所示：

\[\arg \min_{\beta} \sum_{i=1}^{n}|y_{i} - f_{i}(\beta)|^{p} \]

这个目标函数可以通过迭代的方法求解。在每次迭代中，解决一个带权最小二乘问题。在 t +1 步采用Normal Equation求解加权线性最小二乘法 问题。形式如下：

\[\beta ^{t+1} = \arg \min_{\beta} \sum_{i=1}^{n} w_{i}(\beta^{(t)}))|y_{i} - f_{i}(\beta)|^{2} = (X^{T}W^{(t)}X)^{-1}X^{T}W^{(t)}y \]

在这个公式中，\(W^{(t)}\)是对角矩阵（参考加权最小二乘Normal equation中的 \(W^2\)），它的所有元素都初始化为1。每次迭代中，通过下面的公式更新：

\[W_{i}^{(t)} = |y_{i} - X_{i}\beta^{(t)}|^{p-2} \]

原理也比较简单, 权重是p-2次项, 作用在线性回归2次项之上, 得到的便是p次项。每次迭代根据前一次迭代得到的β值计算当前轮的W，继而计算当前轮的β值。

IRLS for LR

逻辑回归(logistic regression)不像线性回归那样有闭解（closed-form solution）, 因为sigmoid函数带来了非线性。

LR也可用IRLS来求解，参考 《PRML》Logistic回归的IRLS求解

IRLS 优点:

IRLS 常用于解决介于L1-norm 与 L2-norm 之间的问题，即范数p取值 \(1\le p \le 2\). IRLS 是一种计算近似 L1 问题的简单方法, L1-norm solution 通常相比于 L2, 对高振幅的噪声更不敏感。

参考：

线性回归-梯度下降

随机初始化θ；
沿着负梯度方向迭代更新θ，为了使得 J(θ) 更小:
\(\theta_{j} = \theta_{j} - \alpha \frac{\partial J(\theta)}{\partial \theta}\)
梯度方向为:

\[\begin{align} \frac{\partial J(\theta)}{\partial\theta_{j}} &= \frac{\partial}{\partial\theta_{j}}\frac{1}{2}(h_{\theta}(x) - y)^{2} \\ &= (h_{\theta}(x) - y)\frac{\partial}{\partial\theta_{j}}(h_{\theta}(x) - y) \\ &= (h_{\theta}(x) - y)x_{j} \end{align} \]

梯度下降和Normal equation的区别

假设我们有 m 组样本，每个样本有 n 个特征

对比 Gradient descent Normal equation
超参数 需要选取合适的学习率α 不需要学习率
迭代计算 需要很多次的迭代计算 不需要迭代计算
计算量 在n很大时，能工作得很好 需要计算\((X^TX)^{−1}\)，其时间复杂度为\(O(n^3)\);
当n非常大时，速度非常慢

我们如何来选择使用哪种方法么呢？这个就要跟着感觉走了。因为我们无法精确的得到选择Gradient descent还是Normal equation的阈值。一般如果 n 不超过1000的话，选用Normal equation还是挺好的，如果要是 n 超过10000的话，那么矩阵 X 的规模将是10000×10000，这时候Gradient descent将是一个更好的选择。

牛顿法求最大似然估计

解最大似然估计的方法可以是求导迭代的方法,但是收敛速度慢.而使用牛顿下降法可以加快收敛速度.

当要求解\(f(\theta)=0\)时,如果f可导,那么可以通过迭代公式:

\[\theta:=\theta-\frac{f(\theta)}{f^{'}(\theta)} \]

来求解最小值.(\(\theta_{k+1}\)为点(\(\theta_k,f(\theta_k)\))处的切线与横轴的交点)

当应用于求解最大似然估计的最大值时,变成求解 \(l^{'}(\theta)=0\) 的问题,迭代公式变为:

\[\theta:=\theta-\frac{l^{'}(\theta)}{l^{''}(\theta)} \]

当\(\theta\)是向量时,牛顿法可以使用下面的式子表示:

\[\theta:=\theta-H^{-1}\nabla_\theta l(\theta) \]

其中\(H_{ij}=\frac{\partial^2l(\theta)}{\partial\theta_i \partial\theta_j}\)是\(n\times n\)的Hessian矩阵.

牛顿法收敛速度虽然很快,但求 Hessian 矩阵的逆的时候比较耗费时间。
当初始点 \(X_0\) 靠近极小值 X 时,牛顿法的收敛速度是最快的。但是当 \(X_0\) 远离极小值时,
牛顿法可能不收敛,甚至连下降都保证不了。原因是迭代点 \(X_{k+1}\) 不一定是目标函数 f 在牛顿方向上的极小点。

正则化

对误差项(损失函数)加惩罚项是防止过拟合的一种手段.
加L1正则的线性回归有lasso回归,趋向于迫使系数变为0,这使得它更适合于学习稀疏模型.
加L2正则的回归有岭回归(ridge regression),迫使系数缩小.

系数非负的线性回归

非负最小二乘回归, 其实是带不等式条件的二次规划（Quadratic Programming），可以用二次规划来求解：

\(\arg\min_w|Xw−y|^2\ \ s.t.\ w>0\)，

二次规划参考：

在sklearn Lasso回归 sklearn.linear_model.Lasso 里可以设置 positive=True 来使得系数都为正。如果不需要正则化将 alpha 设成0或者很小的数。

from sklearn.linear_model import Lasso
model = Lasso(alpha=0.00001, positive=True)
model.fit(X,y)

参考：sofasofa

posted @ 2017-10-12 11:12  康行天下  阅读(1088)  评论(0)  编辑  收藏  举报