第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维向量。观察这个假设函数的各项,联想到矩阵/向量乘法:
- 设:\(θ=\begin{bmatrix}θ_0\\ θ_1\\ \vdots \\ θ_n\end{bmatrix}, x=\begin{bmatrix}x_0\\ x_1\\ \vdots \\ x_n\end{bmatrix}\);
- 因为:\(θ^T \cdot x = θ_0 x_0 + θ_1 x_1 + θ_2 x_2 \dots + θ_n x_n\);
- 所以:\(h_θ(x)=θ^T \cdot x\)。这是将算法向量化的表述,在实践时十分重要;
- 另外,如果两个矩阵相乘后会得到一个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
。然而这是如何得到这条公式的呢……
- 由于:\(h_θ(x)=\vec{θ}^T \cdot \vec{x} = \vec{x}^T \cdot \vec{θ}\);
- 所以:\(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}\);
- 又由于向量有这个等值运算:\(\sum \vec{v}_i^2 = \vec{v}^T \cdot \vec{v}\),其实就是向量v与自己的内积公式;
- 所以代价函数便可以改写为:\(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})\);
- 由于\(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})]\);
- 又由于\(X \vec{θ}\)是一个m维向量,所以\(X \vec{θ}
= X \cdot \vec{θ}\),代价函数进一步改写为\(J(\vec{θ}) = tr[\frac{1}{2m} (X \vec{θ} – \vec{y})^T (X \vec{θ} – \vec{y})]\); - 线性回归的情况下,我们直接可以认为函数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}
\); - 这里用到了矩阵求迹、求导等技巧,详细请查阅矩阵求导文章。
继续展开方程:
\(
\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 更新