线性模型、最优化方法（二）

一、线性回归

一般的，线性回归模型表示为

\[h_{\theta}(x)=\theta_0+\theta_1x_1+...+\theta_nx_n=\sum_{i=0}^{n}\theta_ix_i=\theta^Tx \]

上式中令\(x_0=1\),这样\(x\)实际上是\(n+1\)维，\(x=[1,x_1,x_2,…,x_n]^T\).当我们知道训练数据集后怎样得到参数\(\theta\)的值呢？一个比较合理的想法是尽可能的让\(h_{\theta}(x)\)接近真实值\(y\),下面我们来定义一个函数来度量\(h_{\theta}\)与\(y\)的距离：

\[J(\theta)=\frac{1}{2}\sum_{i=1}^m(h_{\theta}(x^{(i)})-y^{(i)})^2 \]

其中\(m\)表示训练集中样本个数，\(x^{(i)}\)和\(y^{(i)}\)分别表示第\(i\)个样例的特征值和真实结果。我们希望找到一个\(\theta\)使的\(J(\theta)\)最小。下面分别采用梯度下降算法和最小二乘法来求解\(\theta\)。

二、梯度下降算法

我们希望能找到\(\theta\)使得\(J(\theta)\)达到最小，于是我们可以使一个搜素算法，初始化\(\theta\)为一个任意值，在迭代过程中不断的更新\(\theta\)使得\(J(\theta)\)更小，直到收敛为止（可以认为\(J(\theta)\)不再改变）。下面我们来考虑下梯度下降算法，它先给\(\theta\)一个初始化值，然后通过下面公式来更新\(\theta\)

\[\theta_j:=\theta_j-\alpha\frac{\partial}{\partial \theta_{j}}J(\theta)\quad 其中 \quad j = 0,1,...,n \]

上式中\(\alpha\)称为学习速率。为了理解上面的算法，我们需要推导一下偏微分部分，推导过程中我们只考虑一个样例(x,y),然后只需要把所有的加起来就可以了

\[\frac{\partial}{\partial \theta_{j}}J(\theta)=\frac{1}{2}\frac{\partial}{\partial \theta_{j}}(h_{\theta}(x)-y)^2 \\ =(h_{\theta}(x)-y)\frac{\partial}{\partial \theta_{j}}(\sum_{i=0}^n\theta_i x_i-y) \\ =(h_{\theta}(x)-y)x_j \quad \quad \quad \quad \quad \quad \quad \]

对于单独的第\(i\)个样本，更新规则为

\[\theta_j:=\theta_j-\alpha(h_{\theta}(x^{(i)})-y^{(i)})x_{j}^{(i)}\quad 其中 \quad j = 0,1,...,n \]

当有\(m\)个样本时，更新规则为

\[\theta_j:=\theta_j-\alpha \sum_{i=1}^m(h_{\theta}(x^{(i)})-y^{(i)})x_j^{(i)}\quad 其中 \quad j = 0,1,...,n \]

上式每次迭代将所有样本考虑进来更新\(\theta\),这种方法叫做批量梯度下降。梯度下降容易陷入一个局部的最优解，事实上\(J(\theta)\)是一个凸函数，没有所谓的局部最优解，它只有一个全局最优解。

​ 批量梯度下降一般工作的很好，考虑下面的算法：

该算法每一个样例都更新了\(\theta\)一次，这个算法被称作随机梯度下降（或增量梯度下降）SGD，当\(m\)很大的时候，批量梯度下会很耗时，通常，随机梯度下降比批量梯度下降要快很多，然而它可能永远达不到最小值，并在最小值附近震荡，而实际上这种近似最小值已经够了。在样本量比较大的情况下，随机梯度下降比梯度下降更实用。

三、最小二乘法

梯度下降算法给出了一种最小化\(J(\theta)\)的方法。这里将介绍第二种方法：最小二乘法。我们先求\(J(\theta)\)的偏导数，然后令偏导数为0，用等式解出\(\theta\)，此时的\(\theta\)使得\(J(\theta)\)最小（\(J(\theta)\)是下凸函数）。

​ 先介绍下矩阵迹相关知识，对于一个\(n \times n\)的方阵\(A\),定义它的迹为

\[tr(A)=\sum_{i=1}^nA_{ii} \]

即矩阵对角线的和。实数的迹是它本身，有关迹的一些性质

\[tr(A)=tr(A^T) \quad \quad \quad\quad\quad\quad\quad\quad\\ tr(A+B)=tr(A)+tr(B)\quad \quad \quad \\ tr(kA)=ktr(A)\quad \quad \quad\quad\quad\quad\quad\;\;\\ tr(AB)=tr(BA)\quad \quad \quad\quad\quad\quad\quad\;\; \\ tr(ABC)=tr(CAB)=tr(BCA) \quad\\ tr(ABCD)=tr(BCDA)=tr(CDBA)=tr(DABC) \quad\\ \]

不加证明的给出有关矩阵微分结论

\[\nabla_A tr(AB)=B^T \quad \quad \quad \quad \quad \quad \quad \quad \\ \nabla_{A^T}f(A)=(\nabla_Af(A))^T \quad \quad \quad \quad \quad\\ \nabla_Atr(ABA^TC)=CAB+C^TAB^T \\ \nabla_A |A|=|A|(A^{-1})^T \quad \quad \quad \quad\quad \quad\;\; \]

其中\(|A|\)表示\(A\)的行列式，第四个式子中\(A\)可逆。定义\(X\)为\(m \times n\)（实际上是\(m \times (n+1)\)）的矩阵，用它表示训练的数据特征集

\[X=[x^{(1)},x^{(2)},...,x^{(m)}] \]

定义真实值向量\(y=[y^{(1)},y^{(2)},…,y^{(m)}]^T\),于是

\[X^T\theta - y =[h_{\theta}( x^{(1)})-y^{(1)},h_{\theta}( x^{(2)})-y^{(2)},...,h_{\theta}( x^{(m)})-y^{(m)}]^T \]

我们知道对于一个\(n\)维向量\(z\)有\(z^Tz=\sum_{i=1}^{n}z_{i}^2\)，于是

\[J(\theta)=\frac{1}{2}\sum_{i=1}^m(h_{\theta}(x^{(i)})-y^{(i)})^2=\frac{1}{2}(X^T\theta-y)^T( X^T\theta-y) \]

下面来对\(J(\theta)\)求导

\[\nabla_{\theta} J(\theta)=\frac{1}{2}\nabla_{\theta}(X^T\theta-y)^T( X^T\theta-y) \quad\quad\quad\quad\quad\quad\quad\quad\\ =\frac{1}{2}\nabla_{\theta}(\theta^TXX^T\theta-2\theta^T Xy+y^Ty) \quad\quad \;\;\\ =\frac{1}{2}\nabla_{\theta}(tr\;\theta^TXX^T\theta-2\;tr\;\theta^T Xy+y^Ty) \\ =\frac{1}{2}\nabla_{\theta}tr\;\theta\theta^TXX^T-(\nabla_{\theta^T} tr\;\theta^T Xy)^T \quad \\ =XX^T\theta-Xy \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\;\; \]

令\(\nabla_{\theta} J(\theta)=0\)得\(XX^T\theta=Xy\),求解

\[\theta = (XX^T)^{-1}Xy \]

这样我们直接求出来了\(\theta\)的解析解，这就是最小二乘法。

四、有关概率的解释

当我们求解回归问题时，为什么选择线性回归？为什么代价函数\(J\)是合理的呢？这里我们将给出一个概率性的假设，然后很自然的推导出最小二乘法。

我们结果的目标值和预测值满足以下等式

\[y^{(i)}=\theta^Tx^{(i)}+\epsilon^{(i)} \]

其中\(\epsilon^{(i)}\)表示预测误差,我们再继续假设\(\epsilon^{(i)}\)满足独立同分布，它们服从均值为0，方差为\(\delta\)的高斯分布\(\epsilon ^{(i)}\sim N(0,\delta^{2})\),概率密度函数即为

\[p(\epsilon ^{(i)})=\frac{1}{\sqrt{2 \pi}\delta}exp(-\frac{(\epsilon^{(i)})^2}{2\delta}) \]

于是有\(y^{(i)}\sim N(\theta^Tx^{(i)},\delta^2)\)分布，则

\[p(y^{(i)}|x^{(i)} ;\theta)=\frac{1}{\sqrt{2 \pi}\delta}exp(-\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\delta}) \]

需要注意\(p(y^{(i)}|x^{(i)} ;\theta)\)中德\(\theta\)是一个确定的值，并不是随机变量，这个和\(p(y^{(i)}|x^{(i)} ,\theta)\)是不同的，此处的\(\theta\)是一个随机变量。我们希望观测概率最大，也就是关于观测值的似然函数

\[L(\theta)=L(\theta;x,y)=p(y|x;\theta) \]

注意前面我们假设了随机变量\(\epsilon^{(i)}\) 是相互独立的，所以

\[L(\theta)=\prod_{i=1}^mp(y^{(i)}|x^{(i)})=\prod_{i=1}^m\frac{1}{\sqrt{2 \pi}\delta}exp(-\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\delta}) \]

现在给出\(y^{(i)},x^{(i)}\),如何选择参数\(\theta\)使得\(L(\theta)\)最大呢？这中方法被称极大似然估计。我们取对数似然函数

\[l(\theta)=ln\;L(\theta)=\frac{m}{ln\;(\sqrt{2\pi}\delta)}-\sum_{i=1}^m\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\delta} \]

我们需要最大化\(l(\theta)\)也就是得极小化

\[\sum_{i=1}^m\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2}=J(\theta) \]

我们总结一下，在概率假设前提下，最小二乘回归就是找到\(\theta\)的最大似然估计。

五、局部加权线性回归

在原来的线性回归中，为了预测在点\(x\)处的值，我们会

• 拟合\(\theta\)来最小化\(\sum_i(y^{(i)}-\theta^Tx^{(i)})^2\)
• 输出\(\theta^Tx\)

而在局部线性回归中，我们会

• 拟合\(\theta\)来最小化\(\sum_iw^{(i)}(y^{(i)}-\theta^Tx^{(i)})^2\)
• 输出\(\theta^Tx\)

这里的\(w^{(i)}\)是非负数的权重，直观上看，如果\(w^{(i)}\)很小，将会被忽略\((y^{(i)}-\theta^Tx^{(i)})^2\),下面是\(w^{(i)}\)一个常用的权重

\[w^{(i)}=exp(-\frac{(x^{(i)}-x)^2}{2\tau^2}) \]

这个可以内比高斯分布的图形，但事实上它们没半毛钱关系。距离\(x\)越远的点权重越小，越近的点权重越大，当\(|x^{(i)}-x|\)越小，\(w^{(i)} \to 1\),当\(|x^{(i)}-x|\)很大时，\(w^{(i)} \to 0\).参数\(\tau\)通常称为波长参数，控制权重函数图像的形状。

​ 局部加权线性回归是我们看到的第一个非参数学习模型，我们之前看到的是未加权的参数学习算法，因为它具有固定个数的参数，这些参数只需要学习一遍。相反的，使用局部加权线性回归，我们需要保持整个训练集，每次预测都要用不同的参数。

六、逻辑回归模型

现在我们来考虑下分类问题。这和回归问题一样，唯一的区别是我们要预测的值是离散的。现在我们来讨论二元分类问题，其中\(y\)能只能取\(0\)或者\(1\).例如现实生活中，邮件是否是垃圾邮件，只有是和否两种答案。通常我们也可以称这样的样本为正或负样本。

​ 上述问题，我们忽视\(y\)是离散值的事实，也可以采用线性回归来分类，但是\(y \in \{0,1\}\),线性回归得到的大于1或者小于0的值都是没有意义的。为了解决这个问题，我们修改假设函数\(h_{\theta}(x)\)的形式

\[h_{\theta}(x)=g(\theta^Tx)=\frac{1}{1+exp(-\theta^Tx)} \]

其中

\[g(z)=\frac{1}{1+e^{-z}} \]

称为\(sigmiod\)函数着或logistic函数。sigmoid函数将\(\theta^Tx\)映射为\(0\to 1\)的概率，sigmoid函数还有一个很好的性质

\[g(z)'=g(z)(1-g(z)) \]

那么给出了逻辑回归模型，我们如何得到参数\(\theta\),和前面最小二乘法概率解释做法一样，可以采用极大似然估计方法来求解参数\(\theta\).我们假设

\[P(y=1|x;\theta)=h_{\theta}(x) \quad\quad \\ P(y=0|x;\theta)=1-h_{\theta}(x) \]

注意到上面两个式子可以写成

\[P(y|x;\theta)=h_{\theta}(x)^y(1-h_{\theta}(x))^{1-y} \]

我们假设有\(m\)个训练样本相互独立，似然函数为

\[L(\theta)=\prod_{i=1}^mp(y|x;\theta)\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad \\ =\prod^{m}_{i=1}h_{\theta}(x^{(i)})^{y^{(i)}}(1-h_{\theta}(x))^{1-y^{(i)}} \]

然后再化简为对数似然函数

\[l(\theta)=\sum_{i=1}^m [y^{(i)}ln\;h_\theta(x)+(1-y^{(i)})ln\;(1-h_\theta(x))] \]

要最大化\(l(\theta)\)，我们可以向前面那样，使用梯度上升算法根据下面迭代公式

\[\theta := \theta +\alpha \nabla_\theta l(\theta) \]

更新\(\theta\).我们用一个样例来求导，然后用随机梯度上升算法求解。

\[\frac{\partial}{\partial \theta_j}l(\theta)=[\frac{y}{g(\theta^Tx)}-\frac{1-y}{1-g(\theta^Tx)}]\frac{\partial}{\theta_j}g(\theta^Tx) \\ =\frac{y-g(\theta^Tx)}{g'(\theta^Tx)}g'(\theta^Tx)x_j=(y-h_\theta(x))x_j \]

这个形式和最小二乘法得到的一样，区别是模型\(h_\theta(x)\)不一样。我们得到的随机梯度上升迭代公式为

\[\theta_j := \theta_j +\alpha (y^{(i)}-h_\theta(x^{(i)}))x_j \quad foreach \;j \;from \; 1 \to n \]

七、感知机模型

我们现在来讨论一个具有历史意义的算法，考虑修改逻辑回归算法让结果强制映射为0和1，修改\(sigmoid\)函数为下面形式

\[g(z)=\left \{ \begin{matrix} 0 & z<0\\ 1 & z \ge 0 \end{matrix} \right . \]

于是

\[h_\theta(x)=g(\theta^Tx)=\left \{ \begin{matrix} 0 & \theta^Tx<0\\ 1 & \theta^Tx \ge 0 \end{matrix} \right . \]

我们得到的\(\theta\)学习规则与前面线性模型一样

\[\theta_j=\theta_j+\alpha(y^{(i)}-h_\theta(x^{(i)}))x_j^{(i)} \]

这样，我们就得到了感知机学习算法。在过去，感知机被粗略的认为是大脑中的一个简单的神经元，需要注意尽管感知机算法在形式上和前面讨论的算法相似，但它与逻辑回归和最小二乘回归实际上是不同的模型。特别是将感知机赋予概率的意义，或者是利用极大似然估计来推导它。

八 牛顿法

以逻辑回归为例，来讨论另一个最优化算法：牛顿法。下面使用牛顿法来最大化\(l(\theta)\).牛顿法的思想就是通过迭代来找到函数的零点。特别的，假设我们有一个函数\(f:\textbf{R} \mapsto \textbf{R}\),我们想找到一个\(\theta\)满足\(f(\theta)=0\),这里\(\theta \in R\),牛顿法的迭代公式为

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

这种方法通常有一个直观的解释，我们可以通过一个线性函数来近似\(f\),该线性函数在当前猜测的\(\theta\)点处与\(f\)相切，接下来求解线性函数零点，并让下一个猜测点在线性函数的零点位置，如下图

牛顿迭代法给出了一个求解零点的方法，但是我们怎么来通过牛顿迭代法来优化\(l(\theta)\)呢？求解\(l(\theta)\)的最大值先要求解\(l'(\theta)\)的零点，所以上式中我们令\(f(\theta)=l'(\theta)\),于是

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

在逻辑回归模型中，\(\theta\)是一个向量而不是一个实数，所以我们需要将牛顿法推广到多维的情况，推广公式为

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

其中，\(\nabla_\theta l(\theta)\)是\(l(\theta)\)的梯度，\(H^{-1}\)是\(n\)阶方阵（实际上是\(n+1\)阶），他被称为\(Hessian\)矩阵，它的元素其实就是\(l(\theta)\)二阶偏导数组成的

\[H^{-1}_{i,j}=\frac{\partial^2 \theta}{\partial \theta_i \partial \theta_j} \]

牛顿法是二阶收敛，梯度下降是一阶收敛，牛顿法通常比梯度下降算法收敛速度要快，只需要更少的迭代次数就能获得最小值。然而一次牛顿迭代要比一次梯度下降更昂贵，因为它需要计算\(Hessan\)矩阵并且求它的逆，这将花费不少时间。但是当参数个数\(n\)不是太大时，总体来说速度还是要快很多。