使用梯度下降求解最小二乘

使用梯度下降求解最小二乘

---

简介

在前一篇的总结中，从矩阵的视角回顾了最小二乘，最终还得到了如下求解最小二乘的解的方程：
$x = (W^{T}W)^{-1}W^{T}y \tag{1}$x=(WTW)−1WTy(1)
这个方程本身没有什么问题，只不过也有其局限性，当数据量很大的时候，$W$W也会很大，求解$W^{T}W$WTW本质上是一个不太现实的事情。如果只有1000个样本，那么最终需要求解$1000 \times 1000$1000×1000的矩阵的逆，如果样本有1000个，那么就得求解$10000 \times 10000$10000×10000的矩阵的逆了，对于方阵来说，求解逆矩阵本身就是一个开销比较大的事情，而且，随着矩阵的增大，逆矩阵究竟是否存在都是一个不好说的事情，假如说逆不存在，那么显然这个方法是行不通的了。在前边对于实对称矩阵的总结中说道，实对称矩阵总是可以找到$n$n个标准正交的特征向量，但是也没有说就是一定是满秩的，所以逆矩阵不存在还是很有可能的，而且前边的等式本来就是$W^{T}Wx = W^{T}y \tag{2}$WTWx=WTy(2)而已，所以如果逆矩阵确实不存在，那么也就只能使用等式2来求解了。虽然等式2不是见得有好的实现方法。
在数值计算领域，函数拟合，数据拟合任务，一般都会采用梯度下降算法进行求解，而如果数据量过大，不适合全量梯度下降，也可以使用随机梯度下降算法进行求解。可以说，梯度下降，真是一个万金油方法了。这篇总结文章中，会给出梯度下降求解最小二乘的一个Python实现，同时，给出梯度下降的简易实现和说明。

梯度下降算法步骤说明

在进一步之前，我们先换一种表达方式，前边等式的表达方式换成下边的形式：
$y = Xw \tag{3}$y=Xw(3)这个时候，$X$X矩阵是输入样本矩阵，每一个行向量表示样本，假设每一个样本有$n$n个纬度，而且共有$m$m个样本，那么$X$X就是一个$m \times n$m×n的矩阵了，而需要学习的$w$w权重向量则是一个具有$n$n个代估参数的列向量，可以定义$w_{i}$wi​是$w$w的第$i$i个参数，而$y$y是结果，纬度是$m$m。
梯度下降算法需要定义损失函数，对于最小二乘而言，其实损失函数就是：
$L(w) = \frac{1}{m}\sum_{i=1}^{m}(y'_{i} - y)^{2} \tag{4}$L(w)=m1​i=1∑m​(yi′​−y)2(4)对于梯度下降而言，其实也就是通过计算损失函数对于各个待估参数的偏导，不断迭代更新各个参数，直到损失函数小于某个阈值的过程，对于每一次迭代，其实就是：
$w^{(i+1)} = w^{(i)} - \alpha\nabla_{w}L(w^{(i)}) \tag{5}$w(i+1)=w(i)−α∇w​L(w(i))(5)这里，$\alpha$α是学习率，每次更新权重的比例，而$w^{(i)}$w(i)自然是迭代到第$i$i次，待估权重向量$w$w具体的取值了。每一个损失函数，都是依据不同的权重而不同的，所以每一次更新完权重后，都必须重新计算损失函数的取值。在梯度下降中，会定义一个阈值$tol$tol，当损失少于这个值的时候，就停止了。通常会选择如下的偏导二范数：
$||\nabla_{w}L(w^{(i)})||_{2} &lt; tol \tag{6}$∣∣∇w​L(w(i))∣∣2​<tol(6)当这个条件满足的时候，就可以停止梯度下降了。所以梯度下降的一般过程就是：

• (1) 初始化参数$w$w
• (2) 计算$\nabla_{w}L(w^{(i)})$∇w​L(w(i))
• (3) 判断$||\nabla_{w}L(w^{(i)})||_{2} &lt; tol$∣∣∇w​L(w(i))∣∣2​<tol是否成立，若成立则跳转到步骤(5)，否则执行步骤(4)
• (4) 使用等式$ w^{(i+1)} = w^{(i)} - \alpha\nabla_{w}L(w^{(i)})$更新参数，并跳转到步骤(2)
• (5) 输出待估参数$w$w

可以看到，整个过程并不复杂，只不过，不同的模型，损失函数是不一样的。在任何梯度下降过程中，都是要先定义好模型，模型有各个参数，然后再依据一些原则，比如说最大似然等，得到损失函数，然后在不断迭代知道损失小于某一个阈值为止，可以选用二范数等标准衡量阈值是否满足。这一个过程都是通用的，只不过不同的模型，具体的计算方法都是不一样的。而且学习率的设置也是非常关键的。太大太小都是不行的。对于大样本的情况下，一次取所有数据去计算也是不现实的，所以一般会随机的取部分的样本去计算步骤(2)中的值，这个过程称为随机梯度下降，但是本质上，和梯度下降之间的差别也就是怎么取样本的问题，大体的步骤也就只是这样而已。
关于梯度下降，我想总还是有些局限性的。

代码示例

这篇总结文档中的示例仅仅只是作为演示，去实现最小二乘拟合一些伪造数据，也不具备通用性，只是起到示意作用。
首先，来定义模型，假定我们想要使用下边的方程去拟合数据：
$y = w_{0} + w_{1}x + w_{2}x^{2} \tag{7}$y=w0​+w1​x+w2​x2(7)也可以写作：
$y = Xw \tag{8}$y=Xw(8)那么损失函数就是：
$L(w) = \frac{1}{m}\sum_{i=1}^{m}(w_{0} + w_{1}x + w_{2}x^{2} - y)^{2} \tag{9}$L(w)=m1​i=1∑m​(w0​+w1​x+w2​x2−y)2(9)所以：
$\nabla_{w}L(w) = \begin{bmatrix} \frac{\partial L(w)}{\partial w_{0}} \\ \frac{\partial L(w)}{\partial w_{1}} \\ \frac{\partial L(w)}{\partial w_{2}} \end{bmatrix} = \begin{bmatrix} \frac{2}{m}\sum_{i=1}^{m}(w_{0} + w_{1}x + w_{2}x^{2} - y) \\ \frac{2}{m}\sum_{i=1}^{m}x(w_{0} + w_{1}x + w_{2}x^{2} - y) \\ \frac{2}{m}\sum_{i=1}^{m}x^{2}(w_{0} + w_{1}x + w_{2}x^{2} - y) \end{bmatrix} \tag{10}$∇w​L(w)=⎣⎢⎡​∂w0​∂L(w)​∂w1​∂L(w)​∂w2​∂L(w)​​⎦⎥⎤​=⎣⎡​m2​∑i=1m​(w0​+w1​x+w2​x2−y)m2​∑i=1m​x(w0​+w1​x+w2​x2−y)m2​∑i=1m​x2(w0​+w1​x+w2​x2−y)​⎦⎤​(10)
依据上边一些推导结果，就可以进行demo代码的编写了。当然了，实际上，这里边的求导直接手算给出了结果，在实际应用中，一般要么使用符号微分方法，要么直接使用自动微分进行推导，这个不在这里阐述。代码如下所示：

# -*- coding: utf-8 -*-

import numpy as np


LEARNING_RATE = 0.2
TOLERANCE = 0.00001


def compute_grad(w, x, y):
    grad = [0, 0, 0]
    diff = w[0] + w[1] * x + w[2] * x * x - y
    grad[0] = 2.0 * np.mean(diff)
    x_diff = x * diff
    grad[1] = 2.0 * np.mean(x_diff)
    x_x_diff = np.square(x) * diff
    grad[2] = 2.0 * np.mean(x_x_diff)
    return np.array(grad)


def update_weight(w, learning_rate, grad):
    new_w = np.array(w) - learning_rate * grad
    return new_w


def loss_err(w, x, y):
    loss = np.sqrt(np.mean(np.square(w[0] + w[1] * x + w[2] * np.square(x) - y)))
    return loss


def estimate_weight(x, y):
    w = np.array([1, 1, 1])
    grad = compute_grad(w, x, y)
    loss = loss_err(w, x, y)
    w = update_weight(w, LEARNING_RATE, grad)
    loss_new = loss_err(w, x, y)

    i = 1
    while np.abs(loss_new) > TOLERANCE:
        grad = compute_grad(w, x, y)
        w = update_weight(w, LEARNING_RATE, grad)
        loss = loss_new
        loss_new = loss_err(w, x, y)
        i += 1
        if i % 5 == 0:
            print("Iterate %s: loss error is %s" % (i, loss_new))

    print("Final loss error %s after %s iterations" % (loss_new, i))

    return w


def generate_sample():
    x_list = [i for i in range(5000)]
    max_x = max(x_list)  # normalized to one to prevent overflow
    x_list = [x / max_x for x in x_list]
    y_list = [1 + 2 * x + x * x for x in x_list]

    for y in y_list:
        y += np.random.random(1) * 0.1  # add some noise

    return np.array(x_list), np.array(y_list)


if __name__ == '__main__':
    x_data, y_data = generate_sample()
    weight = estimate_weight(x_data, y_data)
    print(weight)

上边的代码没有什么好说的，只不过，在处理输入数据的时候，必须要归一化，否则会溢出，这是非常关键的。在实际应用中，我想也需要这么做。输出的结果：

...
Iterate 7615: loss error is 1.0201223917743257e-05
Iterate 7620: loss error is 1.0146469766079486e-05
Iterate 7625: loss error is 1.009200950241291e-05
Iterate 7630: loss error is 1.0037841549308142e-05
Final loss error 9.99471659581732e-06 after 7634 iterations
[1.0000246  1.99986244 1.00013272]

@fsfzp888
2018 年 05月 20日