(十)实例二：使用梯度下降求解最小二乘

线性代数 机器学习基础

---

简介

在前一篇的总结中，从矩阵的视角回顾了最小二乘，最终还得到了如下求解最小二乘的解的方程：

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

梯度下降算法步骤说明

在进一步之前，我们先换一种表达方式，前边等式的表达方式换成下边的形式：

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

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

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

代码示例

这篇总结文档中的示例仅仅只是作为演示，去实现最小二乘拟合一些伪造数据，也不具备通用性，只是起到示意作用。
首先，来定义模型，假定我们想要使用下边的方程去拟合数据：

$$y = w_{0} + w_{1}x + w_{2}x^{2} \tag{7}$$ 也可以写作：
$$y = Xw \tag{8}$$ 那么损失函数就是：
$$L(w) = \frac{1}{m}\sum_{i=1}^{m}(w_{0} + w_{1}x + w_{2}x^{2} - y)^{2} \tag{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}$$
依据上边一些推导结果，就可以进行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日