第四课 牛顿方法、指数族和广义线性模型

1. 摘要

本课有三块内容,首先补充讲述了上次没有完成的牛顿方法,再接着介绍了指数族的概念,然后推出了广义线性模型。

2. 牛顿方法(Newton's method)

2.1 简单回顾

上次课的最后我们讲到了一种分类算法——Logistic回归,并且给出了它的假设函数和似然函数:

h_{\theta}(x) = P(y=1|x;\theta) = \frac{1}{1 + e^{- \theta x}} \\\\ l(\theta) = \sum_{i=1}^{m}(y^{(i)}ln(h_{\theta}(x^{(i)})) + (1-y^{(i)})ln(1-h_{\theta}(x^{(i)})))

对上述似然函数求导,并且应用梯度上升算法,我们可以推出Logistic模型的参数\theta的递推公式:

\theta_{j} := \theta_{j} + \alpha \sum_{i=1}^{m}(y^{(i)} - h_{\theta}(x^{(i)}))x_{j}^{(i)}

上次课我们使用的是批梯度上升,其算法执行效率可能不高,我们可以采取随机梯度下降来提高效率:

\theta_{j} := \theta_{j} + \alpha (y^{(i)} - h_{\theta}(x^{(i)}))x_{j}^{(i)}

2.2 牛顿方法的思想

实际上,牛顿方法的作用和梯度上升是一样的,但是它的运行速度通常比梯度上升快很多,这一点我们后续还会说到。下面先来了解一下牛顿方法的思路。
首先我们考虑这样一个问题,给定一个函数f(x),假设其图像如下,我们如何求它的零点,即求得使f(x) = 0x的值。

我们随意选择函数上一点x_0,做这一点的切线,该切线通常会和x轴有一个交点(这里暂不考虑那些特殊情况),记为x_1

然后对于x_1,我们重复上述操作,可以得到x_2,以此类推,我们就可以求得函数的零点。
那么结合之前内容,我们为了使得似然函数达到最大值,首先想到的方法就是求其导函数,然后找到导函数值为0的那一点对应的\theta值,就是我们要的结果。
所以我们可以用上面的算法,求得导函数值的零点,也同样可以得到我们要的结果。

2.3 牛顿方法具体使用

了解过基本思路,下面我们具体来说说如何使用牛顿方法。
我们还是先选取x_0,从而我们可以求得这一点的函数值f(x_0)和导数值f'(x_0),那我们如何找到x_1呢?
其实x_0x_1的关系从下图就可以很直观地看出来。

x_0导数值就是这一点切线的正切值,所以,x_0x_1的变化量\Delta就可以这样算出:

\Delta = \frac{f(x_0)}{f'(x_0)}

那么x_1就等于:

x_1 = x_0 - \Delta = x_0 - \frac{f(x_0)}{f'(x_0)}

接着就可以写出递推公式:

x_t = x_{t-1} - \Delta = x_{t-1} - \frac{f(x_{t-1})}{f'(x_{t-1})}

按照上面提到的思路,我们要求的是似然函数倒数的零点,即l'(\theta) = 0\theta值,因此上述递推公式就变为:

\theta_t = \theta_{t-1} - \frac{l'(\theta_{t-1})}{l''(\theta_{t-1})}

这就是牛顿方法应用的过程。

2.4 牛顿方法推广

刚才我们提到牛顿方法的执行效率是很高的,事实上,牛顿方法的收敛速度被称为二次收敛,而梯度下降算法的收敛速度为一次收敛,我们可以想象一下,二次函数和一次函数发散到正无穷的速度,就可以对这两个收敛效率有一个直观的感受了。
我们刚才提到的\theta只是一个实数,而我们之前遇到的问题,\theta往往是一个向量。因此,我们需要把上述递推公式修改为如下形式:

\theta_{t} = \theta_{t-1} - H^{-1}\nabla_{\theta} l

其中H被成为Hessian矩阵,可以理解为对函数求二阶导,具体定义如下:

H_{ij} = \frac{\partial^2l}{\partial \theta_{i} \partial \theta_{j}}

2.5 牛顿方法的原理

这一部分的内容在视频中并没有提到,但我相信每一位愿意追根究底的同学都会问这样一个问题:为什么牛顿法可以找到极值点?它对函数又有什么要求?
下面的论述并不严格,主要考虑现实情况,如有硬伤,欢迎各位读者进行指正。
首先,我们假设目标函数f(x)在其极值点x^*的邻域内二阶可导。接着,\forall x_0 \in N(x^*, \delta),我们在x_0处对f(x)进行二阶泰勒展开:

f(x) = f(x_0) + \frac{f'(x_0)}{1!}(x - x_0) + \frac{f''(x_0)}{2!}(x - x_0)^2 + o(x^3)

对上式两边关于x求导:

f'(x) = f'(x_0) + f''(x_0)(x - x_0)

由于x^*是极值点,从而该点导数值为0,那么:

\begin{aligned}f'(x^*) & = f'(x_0) + f''(x_0)(x^* - x_0) \\\\ 0 & = f'(x_0) + f''(x_0)(x^* - x_0) \\\\ x^* & = x_0 - \frac{f'(x_0)}{f''(x_0)}\end{aligned}

但是,上式是由泰勒展开得到的,只能说上式中的x^*x_0更接近极值点,于是我们可以记上式中的x^*x_1,继续用x_1进行泰勒展开,最终得到迭代公式:

x_n = x_{n-1} - \frac{f'(x_n)}{f''(x_n)}

到这里,我们开始提出的两个问题也就算解答了,对于一个二阶可导的函数,我们不断用泰勒展开将非线性的函数用二次多项式进行逼近,从而得到我们要的结果。

3. 指数族

3.1 指数族分布的概念

到现在为止,我们已经见过了一些概率分布,比如正态分布和伯努利分布。先以伯努利分布为例,我们只要定义P(y = 1; \Phi) = \Phi就可以得到一个以\Phi为参数的伯努利分布,所以可以很自然地想到,我们不断改变\Phi的值,就可以得到一组伯努利分布。同理,如果我们改变正态分布中\sigma\mu的值,同样可以得到一组不同的正态分布。然而事实上,我们可以进一步推广这个过程,定义一种更加广泛的分布形式,只要我们改变其中的参数,我们可以得到不同种类的分布,这就是指数族分布的思想。其具体定义如下:

\begin{aligned} \mathbb{如果一种分布可以写为如下形式:} & \\ & P(y; \eta) = b(y)exp(\eta^Ty - a(\eta)) \\ \mathbb{我们就称该种分布是指数族分布} & \end{aligned}

接着解释一下其中两个符号:

\begin{aligned}\eta & : \text{自然参数} \\ T & : \text{充分统计量。在很多时候,T(y) = y} \end{aligned}

从而,当我们给定一组a\text{、}b和T的时候,我们就可以得到一组以\eta为参数的概率分布。

3.2 伯努利分布写成指数族分布形式

回忆我们之前讲过的,伯努利分布可以定义成如下形式:

P(y; \Phi) = \Phi^y(1 - \Phi)^{1 - y}

还是使用先取对数再取指数的技巧:

\begin{aligned}P(y; \Phi) & = exp(ln(\Phi^y(1 - \Phi)^{1 - y})) \\ & = exp(yln(\Phi) + (1 - y)log(1 - \Phi)) \\ & = exp(yln(\Phi) + log(1 - \Phi) - ylog(1 - \Phi)) \\ & = exp(yln(\frac{\Phi}{1 - \Phi}) + log(1 - \Phi)) \\ \end{aligned}

到此,我们就可以定义:

\begin{aligned}\eta & = ln(\frac{\Phi}{1 - \Phi}) \\ T(y) & = y \\ a(\eta) & = log(1 - \Phi) \\ b(y) & = 1 \\ \end{aligned}

于是伯努利分布化为了指数族分布的形式。

3.3 正态分布写成指数族分布形式

正态分布有两个参数,\mu\sigma,其参数本应该是一个二元组或向量,但在我们之前使用到的模型中,只考虑了期望\mu而并没有用到\sigma,为了简化推导,我们直接指定\sigma = 1,因此,正态分布的密度函数变成如下形式:

\begin{aligned} P(y; \mu) & = \frac{1}{\sqrt{2\pi}}exp(- \frac{(y - \mu)^2}{2}) \\ & = \frac{1}{\sqrt{2\pi}}exp(- \frac{y^2 + \mu^2 - 2y\mu}{2}) \\ & = \frac{1}{\sqrt{2\pi}}exp(- \frac{1}{2}y^2)exp(\mu y - \frac{1}{2}\mu^2) \end{aligned}

其实写到这里,童鞋们就已经能看出来了,在上面的式子中,我们定义:

\begin{aligned}\eta & = \mu \\ T(y) & = y \\ a(\eta) & = \frac{1}{2}\mu^2 \\ b(y) & = \frac{1}{\sqrt{2\pi}}exp(- \frac{1}{2}y^2) \end{aligned}

至此,正态分布也化为了指数族分布形式。

4. 广义线性模型(GLM)

4.1 广义线性模型的基本概念

首先,广义线性模型有以下三个假设:

  1. 模型分布是以\eta为参数指数族分布,
    即:y|x;\theta \sim ExponentialFamily(\eta)
  2. 对于给定的x,我们的目标是输出E[T(y)|x],或称:h(x) = E[T(y)|x]
  3. 参数\etax是线性关系,\eta = \theta^Tx。通常来说,\eta是一个实数,当然对于\eta是向量这种更一般的情况,我们可以定义\eta_{i} = \theta_{i}x

4.2 用伯努利分布的广义线性模型推导Logistic模型

结合以上假设,我们回忆一下伯努利分布:
首先,伯努利分布是一个指数族分布,而我们的目标是:

h_{\theta}(x) = E[y|x; \theta] = P(y | x; \theta) = \Phi

按照我们之前推导过的\eta\Phi的关系,上面的式子又等于:

h_{\theta}(x) = \Phi = \frac{1}{1 + e^{-\eta}}

再加上第三条假设,继续得到

h_{\theta}(x) = \frac{1}{1 + e^{-\theta^Tx}}

至此,我们就用广义线性模型生成了一个Logistic模型。

此外,视频里还提到了正则响应函数和正则关联函数,但是并没有展开讨论,这里也不多说了。

4.3 多项式分布的广义线性模型

刚才我们推导模型,可以进行一个非黑即白的二值判断,我们可以把它理解为一个分类问题。而更进一步的情况是,一个事物通常是k个类别中的一个,此时,我们就需要多项式分布来对这样的问题进行建模。

4.3.1 基本概念

对于一个k类的分类问题,我们的目标变量y \in {1, \cdots, k},因此我们可以定义\Phi_{1}, \cdots, \Phi_{k},满足:

P(y = i) = \Phi_{i}, i \in {1, \cdots, k}

我们注意到,这k类应该可以构成一个完全事件,因此对于第k种情况,我们可以表示为:

\Phi_{k} = 1 - \sum_{i = 1}^{k - 1}\Phi_{i}

因此,在这个模型中,我们只需要定义k-1个参数,即\Phi_{1}, \cdots, \Phi_{k-1}就可以了。

4.3.2 软软易推导

我们首先将多项式分布写为指数族分布形式,多项式分布的指数族形式与我们上面提到的伯努利分布和正态分布略有区别,其T(y)并不等于y,而是如下定义:

T(1) = \begin{bmatrix} 1 \\ 0 \\ \vdots \\ 0 \\ \end{bmatrix} , T(2) = \begin{bmatrix} 0 \\ 1 \\ \vdots \\ 0 \\ \end{bmatrix} , \cdots , T(k-1) = \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 1 \\ \end{bmatrix} , T(k) = \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ \end{bmatrix} \in \reals^{k-1}

我们可以定义这样一个符号 1\lbrace judgement \rbrace ,当括号中的表达式为真时,其值为1,为假时则值为0。例如

1\lbrace1=3\rbrace=0\, , \quad 1\lbrace\text{PHP是最好的语言}\rbrace=1\, , \quad 1\lbrace\text{凌波丽是女神}\rbrace=1

有了这样一个符号表示,我们可以把最开始的P(y)表示为:

P(y) = \Phi_{1}^{1\lbrace y=1\rbrace}\times\Phi_{2}^{1\lbrace y=2\rbrace}\times\cdots\times\Phi_{k}^{1\lbrace y=k\rbrace}

这里的\Phi_{k}如前所述,是由前k-1个参数表示出来的。同时,按照我们之前定义的T(y),其第i个分量可以表示为T(y)_i = 1\lbrace y=i \rbrace,这样我们继续推导

\begin{aligned} P(y) & = \Phi_{1}^{T(y)_1}\Phi_{2}^{T(y)_2}\cdots\Phi_{k-1}^{T(y)_{k-1}}\Phi_{k}^{1 - \sum_{j=1}^{k-1}T(y)_j} \\ & = exp(T(y)_1ln\Phi_{1} + T(y)_2ln\Phi_2 + \cdots + (1 - \sum_{i=1}^{k-1})ln\Phi_k) \\ & = exp(T(y)_1ln\frac{\Phi_1}{\Phi_k} + T(y)_2ln\frac{\Phi_2}{\Phi_k} + \cdots + T(y)_{k-1}ln\frac{\Phi_{k-1}}{\Phi_k} + ln\Phi_k) \\ & = b(y)exp((\eta^TT(y) - a(\eta))) \end{aligned}

其中:

\begin{aligned} \eta & = \begin{bmatrix} ln(\Phi_1 / \Phi_k) \\ ln(\Phi_2 / \Phi_k) \\ \vdots \\ ln(\Phi_{k-1} / \Phi_k) \end{bmatrix} \\ b(y) & = 1 \\ a(\eta) & = -ln(\Phi_k) \end{aligned}

从而,我们可以从上面的式子中,反解出\Phi_i\eta_i的关系

\begin{aligned} \Phi_i & = \frac{e^{\eta_i}}{1 + \sum_{j = 1}^{k-1}e^{\eta_j}} \\ & = \frac{e^{\theta_i ^Tx}}{1 + \sum_{j = 1}^{k-1}e^{\theta_j ^Tx}} \qquad (use \enspace the \enspace third \enspace assumption, \medspace \eta_i = \theta_i ^Tx) \end{aligned}

再回看我们的拟合函数,可以得到:

\begin{aligned} h_{\theta}(x) & = E[T(y) | x; \theta] \\ & = E \Bigg \lbrack \begin{pmatrix} 1( y = 1 ) \\ 1 ( y = 2 ) \\ \vdots \\ 1( y = k-1 ) \\ \end{pmatrix} \Bigg | x ; \theta \Bigg \rbrack \\ & = \begin{bmatrix} \Phi_1 \\ \Phi_2 \\ \vdots \\ \Phi_{k-1} \\ \end{bmatrix} \\ & = \begin{bmatrix} \frac{e^{\theta_1 ^Tx}}{1 + \sum_{j = 1}^{k-1}e^{\theta_j ^Tx}} \\ \frac{e^{\theta_2 ^Tx}}{1 + \sum_{j = 1}^{k-1}e^{\theta_j ^Tx}} \\ \vdots \\ \frac{e^{\theta_{k-1} ^Tx}}{1 + \sum_{j = 1}^{k-1}e^{\theta_j ^Tx}} \\ \end{bmatrix} \end{aligned}

至此,我们得到了一种判断k类结果的算法,它有一个软软的名字叫Softmax Regression

4.3.3 软软的正确使用姿势

其实Softmax Regression的使用方法跟我们之前用到的回归算法完全一样,基本思路还是通过最大化训练集的似然函数,以此来得到参数\theta
具体过程如下:
假设我们有个分类问题,对于目标变量y,可能取\lbrace 1, 2, 3, \ldots, k\rbrace里面的某个值。
于是,我们收集了训练集\lbrace (x^{(1)}, y^{(1)}), (x^{(2)}, y^{(2)}), \ldots, (x^{(m)}, y^{(m)}) \rbrace,其似然函数表示为:

\begin{aligned} L(\theta) & = \prod_{i = 1}^{m}P(y^{(i)}|x^{(i)}; \theta) \\ & = \prod_{i = 1}^{m}(\Phi_{1}^{1\lbrace y^{(i)} = 1\rbrace}\Phi_{2}^{1\lbrace y^{(i)} = 2\rbrace}\cdots \Phi_{k}^{1\lbrace y^{(i)} = k\rbrace}) \end{aligned}

然后,把\Phi用我们之前计算出来的公式替换,配上前几节我们一直使用的数学技巧,取对数做简化,再求导。

\begin{aligned} l(\theta) & = logL(\theta) \\ & = \sum_{i = 1}^{m}log\prod_{l=1}^{k}\Phi_l^{1\lbrace y^{(i)} = l\rbrace} \\ & = \sum_{i = 1}^{m}log\prod_{l=1}^{k}(\frac{e^{\theta_{l}^Tx^{(i)}}}{\sum_{j=1}^{k}e^{\theta_j^Tx^{(i)}}})^{1\lbrace y^{(i)} = l\rbrace} \end{aligned}