我的机器学习主线「优化方法」

我的机器学习主线「优化方法」,第1张

文章目录
    • 优化算法
      • 凸性
        • 凸集
        • 凸函数
          • 凸函数的局部极小值
          • 凸函数的下水平集
          • 凸性和二阶导数
        • 詹森不等式
        • 约束
          • 拉格朗日函数
          • 惩罚
          • 投影
        • 小结
      • 梯度下降算法
        • 一维梯度下降
        • 多元梯度下降
        • 自适应方法
          • 牛顿法
          • 收敛性分析
          • 梯度预处理
          • 线搜索
        • 随机梯度下降
          • 随机梯度更新
          • 收敛性分析
        • 小批量随机梯度下降
          • 向量化和缓存
          • 小批量
        • 小结
      • 启发式算法
        • 退火算法
        • 粒子群算法
        • 遗传算法
        • 蚁群算法

优化算法

机器学习中通常会先定义损失函数,一旦我们有了损失函数,就可以使用优化算法来尝试最小化损失。在优化中,损失函数通常被称为优化问题的目标函数

优化提供了一种最大限度地减少机器学习损失函数的方法,但实质上两者的目标是不同的,前者主要关注的是最小化目标,后者则关注在给定有限数据量的情况下寻找合适的模型。例如,训练误差和泛化误差通常不同,由于优化算法的目标函数通常是基于训练数据集的损失函数,因此优化的目标是减少训练误差。但是机器学习(或更广义的统计推断)的目标是减少泛化误差,为了实现后者,除了使用优化算法来减少训练误差之外还需要注意过拟合

  • 优化有许多挑战

    例如,局部最小值,对于任何目标函数 f ( x ) f(x) f(x),如果在 x x x 处对应的 f ( x ) f(x) f(x) 值小于在 x x x 附近任何其他点的 f ( x ) f(x) f(x) 值,那么 f ( x ) f(x) f(x) 可能是局部最小值。如果 f ( x ) f(x) f(x) x x x 处的值是整个域上目标函数的最小值,那么 f ( x ) f(x) f(x) 是全局最小值,目标函数通常有许多局部最优解。当优化问题的数值解接近局部最优值时,随着目标函数解的梯度接近或变为零,通过最终迭代获得的数值解可能仅使目标函数局部最优,而不是全局最优

    例如,鞍点,指函数的所有梯度都消失但既不是全局最小值也不是局部最小值的任何位置。对于函数 f ( x ) = x 3 f(x)=x3 f(x)=x3。它的一阶和二阶导数在 x = 0 x=0 x=0 时消失,此时优化可能会停止,但不是最小值,对于高维度问题,当函数在零梯度位置处的 H e s s i a n Hessian Hessian 矩阵的非正定时即为鞍点,这使得鞍点比局部最小值更有可能

    例如,梯度消失,最小化函数 f ( x ) = t a n h ⁡ ( x ) f(x)=tanh⁡(x) f(x)=tanh(x),若初值为 x = 4 x=4 x=4 f f f 的梯度接近零, f ′ ( 4 ) = 1 − t a n h 2 ⁡ ( 4 ) = 0.0013 f^′(4)=1−tanh^2⁡(4)=0.0013 f(4)=1tanh2(4)=0.0013 优化将会停滞很长一段时间

凸性

凸性在优化算法的设计中起到至关重要的作用,主要是由于在这种情况下对算法进行分析和测试要容易得多。换言之,如果该算法甚至在凸性条件设定下的效果很差,通常我们很难在其他条件下看到好的结果。此外,即使机器学习中的优化问题通常是非凸的,它们也经常在局部极小值附近表现出一些凸性

凸集

凸集是凸性的基础。简单地说,如果对于任何 a , b ∈ X a,b \in \mathcal{X} abX,连接 a a a b b b 的线段也位于 X \mathcal{X} X 中,则向量空间中的一个集合 X \mathcal{X} X 是 凸的。在数学术语上,这意味着对于所有 λ ∈ [ 0 , 1 ] \lambda \in [0, 1] λ[0,1],我们得到
λ a + ( 1 − λ ) b ∈ X  当  a , b ∈ X \lambda a + (1-\lambda) b \in \mathcal{X} \text{ 当 } a, b \in \mathcal{X} λa+(1λ)bX  a,bX

下第一组存在不包含在集合内部的线段,所以该集合是非凸的,而下另外两组则没有这样的问题,所以是凸的
假设下第一组 X \mathcal{X} X 和下第二组 Y \mathcal{Y} Y 是凸集,那么 X ∩ Y \mathcal {X} \cap \mathcal{Y} XY 也是凸集的。现在考虑任意 a , b ∈ X ∩ Y a, b \in \mathcal{X} \cap \mathcal{Y} a,bXY,因为 X \mathcal{X} X Y \mathcal{Y} Y 是凸集,所以连接 a a a b b b 的线段包含在 X \mathcal{X} X Y \mathcal{Y} Y 中。鉴于此,它们也需要包含在 X ∩ Y \mathcal {X} \cap \mathcal{Y} XY 中,从而证明的定理

可以毫不费力地进一步得到,给定凸集 X i \mathcal{X}_i Xi,它们的交集 ∩ i X i \cap_{i} \mathcal{X}_i iXi 是凸的。但是反向是不正确的,考虑两个不相交的集合 X ∩ Y = ∅ \mathcal{X} \cap \mathcal{Y} = \emptyset XY=,取 a ∈ X a \in \mathcal{X} aX b ∈ Y b \in \mathcal{Y} bY。假设 X ∩ Y = ∅ \mathcal{X} \cap \mathcal{Y} = \emptyset XY=,连接 a a a b b b 的线段需要包含一部分既不在 X \mathcal{X} X 也不在 Y \mathcal{Y} Y 中。因此线段也不在 X ∪ Y \mathcal{X} \cup \mathcal{Y} XY 中,因此证明了凸集的并集不一定是凸的,即非凸的

通常机器学习中的问题是在凸集上定义的。例如, R d \mathbb{R}^d Rd 即实数的 d d d 维向量的集合是凸集(毕竟 R d \mathbb{R}^d Rd 中任意两点之间的线存在 R d \mathbb{R}^d Rd)中。在某些情况下使用有界长度的变量,例如球的半径定义为 { x ∣ x ∈ R d  and  ∥ x ∥ ≤ r } \{\mathbf{x} | \mathbf{x} \in \mathbb{R}^d \text{ and } \| \mathbf{x} \| \leq r\} {xxRd and xr}

凸函数

有了凸集就可以引入凸函数 f f f。给定一个凸集 X \mathcal{X} X,如果对于所有 x , x ′ ∈ X x, x' \in \mathcal{X} x,xX 和所有 λ ∈ [ 0 , 1 ] \lambda \in [0, 1] λ[0,1],一个函数 f : X → R f: \mathcal{X} \to \mathbb{R} f:XR 是凸的,可以得到
λ f ( x ) + ( 1 − λ ) f ( x ′ ) ≥ f ( λ x + ( 1 − λ ) x ′ ) \lambda f(x) + (1-\lambda) f(x') \geq f(\lambda x + (1-\lambda) x') λf(x)+(1λ)f(x)f(λx+(1λ)x)
为了说明这一点,下面定义一些函数,包括凸函数和非凸函数
f ( x ) = 0.5 ∗ x 2 凸函数 g ( x ) = c o s ( π x ) 非凸函数 h ( x ) = exp ⁡ ( 0.5 ∗ x ) 凸函数 \begin{aligned} f(x)&=0.5\ast x^2 &\text{凸函数}\ g(x)&=cos(\pi x)&\text{非凸函数}\ h(x)&=\exp(0.5\ast x)&\text{凸函数} \end{aligned} f(x)g(x)h(x)=0.5x2=cos(πx)=exp(0.5x)凸函数非凸函数凸函数

余弦函数为非凸的,而抛物线函数和指数函数为凸的,且 X \mathcal{X} X 是凸集的要求是必要的,否则可能无法很好地界定 f ( λ x + ( 1 − λ ) x ′ ) f(\lambda x + (1-\lambda) x') f(λx+(1λ)x) 的结果

凸函数的局部极小值

假设 x ∗ ∈ X x^{\ast} \in \mathcal{X} xX 是一个局部最小值,则存在一个很小的正值 p p p,使得当 x ∈ X x \in \mathcal{X} xX 满足 0 < ∣ x − x ∗ ∣ ≤ p 0 < |x - x^{\ast}| \leq p 0<xxp 时,有 f ( x ∗ ) < f ( x ) f(x^{\ast}) < f(x) f(x)<f(x)

现在假设局部极小值 x ∗ x^{\ast} x 不是 f f f 的全局极小值,存在 x ′ ∈ X x' \in \mathcal{X} xX 使得 f ( x ′ ) < f ( x ∗ ) f(x') < f(x^{\ast}) f(x)<f(x)。则存在 λ ∈ [ 0 , 1 ) \lambda \in [0, 1) λ[0,1),比如 λ = 1 − p ∣ x ∗ − x ′ ∣ \lambda = 1 - \frac{p}{|x^{\ast} - x'|} λ=1xxp,使得 0 < ∣ λ x ∗ + ( 1 − λ ) x ′ − x ∗ ∣ ≤ p 0 < |\lambda x^{\ast} + (1-\lambda) x' - x^{\ast}| \leq p 0<λx+(1λ)xxp

然而,根据凸性的性质有
f ( λ x ∗ + ( 1 − λ ) x ′ ) ≤ λ f ( x ∗ ) + ( 1 − λ ) f ( x ′ ) < λ f ( x ∗ ) + ( 1 − λ ) f ( x ∗ ) = f ( x ∗ ) \begin{aligned} f(\lambda x^{\ast} + (1-\lambda) x') &\leq \lambda f(x^{\ast}) + (1-\lambda) f(x') \ &< \lambda f(x^{\ast}) + (1-\lambda) f(x^{\ast}) \ &= f(x^{\ast}) \ \end{aligned} f(λx+(1λ)x)λf(x)+(1λ)f(x)<λf(x)+(1λ)f(x)=f(x)
这与 x ∗ x^{\ast} x 是局部最小值相矛盾。因此不存在 x ′ ∈ X x' \in \mathcal{X} xX 满足 f ( x ′ ) < f ( x ∗ ) f(x') < f(x^{\ast}) f(x)<f(x)。综上所述,局部最小值 x ∗ x^{\ast} x 也是全局最小值。例如,对于凸函数 f ( x ) = ( x − 1 ) 2 f(x) = (x-1)^2 f(x)=(x1)2,有一个局部最小值 x = 1 x=1 x=1,它也是全局最小值

凸函数的局部极小值同时也是全局极小值这一性质是很方便的。这意味着如果最小化函数就不会卡住。但是这并不意味着不能有多个全局最小值,或者不存在一个全局最小值。例如,函数 f ( x ) = m a x ( ∣ x ∣ − 1 , 0 ) f(x) = \mathrm{max}(|x|-1, 0) f(x)=max(x1,0) [ − 1 , 1 ] [-1,1] [1,1] 区间上都是最小值。相反,函数 f ( x ) = exp ⁡ ( x ) f(x) = \exp(x) f(x)=exp(x) R \mathbb{R} R 上没有取得最小值。对于 x → − ∞ x \to -\infty x,它趋近于 0 0 0,但是没有 f ( x ) = 0 f(x) = 0 f(x)=0 x x x

凸函数的下水平集

通过凸函数的下水平集定义凸集,给定一个定义在凸集 X \mathcal{X} X 上的凸函数 f f f,其任意一个下水平集
S b : = { x ∣ x ∈ X  and  f ( x ) ≤ b } \mathcal{S}_b := \{x | x \in \mathcal{X} \text{ and } f(x) \leq b\} Sb:={xxX and f(x)b}
是凸的,对于任何 x , x ′ ∈ S b x, x' \in \mathcal{S}_b x,xSb,需要证明,当 λ ∈ [ 0 , 1 ] \lambda \in [0, 1] λ[0,1]时, λ x + ( 1 − λ ) x ′ ∈ S b \lambda x + (1-\lambda) x' \in \mathcal{S}_b λx+(1λ)xSb。因为 f ( x ) ≤ b f(x) \leq b f(x)b f ( x ′ ) ≤ b f(x') \leq b f(x)b,所以
f ( λ x + ( 1 − λ ) x ′ ) ≤ λ f ( x ) + ( 1 − λ ) f ( x ′ ) ≤ b f(\lambda x + (1-\lambda) x') \leq \lambda f(x) + (1-\lambda) f(x') \leq b f(λx+(1λ)x)λf(x)+(1λ)f(x)b

凸性和二阶导数

当一个函数的二阶导数 f : R n → R f: \mathbb{R}^n \rightarrow \mathbb{R} f:RnR 存在时,很容易检查这个函数的凸性。需要做的就是检查 ∇ 2 f ⪰ 0 \nabla^2f \succeq 0 2f0,即对于所有 x ∈ R n \mathbf{x} \in \mathbb{R}^n xRn x ⊤ H x ≥ 0 \mathbf{x}^\top \mathbf{H} \mathbf{x} \geq 0 xHx0。例如,函数 f ( x ) = 1 2 ∥ x ∥ 2 f(\mathbf{x}) = \frac{1}{2} \|\mathbf{x}\|^2 f(x)=21x2 是凸的,因为 ∇ 2 f = 1 \nabla^2 f = \mathbf{1} 2f=1,即其导数是单位矩阵

更严谨的 f f f 为凸函数,当且仅当任意二次可微一维函数 f : R n → R f: \mathbb{R}^n \rightarrow \mathbb{R} f:RnR 是凸的。对于任意二次可微多维函数 f : R n → R f: \mathbb{R}^{n} \rightarrow \mathbb{R} f:RnR,它是凸的当且仅当它的 H e s s i a n Hessian Hessian 矩阵即 ∇ 2 f ⪰ 0 \nabla^2f\succeq 0 2f0

首先,证明一维情况。为了证明凸函数的 f ′ ′ ( x ) ≥ 0 f''(x) \geq 0 f(x)0,使用
1 2 f ( x + ϵ ) + 1 2 f ( x − ϵ ) ≥ f ( x + ϵ 2 + x − ϵ 2 ) = f ( x ) \frac{1}{2} f(x + \epsilon) + \frac{1}{2} f(x - \epsilon) \geq f\left(\frac{x + \epsilon}{2} + \frac{x - \epsilon}{2}\right) = f(x) 21f(x+ϵ)+21f(xϵ)f(2x+ϵ+2xϵ)=f(x)
因为二阶导数是由有限差分的极限给出的,所以遵循
f ′ ′ ( x ) = lim ⁡ ϵ → 0 f ( x + ϵ ) + f ( x − ϵ ) − 2 f ( x ) ϵ 2 ≥ 0 f''(x) = \lim_{\epsilon \to 0} \frac{f(x+\epsilon) + f(x - \epsilon) - 2f(x)}{\epsilon^2} \geq 0 f(x)=ϵ0limϵ2f(x+ϵ)+f(xϵ)2f(x)0
为了证明 f ′ ′ ≥ 0 f'' \geq 0 f0 可以推导 f f f 是凸的,使用一个事实 f ′ ′ ≥ 0 f'' \geq 0 f0 意味着 f ′ f' f 是一个单调的非递减函数。假设 a < x < b a < x < b a<x<b R \mathbb{R} R 中的三个点,其中 x = ( 1 − λ ) a + λ b x = (1-\lambda)a + \lambda b x=(1λ)a+λb λ ∈ ( 0 , 1 ) \lambda \in (0, 1) λ(0,1) 根据中值定理,存在 α ∈ [ a , x ] \alpha \in [a, x] α[a,x] β ∈ [ x , b ] \beta \in [x, b] β[x,b],使得
f ′ ( α ) = f ( x ) − f ( a ) x − a  且  f ′ ( β ) = f ( b ) − f ( x ) b − x f'(\alpha) = \frac{f(x) - f(a)}{x-a} \text{ 且 } f'(\beta) = \frac{f(b) - f(x)}{b-x} f(α)=xaf(x)f(a)  f(β)=bxf(b)f(x)
通过单调性 f ′ ( β ) ≥ f ′ ( α ) f'(\beta) \geq f'(\alpha) f(β)f(α),因此
x − a b − a f ( b ) + b − x b − a f ( a ) ≥ f ( x ) \frac{x-a}{b-a}f(b) + \frac{b-x}{b-a}f(a) \geq f(x) baxaf(b)+babxf(a)f(x)
由于 x = ( 1 − λ ) a + λ b x = (1-\lambda)a + \lambda b x=(1λ)a+λb,所以
λ f ( b ) + ( 1 − λ ) f ( a ) ≥ f ( ( 1 − λ ) a + λ b ) \lambda f(b) + (1-\lambda)f(a) \geq f((1-\lambda)a + \lambda b) λf(b)+(1λ)f(a)f((1λ)a+λb)
从而证明了凸性

其次,需要一个引理证明多维情况 f : R n → R f: \mathbb{R}^n \rightarrow \mathbb{R} f:RnR 是凸的当且仅当对于所有 x , y ∈ R n \mathbf{x}, \mathbf{y} \in \mathbb{R}^n x,yRn
g ( z ) = d e f f ( z x + ( 1 − z ) y )  where  z ∈ [ 0 , 1 ] g(z) \stackrel{\mathrm{def}}{=} f(z \mathbf{x} + (1-z) \mathbf{y}) \text{ where } z \in [0,1] g(z)=deff(zx+(1z)y) where z[0,1]
是凸的

为了证明 f f f 的凸性意味着 g g g 是凸的,可以证明,对于所有的 a , b , λ ∈ [ 0 , 1 ] a,b,\lambda \in[0,1] abλ[01](这样有 0 ≤ λ a + ( 1 − λ ) b ≤ 1 0 \leq \lambda a + (1-\lambda) b \leq 1 0λa+(1λ)b1
g ( λ a + ( 1 − λ ) b ) = f ( ( λ a + ( 1 − λ ) b ) x + ( 1 − λ a − ( 1 − λ ) b ) y ) = f ( λ ( a x + ( 1 − a ) y ) + ( 1 − λ ) ( b x + ( 1 − b ) y ) ) ≤ λ f ( a x + ( 1 − a ) y ) + ( 1 − λ ) f ( b x + ( 1 − b ) y ) = λ g ( a ) + ( 1 − λ ) g ( b ) . \begin{aligned} &g(\lambda a + (1-\lambda) b)\ =&f\left(\left(\lambda a + (1-\lambda) b\right)\mathbf{x} + \left(1-\lambda a - (1-\lambda) b\right)\mathbf{y} \right)\ =&f\left(\lambda \left(a \mathbf{x} + (1-a) \mathbf{y}\right) + (1-\lambda) \left(b \mathbf{x} + (1-b) \mathbf{y}\right) \right)\ \leq& \lambda f\left(a \mathbf{x} + (1-a) \mathbf{y}\right) + (1-\lambda) f\left(b \mathbf{x} + (1-b) \mathbf{y}\right) \ =& \lambda g(a) + (1-\lambda) g(b). \end{aligned} ===g(λa+(1λ)b)f((λa+(1λ)b)x+(1λa(1λ)b)y)f(λ(ax+(1a)y)+(1λ)(bx+(1b)y))λf(ax+(1a)y)+(1λ)f(bx+(1b)y)λg(a)+(1λ)g(b).
为了证明这一点,可以证明对 [ 0 , 1 ] [0,1] [01] 中所有的 λ \lambda λ
f ( λ x + ( 1 − λ ) y ) = g ( λ ⋅ 1 + ( 1 − λ ) ⋅ 0 ) ≤ λ g ( 1 ) + ( 1 − λ ) g ( 0 ) = λ f ( x ) + ( 1 − λ ) f ( y ) . \begin{aligned} &f(\lambda \mathbf{x} + (1-\lambda) \mathbf{y})\ =&g(\lambda \cdot 1 + (1-\lambda) \cdot 0)\ \leq& \lambda g(1) + (1-\lambda) g(0) \ =& \lambda f(\mathbf{x}) + (1-\lambda) f(\mathbf{y}). \end{aligned} ==f(λx+(1λ)y)g(λ1+(1λ)0)λg(1)+(1λ)g(0)λf(x)+(1λ)f(y).
最后,利用上面的引理和一维情况的结果,可以证明多维情况,多维函数 f : R n → R f:\mathbb{R}^n\rightarrow\mathbb{R} f:RnR 是凸函数,当且仅当 g ( z ) = d e f f ( z x + ( 1 − z ) y ) g(z) \stackrel{\mathrm{def}}{=} f(z \mathbf{x} + (1-z) \mathbf{y}) g(z)=deff(zx+(1z)y) 是凸的,这里 z ∈ [ 0 , 1 ] z \in [0,1] z[0,1] x , y ∈ R n \mathbf{x}, \mathbf{y} \in \mathbb{R}^n x,yRn。根据一维情况,此条成立的条件为,当且仅当对于所有 x , y ∈ R n \mathbf{x}, \mathbf{y} \in \mathbb{R}^n x,yRn g ′ ′ = ( x − y ) ⊤ H ( x − y ) ≥ 0 g'' = (\mathbf{x} - \mathbf{y})^\top \mathbf{H}(\mathbf{x} - \mathbf{y}) \geq 0 g=(xy)H(xy)0 H = d e f ∇ 2 f \mathbf{H} \stackrel{\mathrm{def}}{=} \nabla^2f H=def2f)这相当于根据半正定矩阵的定义, H ⪰ 0 \mathbf{H} \succeq 0 H0

詹森不等式

给定一个凸函数 f f f,最有用的数学工具之一就是詹森不等式,是凸性定义的一种推广
∑ i α i f ( x i ) ≥ f ( ∑ i α i x i )  and  E X [ f ( X ) ] ≥ f ( E X [ X ] ) \sum_i \alpha_i f(x_i) \geq f\left(\sum_i \alpha_i x_i\right) \text{ and } E_X[f(X)] \geq f\left(E_X[X]\right) iαif(xi)f(iαixi) and EX[f(X)]f(EX[X])
其中 α i \alpha_i αi 是满足 ∑ i α i = 1 \sum_i \alpha_i = 1 iαi=1 的非负实数, X X X 是随机变量,即凸函数的期望不小于期望的凸函数,其中后者通常是一个更简单的表达式

詹森不等式的一个常见应用,用一个较简单的表达式约束一个较复杂的表达式,可以应用于部分观察到的随机变量的对数似然。具体地,由于 ∫ P ( Y ) P ( X ∣ Y ) d Y = P ( X ) \int P(Y) P(X \mid Y) dY = P(X) P(Y)P(XY)dY=P(X),所以
E Y ∼ P ( Y ) [ − log ⁡ P ( X ∣ Y ) ] ≥ − log ⁡ P ( X ) E_{Y \sim P(Y)}[-\log P(X \mid Y)] \geq -\log P(X) EYP(Y)[logP(XY)]logP(X)
这里 Y Y Y是典型的未观察到的随机变量, P ( Y ) P(Y) P(Y) 是它可能如何分布的最佳猜测, P ( X ) P(X) P(X) 是将 Y Y Y 积分后的分布。例如,在聚类中 Y Y Y 可能是簇标签,而在应用簇标签时, P ( X ∣ Y ) P(X \mid Y) P(XY) 是生成模型

约束

凸优化的一个很好的特性是能够让我们有效地处理约束,即它使我们能够解决以下形式的约束优化问题
m i n i m i z e   x f ( x )  subject to  c i ( x ) ≤ 0  for all  i ∈ { 1 , … , N } . \begin{aligned} \mathop{\mathrm{minimize~}}_{\mathbf{x}} & f(\mathbf{x}) \ \text{ subject to } & c_i(\mathbf{x}) \leq 0 \text{ for all } i \in \{1, \ldots, N\}. \end{aligned} minimize x subject to f(x)ci(x)0 for all i{1,,N}.
这里 f f f 是目标函数, c i c_i ci 是约束函数。例如第一个约束 c 1 ( x ) = ∥ x ∥ 2 − 1 c_1(\mathbf{x}) = \|\mathbf{x}\|_2 - 1 c1(x)=x21,则参数 x \mathbf{x} x 被限制为单位球。如果第二个约束 c 2 ( x ) = v ⊤ x + b c_2(\mathbf{x}) = \mathbf{v}^\top \mathbf{x} + b c2(x)=vx+b,那么这对应于半空间上所有的 x \mathbf{x} x。同时满足这两个约束等于选择一个球的切片作为约束集

拉格朗日函数

求解一个有约束的优化问题是困难的,解决这个问题的一种方法来自物理中相当简单的直觉。想象一个球在一个盒子里,球会滚到最低的地方,重力将与盒子两侧对球施加的力平衡。简而言之,目标函数(即重力)的梯度将被约束函数的梯度所抵消(由于墙壁的推回作用,需要保持在盒子内)任何不起作用的约束(即球不接触壁)都将无法对球施加任何力

使用拉格朗日函数 L L L,上述推理可以通过以下鞍点优化问题来表示
L ( x , α 1 , … , α n ) = f ( x ) + ∑ i = 1 n α i c i ( x )  where  α i ≥ 0 L(\mathbf{x}, \alpha_1, \ldots, \alpha_n) = f(\mathbf{x}) + \sum_{i=1}^n \alpha_i c_i(\mathbf{x}) \text{ where } \alpha_i \geq 0 L(x,α1,,αn)=f(x)+i=1nαici(x) where αi0
这里的变量 α i \alpha_i αi i = 1 , … , n i=1,\ldots,n i=1,,n)是拉格朗日乘数它确保约束被正确地执行。选择它们的大小足以确保所有 i i i c i ( x ) ≤ 0 c_i(\mathbf{x}) \leq 0 ci(x)0。例如,对于 c i ( x ) < 0 c_i(\mathbf{x}) < 0 ci(x)<0 中任意 x \mathbf{x} x,最终会选择 α i = 0 \alpha_i = 0 αi=0。此外,这是一个鞍点优化问题。在这个问题中,想要使 L L L 相对于 α i \alpha_i αi 最大化,同时使它相对于 x \mathbf{x} x 最小化, L L L 的鞍点是原始约束优化问题的最优解

惩罚

一种至少近似地满足约束优化问题的方法是采用拉格朗日函数 L L L。除了满足 c i ( x ) ≤ 0 c_i(\mathbf{x}) \leq 0 ci(x)0 之外,我们只需将 α i c i ( x ) \alpha_i c_i(\mathbf{x}) αici(x) 添加到目标函数 f ( x ) f(x) f(x)。这样可以确保不会严重违反约束

机器学习中一直在使用这个技巧。比如权重衰减,在目标函数中加入 λ 2 ∣ w ∣ 2 \frac{\lambda}{2} |\mathbf{w}|^2 2λw2,以确保 w \mathbf{w} w 不会增长太大。使用约束优化的观点,可以看到,对于若干半径 r r r,这将确保 ∣ w ∣ 2 − r 2 ≤ 0 |\mathbf{w}|^2 - r^2 \leq 0 w2r20。通过调整 λ \lambda λ 的值,可以改变 w \mathbf{w} w 的大小

通常,添加惩罚是确保近似满足约束的一种好方法。在实践中证明虽然牺牲了一些精确但结果更满意更可靠。此外,对于非凸问题,许多使精确方法在凸情况下的性质(例如,可求最优解)是不成立的

投影

满足约束条件的另一种策略是投影,同样,例如在处理梯度截断时,通过
g ← g ⋅ m i n ( 1 , θ / ∥ g ∥ ) \mathbf{g} \leftarrow \mathbf{g} \cdot \mathrm{min}(1, \theta/\|\mathbf{g}\|) ggmin(1,θ/g)
确保梯度的长度以 θ \theta θ 为界限

这就是 g \mathbf{g} g 在半径为 θ \theta θ 的球上的投影,更泛化地说,在凸集 X \mathcal{X} X 上的投影被定义为
P r o j X ( x ) = a r g m i n x ′ ∈ X ∥ x − x ′ ∥ \mathrm{Proj}_\mathcal{X}(\mathbf{x}) = \mathop{\mathrm{argmin}}_{\mathbf{x}' \in \mathcal{X}} \|\mathbf{x} - \mathbf{x}'\| ProjX(x)=argminxXxx
它是 X \mathcal{X} X 中离 X \mathbf{X} X 最近的点

投影的数学定义就如上,两个凸集,一个圆和一个菱形。两个集合内的点(黄色)在投影期间保持不变。两个集合(黑色)之外的点投影到集合中接近原始点(黑色)的点(红色)

凸投影的一个用途是计算稀疏权重向量。例如,将权重向量投影到一个 L 1 L_1 L1 的球(菱形的一个广义版本)上

小结
  • 凸集的交点是凸的,并集不是
  • 一个二次可微函数是凸函数,当且仅当其 H e s s i a n Hessian Hessian 矩阵是半正定的
  • 根据詹森不等式「一个多变量凸函数的总期望值」大于或等于「用每个变量的期望值计算这个函数的总值」
  • 凸约束可以通过拉格朗日函数来实现,只需在目标函数中加上一个惩罚就可以了
  • 投影映射到凸集中最接近原始点的点
梯度下降算法

梯度下降是机器学习的关键优化算法

一维梯度下降

连续可微实值函数 f : R → R f: \mathbb{R} \rightarrow \mathbb{R} f:RR 利用泰勒展开
f ( x + ϵ ) = f ( x ) + ϵ f ′ ( x ) + O ( ϵ 2 ) f(x + \epsilon) = f(x) + \epsilon f'(x) + \mathcal{O}(\epsilon^2) f(x+ϵ)=f(x)+ϵf(x)+O(ϵ2)
一阶近似 f ( x + ϵ ) f(x+\epsilon) f(x+ϵ) 可通过 x x x 处的函数值 f ( x ) f(x) f(x) 和一阶导数 f ′ ( x ) f'(x) f(x) 得出。假设在负梯度方向上移动的 ϵ \epsilon ϵ 会减少 f f f 选择固定步长 η > 0 \eta > 0 η>0 然后取 ϵ = − η f ′ ( x ) \epsilon = -\eta f'(x) ϵ=ηf(x) 将其代入泰勒展开式可以得到
f ( x − η f ′ ( x ) ) = f ( x ) − η f ′ 2 ( x ) + O ( η 2 f ′ 2 ( x ) ) f(x - \eta f'(x)) = f(x) - \eta f'^2(x) + \mathcal{O}(\eta^2 f'^2(x)) f(xηf(x))=f(x)ηf2(x)+O(η2f2(x))
η f ′ 2 ( x ) > 0 \eta f'^2(x)>0 ηf2(x)>0 梯度下降,如果导数 f ′ ( x ) ≠ 0 f'(x) \neq 0 f(x)=0 即梯度没有消失则可继续展开。此外总有 η \eta η 小到足以使高阶项与结果不相关
f ( x − η f ′ ( x ) ) ⪅ f ( x ) f(x - \eta f'(x)) \lessapprox f(x) f(xηf(x))f(x)
意味着,如果使用
x ← x − η f ′ ( x ) x \leftarrow x - \eta f'(x) xxηf(x)
迭代 x x x,函数 f ( x ) f(x) f(x) 的值可能会下降。因此,在梯度下降中,首先选择初始值 x x x 和常数 η > 0 \eta > 0 η>0,然后使用它们连续迭代 x x x,直到梯度 ∣ f ′ ( x ) ∣ |f'(x)| f(x) 的幅度足够小或迭代次数达到规定值

通过目标函数 f ( x ) = x 2 f(x)=x^2 f(x)=x2 实践梯度下降。求导可知 x = 0 x=0 x=0 f ( x ) f(x) f(x) 取最小值,令 x = 10 x=10 x=10 作为初始值,设 η = 0.2 \eta=0.2 η=0.2 迭代 10 10 10 x x x 可以发现 x x x 的值最终将接近最优解

# 目标函数
def f(x):  
    return x ** 2

# 目标函数的梯度
def f_grad(x):  
    return 2 * x
  
# 梯度下降
def gd(eta, f_grad):
    x = 10.0
    results = [x]
    for i in range(10):
        x -= eta * f_grad(x)
        results.append(float(x))
    print(f'epoch 10, x: {x:f}')
    return results

results = gd(0.2, f_grad)


学习率 η \eta η 决定目标函数能否收敛到局部最小值,以及何时收敛到最小值。如果学习率太小,将导致 x x x 的更新非常缓慢,需要更多的迭代。在优化 f ( x ) = x 2 f(x)=x^2 f(x)=x2 问题中 η = 0.05 \eta = 0.05 η=0.05 尽管经过了 10 10 10 次迭代,仍然离最优解很远

相反,如果学习率太大则 ∣ η f ′ ( x ) ∣ \left|\eta f'(x)\right| ηf(x) f ( x ) f(x) f(x) 的一阶泰勒近似而言可能不够小,导致 O ( η 2 f ′ 2 ( x ) ) \mathcal{O}(\eta^2 f'^2(x)) O(η2f2(x)) 不是可以忽略的无穷小项。此时 x x x 的迭代不能保证降低 f ( x ) f(x) f(x) 的值。当学习率为 η = 1.1 \eta=1.1 η=1.1 时, x x x 跨过了最优解 x = 0 x=0 x=0 并逐渐发散

非凸函数的梯度下降,考虑函数 f ( x ) = x ⋅ cos ⁡ ( c x ) f(x) = x \cdot \cos(cx) f(x)=xcos(cx),其中 c c c 为常数。该函数有无穷多个局部最小值。根据选择的学习率,可能只会得到一个局部最小值且不切实际的高学习率可能导致较差的局部最小值

多元梯度下降

考虑 x = [ x 1 , x 2 , … , x d ] ⊤ \mathbf{x} = [x_1, x_2, \ldots, x_d]^\top x=[x1,x2,,xd] 的情况,目标函数 f : R d → R f: \mathbb{R}^d \to \mathbb{R} f:RdR 将向量映射成标量,相应它的梯度也是多元的,是一个由 d d d 个偏导数组成的向量
∇ f ( x ) = [ ∂ f ( x ) ∂ x 1 , ∂ f ( x ) ∂ x 2 , … , ∂ f ( x ) ∂ x d ] ⊤ \nabla f(\mathbf{x}) = \bigg[\frac{\partial f(\mathbf{x})}{\partial x_1}, \frac{\partial f(\mathbf{x})}{\partial x_2}, \ldots, \frac{\partial f(\mathbf{x})}{\partial x_d}\bigg]^\top f(x)=[x1f(x),x2f(x),,xdf(x)]
梯度中的每个偏导数元素 ∂ f ( x ) / ∂ x i \partial f(\mathbf{x})/\partial x_i f(x)/xi 代表了当输入 x i x_i xi f f f x \mathbf{x} x 处的变化率,与单变量的情况一样,可以对多变量函数使用多元泰勒近似
f ( x + ϵ ) = f ( x ) + ϵ ⊤ ∇ f ( x ) + O ( ∣ ϵ ∣ 2 ) f(\mathbf{x} + \boldsymbol{\epsilon}) = f(\mathbf{x}) + \mathbf{\boldsymbol{\epsilon}}^\top \nabla f(\mathbf{x}) + \mathcal{O}(|\boldsymbol{\epsilon}|^2) f(x+ϵ)=f(x)+ϵf(x)+O(ϵ2)
ϵ \boldsymbol{\epsilon} ϵ 的二阶项中,最陡下降的方向由负梯度 − ∇ f ( x ) -\nabla f(\mathbf{x}) f(x) 得出,选择合适的学习率 η > 0 \eta > 0 η>0 实现梯度下降算法
x ← x − η ∇ f ( x ) \mathbf{x} \leftarrow \mathbf{x} - \eta \nabla f(\mathbf{x}) xxηf(x)
构造一个目标函数 f ( x ) = x 1 2 + 2 x 2 2 f(\mathbf{x})=x_1^2+2x_2^2 f(x)=x12+2x22 以二维向量 x = [ x 1 , x 2 ] ⊤ \mathbf{x} = [x_1, x_2]^\top x=[x1,x2] 作为输入,输出标量。梯度由 ∇ f ( x ) = [ 2 x 1 , 4 x 2 ] ⊤ \nabla f(\mathbf{x}) = [2x_1, 4x_2]^\top f(x)=[2x1,4x2] 给出,从初始位置 [ − 5 , − 2 ] [-5, -2] [5,2] 开始使用学习率 η = 0.1 \eta = 0.1 η=0.1 通过梯度下降优化,可以看到经过 20 20 20 步之后, x \mathbf{x} x 的值缓慢接近其最小值 [ 0 , 0 ] [0, 0] [0,0]

# 目标函数
def f_2d(x1, x2):  
    return x1 ** 2 + 2 * x2 ** 2

# 目标函数的梯度
def f_2d_grad(x1, x2):  
    return (2 * x1, 4 * x2)

# 梯度下降
def gd_2d(x1, x2, s1, s2, f_grad, eta=0.1):
    g1, g2 = f_grad(x1, x2)
    return (x1 - eta * g1, x2 - eta * g2, 0, 0)

自适应方法

选择恰到好处的学习率 η \eta η 很棘手。如果把它选得太小优化太慢,如果太大则会振荡,甚至可能发散。如果有可以自动确定 η \eta η 或者完全不必选择学习率的方法就会方便很多,这种方法是二阶的,除了考虑目标函数的值和梯度、还需要考虑它的曲率等,计算代价较大

牛顿法

事实上函数 f : R d → R f: \mathbb{R}^d \rightarrow \mathbb{R} f:RdR 的泰勒展开式可以进一步写成
f ( x + ϵ ) = f ( x ) + ϵ ⊤ ∇ f ( x ) + 1 2 ϵ ⊤ ∇ 2 f ( x ) ϵ + O ( ∣ ϵ ∣ 3 ) f(\mathbf{x} + \boldsymbol{\epsilon}) = f(\mathbf{x}) + \boldsymbol{\epsilon}^\top \nabla f(\mathbf{x}) + \frac{1}{2} \boldsymbol{\epsilon}^\top \nabla^2 f(\mathbf{x}) \boldsymbol{\epsilon} + \mathcal{O}(|\boldsymbol{\epsilon}|^3) f(x+ϵ)=f(x)+ϵf(x)+21ϵ2f(x)ϵ+O(ϵ3)
H = d e f ∇ 2 f ( x ) \mathbf{H} \stackrel{\mathrm{def}}{=} \nabla^2 f(\mathbf{x}) H=def2f(x) 定义为 f f f H e s s i a n Hessian Hessian d × d d \times d d×d 矩阵 H \mathbf{H} H 即多元函数的二阶偏导数。当 d d d 很小且问题简单时 H \mathbf{H} H 容易计算,但是对于深度神经网络而言 H \mathbf{H} H 一般非常大很难计算, O ( d 2 ) \mathcal{O}(d^2) O(d2) 的存储代价也会很高,且通过反向传播进行计算更是雪上加霜,忽略这些考量与高阶项,取函数对 ϵ \boldsymbol{\epsilon} ϵ 的导数得
∇ f ( x ) + H ϵ = 0  and hence  ϵ = − H − 1 ∇ f ( x ) x ← x − H − 1 ∇ f ( x ) \nabla f(\mathbf{x}) + \mathbf{H} \boldsymbol{\epsilon} = 0 \text{ and hence } \boldsymbol{\epsilon} = -\mathbf{H}^{-1} \nabla f(\mathbf{x})\ \mathbf{x} \leftarrow \mathbf{x} - \mathbf{H}^{-1} \nabla f(\mathbf{x}) f(x)+Hϵ=0 and hence ϵ=H1f(x)xxH1f(x)
解优化问题需要将 H e s s i a n Hessian Hessian 矩阵 H \mathbf{H} H 求逆。例如 f ( x ) = 1 2 x 2 f(x) = \frac{1}{2} x^2 f(x)=21x2 泰勒展开式为 f ( x + ϵ ) = 1 2 x 2 + ϵ x + 1 2 ϵ 2 f(x+\epsilon)= \frac{1}{2} x^2 + \epsilon x + \frac{1}{2} \epsilon^2 f(x+ϵ)=21x2+ϵx+21ϵ2 ∇ f ( x ) = x \nabla f(x) = x f(x)=x H = 1 \mathbf{H} = 1 H=1。对于任何 x x x 可得 ϵ = − x \epsilon = -x ϵ=x f f f 的最小值满足 ∇ f = 0 \nabla f = 0 f=0,此时优化一步就完美地收敛了。对于凸双曲余弦函数 c o s h ( c x ) cosh(cx) cosh(cx) 其中 c c c 为常数,经过几次迭代后,得到了 x = 0 x=0 x=0 处的全局最小值

# 目标函数
def f(x,c=0.5):  
    return torch.cosh(c * x)

# 目标函数的梯度	
def f_grad(x):  
    return c * torch.sinh(c * x)

# 目标函数的Hessian
def f_hessian(x):  
    return c**2 * torch.cosh(c * x)

# 牛顿法
def nt(eta=1):
    x = 10.0
    results = [x]
    for i in range(10):
        x -= eta * f_grad(x) / f_hess(x)
        results.append(float(x))
    return results


考虑非凸函数,如 f ( x ) = x cos ⁡ ( c x ) f(x) = x \cos(c x) f(x)=xcos(cx) 其中 c c c 为常数。如果 H e s s i a n Hessian Hessian 矩阵非正定即二阶导数是负的 f f f 的值可能会趋于增加

# 目标函数
def f(x, c=0.15*np.pi):  
    return x * torch.cos(c * x)
  
# 目标函数的梯度
def f_grad(x):  
    return torch.cos(c * x) - c * x * torch.sin(c * x)

# 目标函数的Hessian
def f_hessian(x):  
    return - 2 * c * torch.sin(c * x) - x * c**2 * torch.cos(c * x)


可以通过取 H e s s i a n Hessian Hessian 的绝对值修正或重新引入学习率。虽然违背了牛顿法的初衷,但因为有二阶信息使曲率较大时保持优化缓慢,而在目标函数较平坦时则优化迅速

收敛性分析

以凸函数 f f f 为例,分析牛顿法收敛速度。目标凸函数三次可微,而且二阶导数不为零,即 f ′ ′ > 0 f'' > 0 f>0

x ( k ) x^{(k)} x(k) 表示 x x x 在第 k t h k^\mathrm{th} kth 次迭代时的值,令 e ( k ) = d e f x ( k ) − x ∗ e^{(k)} \stackrel{\mathrm{def}}{=} x^{(k)} - x^* e(k)=defx(k)x 表示 k t h k^\mathrm{th} kth 迭代时与最优值的距离。通过泰勒展开得条件 f ′ ( x ∗ ) = 0 f'(x^*) = 0 f(x)=0 可以写成
0 = f ′ ( x ( k ) − e ( k ) ) = f ′ ( x ( k ) ) − e ( k ) f ′ ′ ( x ( k ) ) + 1 2 ( e ( k ) ) 2 f ′ ′ ′ ( ξ ( k ) ) 0 = f'(x^{(k)} - e^{(k)}) = f'(x^{(k)}) - e^{(k)} f''(x^{(k)}) + \frac{1}{2} (e^{(k)})^2 f'''(\xi^{(k)}) 0=f(x(k)e(k))=f(x(k))e(k)f(x(k))+21(e(k))2f(ξ(k))
对于 ξ ( k ) ∈ [ x ( k ) − e ( k ) , x ( k ) ] \xi^{(k)} \in [x^{(k)} - e^{(k)}, x^{(k)}] ξ(k)[x(k)e(k),x(k)] 成立。将上述展开除以 f ′ ′ ( x ( k ) ) f''(x^{(k)}) f(x(k)) 得到
e ( k ) − f ′ ( x ( k ) ) f ′ ′ ( x ( k ) ) = 1 2 ( e ( k ) ) 2 f ′ ′ ′ ( ξ ( k ) ) f ′ ′ ( x ( k ) ) e^{(k)} - \frac{f'(x^{(k)})}{f''(x^{(k)})} = \frac{1}{2} (e^{(k)})^2 \frac{f'''(\xi^{(k)})}{f''(x^{(k)})} e(k)f(x(k))f(x(k))=21(e(k))2f(x(k))f(ξ(k))
在牛顿法中 x ( k + 1 ) = x ( k ) − f ′ ( x ( k ) ) / f ′ ′ ( x ( k ) ) x^{(k+1)} = x^{(k)} - f'(x^{(k)}) / f''(x^{(k)}) x(k+1)=x(k)f(x(k))/f(x(k))。代入更新方程,取两边的绝对值得到
∣ e ( k + 1 ) ∣ = 1 2 ( e ( k ) ) 2 ∣ f ′ ′ ′ ( ξ ( k ) ) ∣ f ′ ′ ( x ( k ) ) \left|e^{(k+1)}\right| = \frac{1}{2}(e^{(k)})^2 \frac{\left|f'''(\xi^{(k)})\right|}{f''(x^{(k)})} e(k+1)=21(e(k))2f(x(k))f(ξ(k))
因此,每当处于有界区域 ∣ f ′ ′ ′ ( ξ ( k ) ) ∣ / ( 2 f ′ ′ ( x ( k ) ) ) ≤ c \left|f'''(\xi^{(k)})\right| / (2f''(x^{(k)})) \leq c f(ξ(k))/(2f(x(k)))c, 就有一个二次递减误差
∣ e ( k + 1 ) ∣ ≤ c ( e ( k ) ) 2 \left|e^{(k+1)}\right| \leq c (e^{(k)})^2 e(k+1)c(e(k))2
称之为线性收敛,而将 ∣ e ( k + 1 ) ∣ ≤ α ∣ e ( k ) ∣ \left|e^{(k+1)}\right| \leq \alpha \left|e^{(k)}\right| e(k+1)αe(k) 这样的条件称为恒定收敛速度。虽然无法估计整体收敛的速度,但是一旦接近极小值,收敛将变得非常快。另外,这种分析要求 f f f 在高阶导数上表现良好,即确保 f f f 值的变化过程没有超常的特性

梯度预处理

计算和存储完整的 H e s s i a n Hessian Hessian 非常昂贵,可以避免计算整个 H e s s i a n Hessian Hessian,而只计算对角线项用来进行梯度预处理
x ← x − η d i a g ( H ) − 1 ∇ f ( x ) \mathbf{x} \leftarrow \mathbf{x} - \eta \mathrm{diag}(\mathbf{H})^{-1} \nabla f(\mathbf{x}) xxηdiag(H)1f(x)
数据的预处理中若一个变量以毫米表示高度,另一个变量以公里表示高度,假设这两种自然尺度都以米为单位,那么参数化后就出现了严重的不匹配,使用预处理如归一化 *** 作有利于优化。梯度下降过程中乘 H e s s i a n Hessian Hessian 对角项相当于为每个变量的梯度选择了不同的学习率

线搜索

梯度下降的一个关键问题是可能会优化过冲或不足,解决这一问题的简单方法是在梯度下降中结合线搜索。通过 ∇ f ( x ) \nabla f(\mathbf{x}) f(x) 给出的方向,然后进行二分搜索,以确定哪个学习率 η \eta η 使 f ( x − η ∇ f ( x ) ) f(\mathbf{x} - \eta \nabla f(\mathbf{x})) f(xηf(x)) 取最小值,虽然收敛迅速但是线搜索的每一步都需要评估整个数据集上的目标函数,实现方式代价太高

随机梯度下降

机器学习的数据量庞大,直接使用梯度下降计算代价较大,通常使用随机梯度下降替代

随机梯度更新

在深度学习中,目标函数通常是训练数据集中每个样本的损失函数的平均值。给定 n n n 个样本的训练数据集,我们假设 f i ( x ) f_i(\mathbf{x}) fi(x) 是关于索引 i i i 的训练样本的损失函数,其中 x \mathbf{x} x 是参数向量,可得目标函数
f ( x ) = 1 n ∑ i = 1 n f i ( x ) f(\mathbf{x}) = \frac{1}{n} \sum_{i = 1}^n f_i(\mathbf{x}) f(x)=n1i=1nfi(x)
x \mathbf{x} x 的目标函数的梯度计算为
∇ f ( x ) = 1 n ∑ i = 1 n ∇ f i ( x ) \nabla f(\mathbf{x}) = \frac{1}{n} \sum_{i = 1}^n \nabla f_i(\mathbf{x}) f(x)=n1i=1nfi(x)
如果使用梯度下降法,则每个自变量迭代的计算代价为 O ( n ) \mathcal{O}(n) O(n) n n n 线性增长。当训练数据集较大时,每次迭代的梯度下降计算代价将较高

随机梯度下降可降低每次迭代时的计算代价。在随机梯度下降的每次迭代中,我们对数据样本随机均匀采样一个索引 i i i,其中 i ∈ 1 , … , n i\in{1,\ldots, n} i1,,n,并计算梯度 ∇ f i ( x ) \nabla f_i(\mathbf{x}) fi(x) 以更新 x \mathbf{x} x
x ← x − η ∇ f i ( x ) \mathbf{x} \leftarrow \mathbf{x} - \eta \nabla f_i(\mathbf{x}) xxηfi(x)
其中 η \eta η 是学习率,每次迭代的计算代价从梯度下降的 O ( n ) \mathcal{O}(n) O(n) 降至常数 O ( 1 ) \mathcal{O}(1) O(1)。随机梯度 ∇ f i ( x ) \nabla f_i(\mathbf{x}) fi(x) 是对完整梯度 ∇ f ( x ) \nabla f(\mathbf{x}) f(x) 的无偏估计
E i ∇ f i ( x ) = 1 n ∑ i = 1 n ∇ f i ( x ) = ∇ f ( x ) \mathbb{E}_i \nabla f_i(\mathbf{x}) = \frac{1}{n} \sum_{i = 1}^n \nabla f_i(\mathbf{x}) = \nabla f(\mathbf{x}) Eifi(x)=n1i=1nfi(x)=f(x)
平均而言,随机梯度是对梯度的良好估计,随机梯度下降可认为是向梯度添加均值为 0 0 0、方差为 1 1 1 的随机噪声的梯度下降

# 目标函数
def f(x1, x2):  
    return x1 ** 2 + 2 * x2 ** 2

# 目标函数的梯度
def f_grad(x1, x2):  
    return 2 * x1, 4 * x2

def sgd(x1, x2, s1, s2, f_grad, eta=0.1, lr=1):
    g1, g2 = f_grad(x1, x2)
    # 模拟有噪声的梯度
    g1 += torch.normal(0.0, 1, (1,))
    g2 += torch.normal(0.0, 1, (1,))
    eta_t = eta * lr()
    return (x1 - eta_t * g1, x2 - eta_t * g2, 0, 0)


随机梯度下降中变量的轨迹比梯度下降中观察到的轨迹嘈杂得多。由于梯度的随机性质,即使接近最小值仍然受到通过 η ∇ f i ( x ) \eta \nabla f_i(\mathbf{x}) ηfi(x) 的瞬间梯度所注入的不确定性的影响。经过 50 50 50 次迭代,结果仍然不那么好。并且经过额外的步骤也不会得到改善

因此唯一的选择就只能改变学习率 η \eta η,但如果选择的学习率太小或太大都不合适,解决这些相互冲突的目标的唯一方法是在优化过程中动态降低学习率,用与时间相关的学习率 η ( t ) \eta(t) η(t) 取代 η \eta η 增加了控制优化算法收敛的复杂性。需要注意 η \eta η 的衰减速度,如果太快将过早停止优化。如果太慢优化将花费太多时间,随着时间推移调整 η \eta η 的一些基本策略是
η ( t ) = η i     if     t i ≤ t ≤ t i + 1 分段常数 η ( t ) = η 0 ⋅ e − λ t 指数衰减 η ( t ) = η 0 ⋅ ( β t + 1 ) − α 多项式衰减 \begin{aligned} \eta(t) & = \eta_i \;\text{ if } \;t_i \leq t \leq t_{i+1} && \text{分段常数} \ \eta(t) & = \eta_0 \cdot e^{-\lambda t} && \text{指数衰减} \ \eta(t) & = \eta_0 \cdot (\beta t + 1)^{-\alpha} && \text{多项式衰减} \end{aligned} η(t)η(t)η(t)=ηi if titti+1=η0eλt=η0(βt+1)α分段常数指数衰减多项式衰减
分段常数是训练深度网络的常见策略。或者通过指数衰减更积极地减低它,往往会导致算法收敛之前过早停止。受欢迎的选择是 α = 0.5 \alpha = 0.5 α=0.5 的多项式衰减。在凸优化的情况下,有许多证据表明表现良好

使用指数衰减

def exponential_lr():
    global t
    t += 1
    return math.exp(-0.1 * t)

t = 1
lr = exponential_lr


参数的方差大大减少。但是很难收敛到最优解 x = ( 0 , 0 ) \mathbf{x} = (0, 0) x=(0,0) 即使经过 1000 1000 1000 个迭代步骤,仍然离最优解很远。事实上,该算法根本无法收敛。使用多项式衰减,其中学习率随迭代次数的平方根倒数衰减,那么仅在 50 50 50 次迭代之后,收敛就会更好

def polynomial_lr():
    global t
    t += 1
    return (1 + 0.1 * t) ** (-0.5)

t = 1
lr = polynomial_lr


另外,还可以从较小的学习率开始,然后使其迅速上涨,再让它降低,尽管这会更慢。甚至可以在较小和较大的学习率之间切换

收敛性分析

假设所有 ξ \boldsymbol{\xi} ξ 的目标函数 f ( ξ , x ) f(\boldsymbol{\xi}, \mathbf{x}) f(ξ,x) x \mathbf{x} x 中都是凸的,考虑随机梯度下降更新
x t + 1 = x t − η t ∂ x f ( ξ t , x ) (1) \tag{1} \mathbf{x}_{t+1} = \mathbf{x}_{t} - \eta_t \partial_\mathbf{x} f(\boldsymbol{\xi}_t, \mathbf{x}) xt+1=xtηtxf(ξt,x)(1)
其中 f ( ξ t , x ) f(\boldsymbol{\xi}_t, \mathbf{x}) f(ξt,x) 是训练样本 f ( ξ t , x ) f(\boldsymbol{\xi}_t, \mathbf{x}) f(ξt,x) 的目标函数 ξ t \boldsymbol{\xi}_t ξt 从第 t t t 步的某个分布中提取, x \mathbf{x} x 是模型参数,用 R ( x ) = E ξ [ f ( ξ , x ) ] R(\mathbf{x}) = E_{\boldsymbol{\xi}}[f(\boldsymbol{\xi}, \mathbf{x})] R(x)=Eξ[f(ξ,x)] 表示期望风险, R ∗ R^* R 表示对于 x \mathbf{x} x 的最低风险。最后让 x ∗ \mathbf{x}^* x 表示最小值(我们假设它存在于定义 x \mathbf{x} x 的域中)在这种情况下,可以跟踪时间 t t t 处的当前参数 x t \mathbf{x}_t xt 和风险最小化器 x ∗ \mathbf{x}^* x 之间的距离,观察是否随着时间的推移而改善
∣ x t + 1 − x ∗ ∣ 2 = ∣ x t − η t ∂ x f ( ξ t , x ) − x ∗ ∣ 2 = ∣ x t − x ∗ ∣ 2 + η t 2 ∣ ∂ x f ( ξ t , x ) ∣ 2 − 2 η t ⟨ x t − x ∗ , ∂ x f ( ξ t , x ) ⟩ . (2) \tag{2} \begin{aligned} &|\mathbf{x}_{t+1} - \mathbf{x}^*|^2 \ =& |\mathbf{x}_{t} - \eta_t \partial_\mathbf{x} f(\boldsymbol{\xi}_t, \mathbf{x}) - \mathbf{x}^*|^2 \ =& |\mathbf{x}_{t} - \mathbf{x}^*|^2 + \eta_t^2 |\partial_\mathbf{x} f(\boldsymbol{\xi}_t, \mathbf{x})|^2 - 2 \eta_t \left\langle \mathbf{x}_t - \mathbf{x}^*, \partial_\mathbf{x} f(\boldsymbol{\xi}_t, \mathbf{x})\right\rangle. \end{aligned} ==xt+1x2xtηtxf(ξt,x)x2xtx2+ηt2xf(ξt,x)22ηtxtx,xf(ξt,x).(2)
假设随机梯度 ∂ x f ( ξ t , x ) \partial_\mathbf{x} f(\boldsymbol{\xi}_t, \mathbf{x}) xf(ξt,x) L 2 L_2 L2 范数受到常数 L L L 的限制,因此
η t 2 ∣ ∂ x f ( ξ t , x ) ∣ 2 ≤ η t 2 L 2 (3) \tag{3} \eta_t^2 |\partial_\mathbf{x} f(\boldsymbol{\xi}_t, \mathbf{x})|^2 \leq \eta_t^2 L^2 ηt2xf(ξt,x)2ηt2L2(3)
对于 x t \mathbf{x}_t xt x ∗ \mathbf{x}^* x 之间的距离如何变化的期望。事实上,对于任何具体的步骤序列,距离可能会增加,这取决于 ξ t \boldsymbol{\xi}_t ξt 因此需要点积的边界。因为对于任何凸函数 f f f,所有 x \mathbf{x} x y \mathbf{y} y 都满足 f ( y ) ≥ f ( x ) + ⟨ f ′ ( x ) , y − x ⟩ f(\mathbf{y}) \geq f(\mathbf{x}) + \langle f'(\mathbf{x}), \mathbf{y} - \mathbf{x} \rangle f(y)f(x)+f(x),yx,按凸性有
f ( ξ t , x ∗ ) ≥ f ( ξ t , x t ) + ⟨ x ∗ − x t , ∂ x f ( ξ t , x t ) ⟩ (4) \tag{4} f(\boldsymbol{\xi}_t, \mathbf{x}^*) \geq f(\boldsymbol{\xi}_t, \mathbf{x}_t) + \left\langle \mathbf{x}^* - \mathbf{x}_t, \partial_{\mathbf{x}} f(\boldsymbol{\xi}_t, \mathbf{x}_t) \right\rangle f(ξt,x)f(ξt,xt)+xxt,xf(ξt,xt)(4)
根据 ( 2 ) , ( 3 ) , ( 4 ) (2),(3),(4) (2),(3),(4) 在时间 t + 1 t+1 t+1 时获得参数之间距离的边界
∣ x t − x ∗ ∣ 2 − ∣ x ∗ t + 1 − x ∗ ∣ 2 ≥ 2 η t ( f ( ξ t , x t ) − f ( ξ t , x ∗ ) ) − η t 2 L 2 (5) \tag{5} |\mathbf{x}_{t} - \mathbf{x}^*|^2 - |\mathbf{x}*{t+1} - \mathbf{x}^*|^2 \geq 2 \eta_t (f(\boldsymbol{\xi}_t, \mathbf{x}_t) - f(\boldsymbol{\xi}_t, \mathbf{x}^*)) - \eta_t^2 L^2 xtx2xt+1x22ηt(f(ξt,xt)f(ξt,x))ηt2L2(5)
这意味着,只要当前损失和最优损失之间的差异超过 η t L 2 / 2 \eta_t L^2/2 ηtL2/2,我们就会取得进展。由于这种差异必然会收敛到零,因此学习率 η t \eta_t ηt 也需要消失,根据 ( 5 ) (5) (5) 取期望得
E [ ∣ x t − x ∗ ∣ 2 ] − E [ ∣ x t + 1 − x ∗ ∣ 2 ] ≥ 2 η t [ E [ R ( x t ) ] − R ∗ ] − η t 2 L 2 (6) \tag{6} E\left[|\mathbf{x}_{t} - \mathbf{x}^*|^2\right] - E\left[|\mathbf{x}_{t+1} - \mathbf{x}^*|^2\right] \geq 2 \eta_t [E[R(\mathbf{x}_t)] - R^*] - \eta_t^2 L^2 E[xtx2]E[xt+1x2]2ηt[E[R(xt)]R]ηt2L2(6)
最后一步是对 t ∈ 1 , … , T t \in {1, \ldots, T} t1,,T 的不等式求和。在求和过程中抵消中间项,然后舍去低阶项,可得
∣ x 1 − x ∗ ∣ 2 ≥ 2 ( ∑ t = 1 T η t ) [ E [ R ( x t ) ] − R ∗ ] − L 2 ∑ t = 1 T η t 2 (7) \tag{7} |\mathbf{x}_1 - \mathbf{x}^*|^2 \geq 2 \left (\sum_{t=1}^T \eta_t \right) [E[R(\mathbf{x}_t)] - R^*] - L^2 \sum_{t=1}^T \eta_t^2 x1x22(t=1Tηt)[E[R(xt)]R]L2t=1Tηt2(7)
利用了给定的 x 1 \mathbf{x}_1 x1,因而可以去掉期望,最后定义
x ˉ = d e f ∑ t = 1 T η t x t ∑ t = 1 T η t (8) \tag{8} \bar{\mathbf{x}} \stackrel{\mathrm{def}}{=} \frac{\sum_{t=1}^T \eta_t \mathbf{x}_t}{\sum_{t=1}^T \eta_t} xˉ=deft=1Tηtt=1Tηtxt(8)
因为有
E ( ∑ t = 1 T η t R ( x t ) ∑ t = 1 T η t ) = ∑ t = 1 T η t E [ R ( x t ) ] ∑ t = 1 T η t = E [ R ( x t ) ] (9) \tag{9} E\left(\frac{\sum_{t=1}^T \eta_t R(\mathbf{x}_t)}{\sum_{t=1}^T \eta_t}\right) = \frac{\sum_{t=1}^T \eta_t E[R(\mathbf{x}_t)]}{\sum_{t=1}^T \eta_t} = E[R(\mathbf{x}_t)] E(t=1Tηtt=1TηtR(xt))=t=1Tηtt=1TηtE[R(xt)]=E[R(xt)](9)
根据詹森不等式中( i = t i=t i=t α i = η t / ∑ t = 1 T η t \alpha_i = \eta_t/\sum_{t=1}^T \eta_t αi=ηt/t=1Tηt)和 R R R 的凸性使其满足的 E [ R ( x t ) ] ≥ E [ R ( x ˉ ) ] E[R(\mathbf{x}_t)] \geq E[R(\bar{\mathbf{x}})] E[R(xt)]E[R(xˉ)],因此
∑ t = 1 T η t E [ R ( x t ) ] ≥ ∑ t = 1 T η t E [ R ( x ˉ ) ] (10) \tag{10} \sum_{t=1}^T \eta_t E[R(\mathbf{x}_t)] \geq \sum_{t=1}^T \eta_t E\left[R(\bar{\mathbf{x}})\right] t=1TηtE[R(xt)]t=1TηtE[R(xˉ)](10)
将其代入不等式 ( 7 ) (7) (7) 得到边界
[ E [ x ˉ ] ] − R ∗ ≤ r 2 + L 2 ∑ t = 1 T η t 2 2 ∑ t = 1 T η t (11) \tag{11} \left[E[\bar{\mathbf{x}}]\right] - R^* \leq \frac{r^2 + L^2 \sum_{t=1}^T \eta_t^2}{2 \sum_{t=1}^T \eta_t} [E[xˉ]]R2t=1Tηtr2+L2t=1Tηt2(11)
其中 r 2 = d e f ∣ x 1 − x ∗ ∣ 2 r^2 \stackrel{\mathrm{def}}{=} |\mathbf{x}_1 - \mathbf{x}^*|^2 r2=defx1x2 是初始选择参数与最终结果之间距离的边界。简而言之,收敛速度取决于随机梯度标准的限制方式 L L L 以及初始参数值与最优结果的距离 r r r 。请注意,边界由 x ˉ \bar{\mathbf{x}} xˉ 而不是 x T \mathbf{x}_T xT 表示。因为 x ˉ \bar{\mathbf{x}} xˉ 是优化路径的平滑版本。只要知道 r , L r, L r,L T T T,就可以选择学习率 η = r / ( L T ) \eta = r/(L \sqrt{T}) η=r/(LT )。这个就是上界 r L / T rL/\sqrt{T} rL/T ,也就是说将按照速度 O ( 1 / T ) \mathcal{O}(1/\sqrt{T}) O(1/T ) 收敛到最优解

小批量随机梯度下降

基于梯度的学习方法中梯度下降中使用完整数据集来计算梯度并更新参数,随机梯度下降中一次处理一个训练样本更新参数。二者各有利弊,每当数据非常相似时,梯度下降并不是非常数据高效。而由于 C P U CPU CPU G P U GPU GPU 无法充分利用向量化,随机梯度下降并不特别计算高效,小批量随机梯度下降则是两者的折中方案

向量化和缓存

使用小批量的决策的核心是计算效率。当考虑与多个 G P U GPU GPU 和多台服务器并行处理时,我们需要向每个 G P U GPU GPU 发送至少一张图像。有了每台服务器 8 8 8 G P U GPU GPU 16 16 16 台服务器,就能得到大小为 128 128 128 的小批量

当只有单个 G P U GPU GPU 甚至 C P U CPU CPU 时,事情会更微妙一些,这些设备有多种类型的内存,此时需要考虑多种类型的计算单元以及在它们之间不同的带宽限制。例如,一个 C P U CPU CPU 有少量寄存器, L 1 L1 L1 L 2 L2 L2 缓存,以及 L 3 L3 L3 缓存(在不同的处理器内核之间共享)随着缓存的大小的增加,延迟也在增加,同时带宽在减少,而处理器能够执行的 *** 作远比主内存接口所能提供的多得多

首先,具有 16 16 16 个内核和 AVX-512 \text{AVX-512} AVX-512 向量化的 2 G H z    C P U 2GHz\;CPU 2GHzCPU 每秒可处理高达 2 ⋅ 1 0 9 ⋅ 16 ⋅ 32 = 1 0 12 2 \cdot 10^9 \cdot 16 \cdot 32 = 10^{12} 21091632=1012 个字节。 G P U GPU GPU 的性能一般是该数值的 100 100 100 倍。而另一方面,中端服务器处理器的带宽可能不超过 100 G b / s 100Gb/s 100Gb/s 即不到处理器满负荷所需的十分之一。且并非所有的内存入口都是相等的,内存接口通常为 64 64 64 位或更宽,因此读取单个字节会导致由于更宽的存取而产生的代价,其次,第一次存取的额外开销很大,而按序存取或突发读取相对开销较小

减轻这些限制的方法是使用足够快的 C P U CPU CPU 缓存层次结构来为处理器提供数据。例如,矩阵-矩阵乘法 A = B C \mathbf{A} = \mathbf{B}\mathbf{C} A=BC 有很多方法来计算 A \mathbf{A} A

  1. 可以通过点积进行逐元素计算 A i j = B i , : C : , j ⊤ \mathbf{A}_{ij} = \mathbf{B}_{i,:} \mathbf{C}_{:,j}^\top Aij=Bi,:C:,j,每次计算一个元素 A i j \mathbf{A}_{ij} Aij 时,都需要将一行和一列向量复制到 C P U CPU CPU 中。由于矩阵元素是按顺序对齐的,因此当从内存中读取时,需要访问两个向量中许多不相交的位置,访问速度低
  2. 可以一次计算一列计算 A : , j = B C : , j ⊤ \mathbf{A}_{:,j} = \mathbf{B} \mathbf{C}_{:,j}^\top A:,j=BC:,j,同样也可以一次计算 A \mathbf{A} A一行 A i , : \mathbf{A}_{i,:} Ai,:,相对更有利,能够在遍历 B \mathbf{B} B 的同时,将列向量 C : , j \mathbf{C}_{:,j} C:,j 保留在 C P U CPU CPU 缓存中,内存带宽需求减半,相应地提高了访问速度
  3. 可以直接计算 A = B C \mathbf{A} = \mathbf{B} \mathbf{C} A=BC,表面上是最可取的,然而大多数矩阵可能不能完全放入缓存中
  4. 可以将 B \mathbf{B} B C \mathbf{C} C 分成较小的区块矩阵,然后一次计算 A \mathbf{A} A 的一个区块,实践上有用的方案,可以将矩阵的区块移到缓存中然后在本地将它们相乘

除此之外,编程语言和学习框架本身带来的额外开销也是相当大

小批量

处理单个观测值需要执行许多单一矩阵-矢量(甚至矢量-矢量)乘法,执行 w ← w − η t g t \mathbf{w} \leftarrow \mathbf{w} - \eta_t \mathbf{g}_t wwηtgt 时消耗大,其中
g t = ∂ w f ( x t , w ) \mathbf{g}_t = \partial_{\mathbf{w}} f(\mathbf{x}_{t}, \mathbf{w}) gt=wf(xt,w)
通过将其应用于一个小批量观测值来提高此 *** 作的计算效率,将梯度 g t \mathbf{g}_t gt 替换为一个小批量而不是单个观测值
g t = ∂ w 1 ∣ B t ∣ ∑ i ∈ B t f ( x i , w ) \mathbf{g}_t = \partial_{\mathbf{w}} \frac{1}{|\mathcal{B}_t|} \sum_{i \in \mathcal{B}_t} f(\mathbf{x}_{i}, \mathbf{w}) gt=wBt1iBtf(xi,w)
这对 g t \mathbf{g}_t gt 的统计属性有什么影响,由于 x t \mathbf{x}_t xt 和小批量 B t \mathcal{B}_t Bt 的所有元素都是从训练集中随机抽出的,因此梯度的期望保持不变。另一方面,方差显著降低。由于小批量梯度由正在被平均计算的 b : = ∣ B t ∣ b := |\mathcal{B}_t| b:=Bt 个独立梯度组成,其标准差降低了 b − 1 2 b^{-\frac{1}{2}} b21。这本身是一件好事,因为这意味着更新与完整的梯度更接近了

直观来说,这表明选择大型的小批量 B t \mathcal{B}_t Bt 将是普遍可行的。然而与计算代价的线性增长相比,标准差的额外减少是微乎其微的。在实践中选择一个足够大的小批量匹配 C P U CPU CPU G P U GPU GPU 的内存即可

小结
  • 学习率的大小很重要,学习率太大会使模型发散,学习率太小会优化不足

  • 梯度下降会可能陷入局部极小值,而得不到全局最小值

  • 在高维模型中,调整学习率是很复杂的,可以使用 H e s s i a n Hessian Hessian 矩阵的对角进行预处理,调节学习率比例

  • 牛顿法在凸问题中一旦开始正常工作,速度就会快得多,但对于非凸问题,不要不作任何调整就使用牛顿法

  • 对于凸问题,可以证明选择合适范围内的学习率,随机梯度下降将收敛到最优解

  • 如果学习率太小或太大都会出现问题,通常只有经过多次实验后才能找到合适的学习率

  • 当训练数据集中有大量样本时,随机梯度下降的每次迭代的代价更低

  • 随机梯度下降在非凸情况下一般不可用,但是其梯度中一定程度的噪声可能会使参数跳出局部最小值

  • 由于减少了学习框架的额外开销,使用更好的内存定位以及 C P U CPU CPU G P U GPU GPU 上的缓存,向量化使代码更加高效

  • 随机梯度下降的统计效率与梯度下降一次处理数据的计算效率之间存在权衡,小批量随机梯度下降提供了均衡的方案

  • 一般来说,在训练期间降低学习率有助于训练,小批量随机梯度下降比随机梯度下降和梯度下降的优化速度快,收敛风险小

启发式算法 退火算法 粒子群算法 遗传算法 蚁群算法

欢迎分享,转载请注明来源:内存溢出

原文地址:https://54852.com/langs/921170.html

(0)
打赏 微信扫一扫微信扫一扫 支付宝扫一扫支付宝扫一扫
上一篇 2022-05-16
下一篇2022-05-16

发表评论

登录后才能评论

评论列表(0条)

    保存