第1周的课程介绍了机器学习的基本概念、单变量的线性回归及其梯度下降算法。接下来要把单变量的线性回归梯度下降算法推广到多重变量，然后介绍如何将算法向量化和Octave的基本操作。

多重变量的线性回归

---

多重变量/特征

对于单变量的线性回归，我们有这个假设函数：\(h_θ(x)=θ_0 + θ_1 x\)。

继续以预测房价为例子，单变量的线性回归以房子的面积为x。为了丰富预测的场景，我们加入了其他额外的因变量，例如房间数量、楼层、楼龄，于是便有面积为\(x_1\)、房间数量为\(x_2\)、楼层为\(x_3\)、楼龄为\(x_4\)、楼价为\(y\)，于是预测房价便是关于因变量向量\(x\)和果变量\(y\)的关系。设n为特征的数量，这里n=4。\(x^{(i)}\)是指第i行训练样本，\(x_j^{(i)}\)是指第i行训练样本的第j个特征：\(x = (x_1, x_2, x_3, x_4) \to y\)。

于是多重变量的假设函数就变为：\(h_θ(x)=θ_0 + θ_1 x_1 + θ_2 x_2 \dots + θ_n x_n\)。

为了追求形式的一致，我们可以设所有训练样本的\(x_0^{(i)} = 1\)，于是假设函数便会变为：\(h_θ(x)=θ_0 x_0 + θ_1 x_1 + θ_2 x_2 \dots + θ_n x_n\)。于是向量x便从一个n维向量变为一个n+1维向量。观察这个假设函数的各项，联想到矩阵/向量乘法：

1. 设：\(θ=\begin{bmatrix}θ_0\\ θ_1\\ \vdots \\ θ_n\end{bmatrix}, x=\begin{bmatrix}x_0\\ x_1\\ \vdots \\ x_n\end{bmatrix}\)；
2. 因为：\(θ^T \cdot x = θ_0 x_0 + θ_1 x_1 + θ_2 x_2 \dots + θ_n x_n\)；
3. 所以：\(h_θ(x)=θ^T \cdot x\)。这是将算法向量化的表述，在实践时十分重要；
4. 另外，如果两个矩阵相乘后会得到一个1*1的矩阵，那么这个1*1矩阵的值其实就是两个矩阵的内积，即：\(h_θ(x)= θ^T \cdot x = θ^T x\)。因为\(θ^T\)的尺寸为1*(n+1)、\(x\)的尺寸为(n+1)*1，相乘后便是尺寸为1*1。

多重变量的梯度下降算法

于是代价函数J便推广为：\(J(θ_0, θ_1, θ_2, \dots, θ_n) = J(θ) = \frac{1}{2m}\sum\limits_{i=1}^m[h_θ(x^{(i)}) – y^{(i)}]^2\)。

同时，梯度下降算法则推广为：

重复 至 收敛 {
    \(θ_j := θ_j – α \frac{∂}{∂ θ_j} J(θ) \) （从j=0至j=n同步更新）
}

重温一下，单变量的梯度下降算法为：

重复 至 收敛 {
    \(θ_0 := θ_0 – α \frac{1}{m} \sum\limits_{i=1}^m [h_θ(x^{(i)}) – y^{(i)}]\)
    \(θ_1 := θ_1 – α \frac{1}{m} \sum\limits_{i=1}^m [h_θ(x^{(i)}) – y^{(i)}] \cdot x^{(i)}\)
}

将\(x_0 = 1\)代入单变量算法，重新展开，便会得到：

重复 至 收敛 {
    \(θ_0 := θ_0 – α \frac{1}{m} \sum\limits_{i=1}^m [h_θ(x^{(i)}) – y^{(i)}] \cdot x_0^{(i)}\)
    \(θ_1 := θ_1 – α \frac{1}{m} \sum\limits_{i=1}^m [h_θ(x^{(i)}) – y^{(i)}] \cdot x_1^{(i)}\)
}

于是推广至j=n，便有多重变量的梯度下降算法：

重复 至 收敛 {
    \(θ_0 := θ_0 – α \frac{1}{m} \sum\limits_{i=1}^m [h_θ(x^{(i)}) – y^{(i)}] \cdot x_0^{(i)}\)
    \(θ_1 := θ_1 – α \frac{1}{m} \sum\limits_{i=1}^m [h_θ(x^{(i)}) – y^{(i)}] \cdot x_1^{(i)}\)
    \(\vdots\)
    \(θ_n := θ_n – α \frac{1}{m} \sum\limits_{i=1}^m [h_θ(x^{(i)}) – y^{(i)}] \cdot x_n^{(i)}\)
}

梯度下降实践1：特征值缩放

在进行梯度下降算法前，要考虑样本特征数据的值范围。继续以预测房价为例，房子的面积可能是从0至2000平方尺，房间的数量则可能是1至5间，这两个特征之间数值范围相差甚大，这会导致特征间与代价函数J构成的那个“碗”的形状十分“瘦/尖”，而两个特征如果数值范围相约，“碗”的形状则十分平缓：

如果“碗”的形状十分平缓，梯度下降算法会下降（收敛）得比较快。但如果“碗”的形状十分“瘦/尖”，则有可能会在下降的时候出现左右“徘徊”的情况，即类似于α值太大的情况，这会导致算法很难收敛、执行得比较慢，甚至会错过最优解，即便调小了α值使其能收敛，同样会因为这个α值太小而执行得很慢。

于是便有将特征的值缩放的需要。

缩放特征有两种计算方法，一种是简易的方法，就是将每个特征值与该特征所有值的期望值（均值μ）相减，再除以值的范围：\(x_1 := \frac{x_1 – μ_1}{S_1}, μ_1 = \frac{\sum_{i=1}^m(x_1^{(i)})}{m}, S_1 = max(x_1) – min(x_1)\)。

第二种方法是更科学、符合统计学意义的计算方法，就是将上式的分母改为计算特征值的标准差：\(S_1 = σ_1 = std(x_1) = \sqrt{\frac{1}{m} \sum_{i=1}^m (x_1^{(i)} – μ_1)^2}\)。如果你的特征值本身服从正态分布，这样缩放过后，缩放后的值便相当于特征值本身偏离期望值有多少个标准差单位，这是有概率意义的。

特征值缩放其实还可以用其他方式进行，只要有理由的都可以进行。吴教授总结了特征值缩放后的范围，最好落在\(-1 \leqslant x_i \leqslant 1\)这个范围。不过实践时放宽到\([-3, 3]\)或者\([-\frac{1}{3}, \frac{1}{3}]\)都是能使梯度下降较快收敛的。

另外，由于\(x_0 = 1\)是人为加进去的，所以不需要对它进行缩放。

梯度下降实践2：学习速率α

剩下来就差学习速率α这个值不知道怎么取值了。

首先假想一下，现在已经有一个恰当的α值，能够使得θ值收敛，那么，代价函数J在通过样本数据执行训练时，代价值的最小值与迭代的次数应该有如下图的关系：

这里收敛的定义，是指J(θ)在某次迭代的时候减少不超过\(10^{-3}\)，也就是再迭代下去，J(θ)的最小值都不会有超过0.001的变化，就没必要在迭代了，相当于收敛了。

如果α值太大，J(θ)便会发散，如下图：

如果α值太小，J(θ)便会来来回回窄幅变化：

吴教授总结了经验，α初值可以从0.001、0.01、0.1、1之间选择，如果初值显得太小，则乘以3，例如0.001显得太小，就选0.003，如果0.003也显得小，就选0.01。反过来如果显得太大，也是类似地除以3进行选择便可以了。

所以要用好梯度下降算法，需要不停地用各种数据图来辅助调整参数，这样才靠谱。

特征与多项式回归

在现实的状况中，有时假设函数是需要仔细考虑的。继续以预测房价为例，有些房子的价格预测可能并不是简单的线性假设就可以满足，例如房子的面积可以拆开为沿路边宽度（frontage）*地皮深度（depth）:

那么假设函数便可以变为：\(h_θ(frontage, depth) = θ_0 + θ_1 frontage + θ_2 depth\)。

其实，在真实数据面前，还可以不妨大胆假设，假设函数便可以相应变为以下各种情况：

虚线即将假设定为二次多项式，可以看到，随着面积递增，房价反而会下降；实线即将假设定为三次多项式，但是后续面积的房价便会升得太快。

这种多次多项式的假设，由于x与\(x_2 = x^2\)、\(x_3 = x^3\)的值范围相差甚大，所以把它们引入作为新的特征时，特征缩放这一步工作就显得尤为重要。

然后又可以将假设函数多项式如下设定：

显然，这个假设是与真实数据拟合得很好的。所以在做回归预测时，不妨大胆假设一下。

正规方程

正规方程是另一种求得θ参数向量的算法，它不需要设定α值，可以十分方便地向量化计算出θ参数向量。

对于m个样本n维特征x和要预测的y值，我们可以设所有\(x_0^{(i)} = 1\)，于是有：
\(\vec{x}^{(i)} = \begin{bmatrix}x_0^{(i)} \\ x_1^{(i)} \\ \vdots \\ x_n^{(i)}\end{bmatrix}, X = \begin{bmatrix}
\vec{x}^{(1)T} \\ \vec{x}^{(2)T} \\ \vdots \\ \vec{x}^{(m)T} \end{bmatrix} = \begin{bmatrix}1 && x_1^{(1)} && x_2^{(1)} && \dots && x_n^{(1)}\\ 1 && x_1^{(2)} && x_2^{(2)} && \dots && x_n^{(2)}\\ \vdots \\ 1 && x_1^{(m)} && x_2^{(m)} && \dots && x_n^{(m)} \end{bmatrix}, \vec{y} =
\begin{bmatrix}y^{(1)} \\ y^{(2)} \\ \vdots \\ y^{(m)}\end{bmatrix}, \vec{θ} \in \mathbb{R}^{n+1,1}\)

我们称X为设计矩阵，它的尺寸为m*(n+1)。

吴教授直接给出了如何计算θ参数向量的正规方程：\(\vec{θ} = (X^T X)^{-1} X^T \vec{y}\)。还直接给出了Octave的代码：pinv(X' * X) * X' * y。然而这是如何得到这条公式的呢……

1. 由于：\(h_θ(x)=\vec{θ}^T \cdot \vec{x} = \vec{x}^T \cdot \vec{θ}\)；
2. 所以：\(X \cdot \vec{θ} – y = \begin{bmatrix} \vec{x}^{(1)T} \cdot \vec{θ} \\ \vec{x}^{(2)T} \cdot \vec{θ} \\ \vdots \\ \vec{x}^{(m)T} \cdot \vec{θ} \end{bmatrix} – \begin{bmatrix}y^{(1)} \\ y^{(2)} \\ \vdots \\ y^{(m)}\end{bmatrix} = \begin{bmatrix} h_θ(x^{(1)}) – y^{(1)} \\ h_θ(x^{(2)}) – y^{(2)} \\ \vdots \\ h_θ(x^{(m)}) – y^{(m)} \end{bmatrix}\)；
3. 又由于向量有这个等值运算：\(\sum \vec{v}_i^2 = \vec{v}^T \cdot \vec{v}\)，其实就是向量v与自己的内积公式；
4. 所以代价函数便可以改写为：\(J(\vec{θ}) = \frac{1}{2m}\sum\limits_{i=1}^m[h_θ(x^{(i)}) – y^{(i)}]^2 = \frac{1}{2m} (X \cdot \vec{θ} – \vec{y})^T \cdot (X \cdot \vec{θ} – \vec{y})\)；
5. 由于\(X \cdot \vec{θ} – \vec{y}\)实际上是一个m维向量，所以\((X \cdot \vec{θ} – \vec{y})^T (X \cdot \vec{θ} – \vec{y})\)是一个1阶方阵，即唯一元素便是它们的内积，即：\(\frac{1}{2m} (X \cdot \vec{θ} – \vec{y})^T \cdot (X \cdot \vec{θ} – \vec{y}) = tr[\frac{1}{2m} (X \cdot \vec{θ} – \vec{y})^T (X \cdot \vec{θ} – \vec{y})]\)；
6. 又由于\(X \vec{θ}\)是一个m维向量，所以\(X \vec{θ}
   = X \cdot \vec{θ}\)，代价函数进一步改写为\(J(\vec{θ}) = tr[\frac{1}{2m} (X \vec{θ} – \vec{y})^T (X \vec{θ} – \vec{y})]\)；
7. 线性回归的情况下，我们直接可以认为函数J的导数等于0处便会找到参数向量θ的解，即解这个方程：\(\begin{aligned}
   \frac{∂}{∂\vec{θ}} J(\vec{θ}) &= \triangledown_\vec{θ} J(\vec{θ}) \\
   &= \triangledown_\vec{θ} tr[\frac{1}{2m} (X\vec{θ} – \vec{y})^T (X\vec{θ} – \vec{y})] \\
   &= \frac{1}{2m} \triangledown_\vec{θ} tr[(\vec{θ}^T X^T – \vec{y}^T) (X\vec{θ} – \vec{y})] \\
   &= \frac{1}{2m} \triangledown_\vec{θ} tr(\vec{θ}^T X^T X \vec{θ} – \vec{y}^T X \vec{θ} – \vec{θ}^T X^T \vec{y} + \vec{y}^T \vec{y}) \\
   &= 0
   \end{aligned}
   \)；
8. 这里用到了矩阵求迹、求导等技巧，详细请查阅矩阵求导文章。

继续展开方程：
\(
\begin{aligned}
\triangledown_\vec{θ} J(\vec{θ}) &= \frac{1}{2m} \triangledown_\vec{θ} tr(\vec{θ}^T X^T X \vec{θ} – \vec{y}^T X \vec{θ} – \vec{θ}^T X^T \vec{y} + \vec{y}^T \vec{y}) \\
&= \frac{1}{2m} [\triangledown_\vec{θ} tr(\vec{θ}^T X^T X \vec{θ}) – \triangledown_\vec{θ} tr(\vec{y}^T X \vec{θ}) – \triangledown_\vec{θ} tr(\vec{θ}^T X^T \vec{y}) + \triangledown_\vec{θ} tr(\vec{y}^T \vec{y})] \\
&= \frac{1}{2m} [(2 X^T X \vec{θ}) – (\vec{y}^T X)^T – X^T \vec{y} + 0] \\
&= \frac{1}{2m}(2 X^T X \vec{θ} – 2 X^T \vec{y}) \\
&= \frac{1}{m}(X^T X \vec{θ} – X^T \vec{y}) = 0
\end{aligned}
\)

由于\(m \gt 0\)，所以即是要求\(X^T X \vec{θ} – X^T \vec{y} = 0\)，即\(X^T X \vec{θ} = X^T \vec{y}\)，最后便有\((X^T X)^{-1} (X^T X) \vec{θ} = I \vec{θ} = \vec{θ} = (X^T X)^{-1} X^T \vec{y}\)。

正规方程虽然能够一下子求得整个θ向量，但是也不是所有线性回归的场景都能用正规方程，它与梯度下降算法相比有如下优缺点：

	梯度下降	正规方程
优点	即使特征的数量n很大，也可以正常运行	不必考虑α值；不必迭代、一次性可计算出结果
缺点	要考虑α值；要多次迭代才可以收敛	需要计算\((X^T X)^{-1}\)，未必有解；当特征的数量n很大时，消耗极大量的内存

不可逆的情况

当矩阵\(X^T X\)不可逆时（奇异矩阵），\((X^T X)^{-1}\)理论上便无法求解，原则上正规方程便不可用于计算θ向量。然而，即使是奇异矩阵，其实也可以通过求得其广义逆矩阵、摩尔-彭诺斯广义逆矩阵获得一个近似的逆矩阵。

对于Matlab或者Octave，便是使用pinv代替inv。

吴教授总结称，如果特征间有线性依赖，例如面积在平方尺与平方米（\(1平方米 = 3.28^2 平方尺\)）同时出现在特征向量，这种情况相当于有冗余的特征，这会导致矩阵\(X^T X\)不可逆，这种情况需要删除多余特征。另外，如果训练样本数量m小于特征数量n，即\(m \leqslant n\)，那么说明特征太多了，也会导致矩阵\(X^T X\)不可逆，也需要删除特征，或者对特征进行正则化（这个会在后面的课程详细说明）。

 

Octave入门

---

Octave的语法貌似和Matlab是一样的，都是比较容易学、面向矩阵运算的语言。只需要记得矩阵、向量的下标是从1开始的便可以了；一个.m文件只能定义一个函数，当前工作目录下的所有.m文件都会被自动搜索、加载；plot等绘图函数用于绘制数据的图表十分有用；其他都是可以随时查随时学的，不赘述了。

 

算法向量化

---

将原本通过循环迭代的算法向量化，是一种数学技巧。有些编程语言自己本身就带有线性代数运算语法，而有些则可以通过常见的第三方库实现。利用线性代数将算法向量化，能够在有特殊系统（例如BLAS、CUDA）或者硬件（例如GPU、APU）支持的情况下，快速、批量地进行原本需要循环迭代才可以进行的运算。将算法向量化，如果算法本身能够通过数学形式从一维、二维推广至多维，那么向量化的算法还能够一次性编写好从一维推广至多维的情况，代码还能够保持简洁、方便维护、减少错误。

以线性回归中的梯度下降算法为例，存在多处可以改写为向量化运算的地方。下面以JavaScript写循环迭代算法、以Octave写向量化算法。

多维特征假设函数\(h_θ(x) = \sum\limits_{j=0}^{n} θ_j x_j = \vec{θ}^T \vec{x}, \vec{θ} \in \mathbb{R}^{n+1,1}, \vec{x} \in \mathbb{R}^{n+1,1}\)：

 // js循环迭代版
function hypothesis(theta, x) {
    const n = x.length
    let sum = 0
    for (let j = 0; j < n; ++j)
        sum += theta[j] * x[j]
    return sum
}

 % Octave向量化版
function h = hypothesis(theta, x)
    h = theta' * x;
end

多维代价函数\(J(θ) = \frac{1}{2m}\sum\limits_{i=1}^m[h_θ(x^{(i)}) – y^{(i)}]^2 = \frac{1}{2m} (X \vec{θ} – \vec{y})^T (X \vec{θ} – \vec{y})\)：

 // js循环迭代版
function computeCost(X, y, theta) {
    let sum = 0
    const m = y.length
    for (let i = 0; i < m; ++i) {
        const err = hypothesis(theta, X[i]) - y[i]
        const sqrErr = err * err
        sum += sqrErr
    }
    return sum / 2 / m
}

 % Octave向量化版
function hs = hypothesises(theta, X)
    hs = X * theta
end

function errs = hypothesisErrors(X, y, theta)
    errs = hypothesises(theta, X) - y
end

function J = computeCost(X, y, theta)
    J = sum(hypothesisErrors(X, y, theta).^2) / 2 / m
end

多重变量梯度下降算法：

 // js循环迭代版
function gradientDescent(X, y, theta, alpha, iters) {
    const m = y.length
    const n = X[0].length
    const J_history = [] // 用以记录每次迭代之后的代价值，方便观察是否收敛
    const coefficient = alpha/m
    for (let itr = 0; itr < iters; ++itr) {
        const hypothesisErrors = []
        for (let i = 0; i < m; ++i)
            hypothesisErrors.push(hypothesis(theta, X[i]) - y[i])
        const temp = []
        for (let j = 0; j < n; ++j) {
            let sum = 0
            for (let i = 0; i < m; ++i)
                sum += hypothesisErrors[i] * x[j]
            temp.push(theta[j] - coefficient * sum)
        }
        for (let j = 0; j < n; ++j)
            theta[j] = temp[j] // 批量、同步更新
        J_history.push(computeCost(X, y, theta))
    }
    return [theta, J_history]
}

先将多重变量梯度下降算法向量化：

重复 至 收敛 {
    \(\begin{bmatrix}θ_0\\ θ_1\\ \vdots \\ θ_n\end{bmatrix} := \begin{bmatrix}θ_0\\ θ_1\\ \vdots \\ θ_n\end{bmatrix} – α \frac{1}{m} \begin{bmatrix}x_0^T\\ x_1^T\\ \vdots \\ x_n^T\end{bmatrix} \begin{bmatrix} h_θ(x^{(1)}) – y^{(1)} \\ h_θ(x^{(2)}) – y^{(2)} \\ \vdots \\ h_θ(x^{(m)}) – y^{(m)} \end{bmatrix}\)
}

即：

重复 至 收敛 {
    \(\vec{θ} := \vec{θ} – α \frac{1}{m} X^T (X \vec{θ} – \vec{y})\)
}

 % Octave向量化版
function [theta, J_history] = gradientDescent(X, y, theta, alpha, iters)
    m = length(y);
    J_history = zeros(iters, 1);
    for itr = 1:iters
        theta = theta - (alpha/m) * X' * hypothesisErrors(X, y, theta); % 批量、同步更新
        J_history(itr) = computeCost(X, y, theta);
    end
end

这样对比就很直白了，不向量化的程序代码需要重新安排一些运算的位置以防性能低下，但也会导致代码比较难以阅读。向量化的程序基本上就是直接将数学公式表达出来，十分易于维护——除非数学基础太差、看不懂。

打赏作者

由 bruce 于 2020-06-8 19:36 更新