最优化学习 数值优化的例子：实现最小二乘法

全部笔记的汇总贴：最优化学习目录

---

上溢和下溢

这里先介绍一个数值优化常出现的一个概念上溢和下溢
下溢（Underflow）：当接近零的数被四舍五⼊为零时发生下溢。
上溢（Overflow）：当⼤量级的数被近似为$\infty$或$-\infty$时发⽣上溢。

必须对上溢和下溢进行数值稳定的⼀个例子是softmax 函数。softmax 函数经常用于预测与范畴分布相关联的概率，定义为:
$$\mathrm{softmax}(\mathbf{x}{)}_i=\frac{\exp \left(x_i\right)}{{\sum }_{j=1}^n\exp \left(x_j\right)}$$

最小二乘法介绍

我们接下来会分别用梯度下降法、牛顿法、约束优化法分别解我们的最小二乘法，那我们首先要明白最小二乘法。

最小二乘法（Least Squares）是回归分析中的一种标准方法，它是用来近似超定系统（Overdetermined System）答案的一种方法。超定系统是指数学中的一种概念，一组包含未知数的方程组中，如果方程的数量大于未知数的数量，那么这个系统就是一个超定系统（超定方程组）。超定系统（超定方程组）一般是无解的，只能求近似解。而最小二乘法就是求超定方程组近似解的一种方法。
举个通俗的例子，如下二维平面图中有很多个点，假设我们想用一条直线来拟合数据，即期望能找到一条直线能最好地穿过这些数据点。

我们这里引入线性最小二乘法：

$$\min f(\mathbf{x})=\frac{1}{2}\Vert \mathbf{A}\mathbf{x}-\mathbf{b}{\Vert }_2^2$$

梯度下降法 Gradient Decent

在我的最速下降法（steepest Descent)中介绍了梯度下降法和最速下降法，梯度下降法就是一阶可微的最速下降法，在这里我们可以对其进行求解。
假设我们希望找到最小化该式的$x$值，那我们需要得到它的梯度：
$${\nabla }_{x}f(x)=A^{\top }(Ax-b)=A^{\top }Ax-A^{\top }b$$
梯度下降法建议我们新的点更新为
$$x^{\prime }=x-ϵ{\nabla }_{x}f(x)$$
其中ϵ 为学习率(learning rate)，是⼀个确定步长大小的正标量。

首先我们给定初始数据

# 导入相应的库
import numpy as np
import numpy.linalg as la

x0 = np.array([1.0, 1.0, 1.0])
A = np.array([[1.0, -2.0, 1.0], [0.0, 2.0, -8.0], [-4.0, 5.0, 9.0]])
b = np.array([0.0, 8.0, -9.0])
epsilon = 0.001
delta = 1e-3
# 给定A，b，真正的解x 为[29, 16, 3]

梯度下降法

"""
梯度下降法
"""
def matmul_chain(*args): # 定义一个矩阵连乘的函数
    if len(args) == 0:
        return np.nan
    result = args[0]
    for x in args[1:]:
        result = result @ x
    return result

def gradient_decent(x, A, b, epsilon, delta):
    while la.norm(matmul_chain(A.T, A, x) - matmul_chain(A.T, b)) > delta: # 默认ord = 2,没有到达边界
        x -= epsilon*(matmul_chain(A.T, A, x) - matmul_chain(A.T, b)) # 更新x
    return x

gradient_decent(x0, A, b, epsilon, delta)

array([27.82277014, 15.34731055, 2.83848939])

牛顿法（Newton’s Method）

具体的推导过程和结果可以看牛顿法（Newton’s method）
我们的牛顿法基于我们的优化函数能够二阶可微，这样我们就可以利用一个二阶泰勒展开来近似$x^0$的$f(x)$:
$$f(x)\approx f\left(x^{(0)}\right)+{\left(x-x^{(0)}\right)}^{\top }{\nabla }_{x}f\left(x^{(0)}\right)+\frac{1}{2}{\left(x-x^{(0)}\right)}^{\top }\mathrm{H}\left(x^{(0)}\right)\left(x-x^{(0)}\right)$$
然后根据对$f(x)$求导可以得到我们临界点：

$${\mathbf{x}}^{\ast }={\mathbf{x}}^{(0)}-\mathrm{H}{\left({\mathbf{x}}^{(0)}\right)}^{-1}{\nabla }_{\mathbf{x}}f\left({\mathbf{x}}^{(0)}\right)$$

牛顿法迭代地更新近似函数和跳到近似函数的最小点可以比梯度下降法更快地到达临界点。这在接近全局极小时是⼀个特别有用的性质，但是在鞍点附近是有害的。

对于我们的最小二乘法，我们可以计算得到，我们的海瑟矩阵其实就是$\mathrm{H}=A^{\top }A$
进一步计算我们的最优解
$$x^{\ast }=x^{(0)}-{\left(A^{\top }A\right)}^{-1}\left(A^{\top }A{x}^{(0)}-A^{\top }b\right)={\left(A^{\top }A\right)}^{-1}{A}^{\top }b$$

"""
牛顿法
"""

def newton(x, A, b, delta):
    x = matmul_chain(np.linalg.inv(matmul_chain(A.T, A)), A.T, b)
    return x

newton(x0, A, b, delta)

array([29., 16., 3.])

我们可以看出来，在我们的牛顿法（Newton’s method）中提过，对于二阶的函数，我们几乎可以达到一步就到的情况，这个就是其中的例子

约束优化

这一部分内容详细推导可以看我的KKT条件(最优解的一阶必要条件) 和约束优化问题
我们希望通过$m$个函数$g^{(i)}$和$n$个函数$h^{(j)}$ 描述$S$,那么$S$可以表示为$\mathrm{S}=\{\mathbf{x}\mid \forall i,g^{(i)}(\mathbf{x})=0$ and $\forall j,h^{(j)}(\mathbf{x})\le 0\}$。其中涉及$g^{(i)}$的等式称为等式约束，涉及$h^{(j)}$的不等式称为不等式约束。
我们为每个约束引⼊新的变量${\lambda }_i$和${\alpha }_j$，这些新变量被称为KKT 乘子。⼴义拉格朗日式可以如下定义：
$$L(\mathbf{x},\lambda ,\alpha )=f(\mathbf{x})+\sum _i{\lambda }_i{g}^{(i)}(\mathbf{x})+\sum _j{\alpha }_j{h}^{(j)}(\mathbf{x})$$
可以通过优化无约束的广义拉格朗日式解决约束最小化问题：
$$\underset{x}{\min }\underset{\lambda }{\max }\underset{\alpha ,\alpha \ge 0}{\max }L(x,\lambda ,\alpha )$$
优化该式与下式等价：
$$\underset{x\in \mathbb{S}}{\min }f(\mathbf{x})$$
针对上述实例，约束优化:${\mathbf{x}}^{\top }\mathbf{x}\le 1$
引⼊广义拉格朗日式：
$$L(\mathbf{x},\lambda )=f(\mathbf{x})+\lambda \left({\mathbf{x}}^{\top }\mathbf{x}-1\right)$$
解决以下问题：
$$\underset{\mathbf{x}}{\min }\underset{\lambda ,\lambda \ge 0}{\max }L(\mathbf{x},\lambda )$$
关于$x$对$Lagrangian$微分，我们得到方程:
$${\mathbf{A}}^{\top }\mathbf{A}\mathbf{x}-{\mathbf{A}}^{\top }\mathbf{b}+2\lambda \mathbf{x}=0$$
得到解的形式是：$$\mathbf{x}={\left({\mathbf{A}}^{\top }\mathbf{A}+2\lambda \mathbf{I}\right)}^{-1}{\mathbf{A}}^{\top }\mathbf{b}$$

$\lambda$的选择必须使结果服从约束，可以对$\lambda$梯度上升找到这个值：
$$\frac{\partial }{\partial \lambda }L(x,\lambda )={\mathbf{x}}^{\top }\mathbf{x}-1$$

"""
约束优化，约束解的大小
"""
def constrain_opti(x, A, b, delta):
    k = len(x)
    lamb = 0
    while np.abs(np.dot(x.T, x) - 1) > 1e-5: # delta设为5e-2, 最优设为0
        x = matmul_chain(np.linalg.inv(matmul_chain(A.T, A)+2*lamb*np.identity(k)), A.T, b)
        lamb += np.dot(x.T, x)-1
    return x

constrain_opti(x0, A, b, delta)

array([ 0.23637902, 0.05135858, -0.94463626])