最优化学习 数值优化的例子：实现最小二乘法
Posted Real&Love
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了最优化学习 数值优化的例子：实现最小二乘法相关的知识，希望对你有一定的参考价值。
数值优化的例子：实现最小二乘法
全部笔记的汇总贴：最优化学习目录
上溢和下溢
这里先介绍一个数值优化常出现的一个概念上溢和下溢
下溢（Underflow）：当接近零的数被四舍五⼊为零时发生下溢。
上溢（Overflow）：当⼤量级的数被近似为
∞
\\infin
∞或
−
∞
-\\infin
−∞时发⽣上溢。
必须对上溢和下溢进行数值稳定的⼀个例子是softmax 函数。softmax 函数经常用于预测与范畴分布相关联的概率，定义为:
softmax
(
x
)
i
=
exp
(
x
i
)
∑
j
=
1
n
exp
(
x
j
)
\\operatorname{softmax}(\\boldsymbol{x})_{i}=\\frac{\\exp \\left(x_{i}\\right)}{\\sum_{j=1}^{n} \\exp \\left(x_{j}\\right)}
softmax(x)i=∑j=1nexp(xj)exp(xi)
最小二乘法介绍
我们接下来会分别用梯度下降法、牛顿法、约束优化法分别解我们的最小二乘法，那我们首先要明白最小二乘法。
最小二乘法（Least Squares）是回归分析中的一种标准方法，它是用来近似超定系统（Overdetermined System）答案的一种方法。超定系统是指数学中的一种概念，一组包含未知数的方程组中，如果方程的数量大于未知数的数量，那么这个系统就是一个超定系统（超定方程组）。超定系统（超定方程组）一般是无解的，只能求近似解。而最小二乘法就是求超定方程组近似解的一种方法。
举个通俗的例子，如下二维平面图中有很多个点，假设我们想用一条直线来拟合数据，即期望能找到一条直线能最好地穿过这些数据点。
我们这里引入线性最小二乘法：
min f ( x ) = 1 2 ∥ A x − b ∥ 2 2 \\min f(\\boldsymbol{x})=\\frac{1}{2}\\|\\boldsymbol{A} \\boldsymbol{x}-\\boldsymbol{b}\\|_{2}^{2} minf(x)=21∥Ax−b∥22
梯度下降法 Gradient Decent
在我的最速下降法（steepest Descent)中介绍了梯度下降法和最速下降法，梯度下降法就是一阶可微的最速下降法，在这里我们可以对其进行求解。
假设我们希望找到最小化该式的
x
x
x值，那我们需要得到它的梯度：
∇
x
f
(
x
)
=
A
⊤
(
A
x
−
b
)
=
A
⊤
A
x
−
A
⊤
b
\\nabla_{x} f(x)=A^{\\top}(A x-b)=A^{\\top} A x-A^{\\top} b
∇xf(x)=A⊤(Ax−b)=A⊤Ax−A⊤b
梯度下降法建议我们新的点更新为
x
′
=
x
−
ϵ
∇
x
f
(
x
)
x^{\\prime}=x-\\epsilon \\nabla_{x} f(x)
x′=x−ϵ∇xf(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
x^{0}
x0的
f
(
x
)
f(x)
f(x):
f
(
x
)
≈
f
(
x
(
0
)
)
+
(
x
−
x
(
0
)
)
⊤
∇
x
f
(
x
(
0
)
)
+
1
2
(
x
−
x
(
0
)
)
⊤
H
(
x
(
0
)
)
(
x
−
x
(
0
)
)
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)≈f(x(0))+(x−x(0))⊤∇xf(x(0))+21(x−x(0))⊤H(x(0))(x−x(0))
然后根据对
f
(
x
)
f(x)
f(x)求导可以得到我们临界点：
x ∗ = x ( 0 ) − H ( x ( 0 ) ) − 1 ∇ x f ( x ( 0 ) ) \\boldsymbol{x}^{*}=\\boldsymbol{x}^{(0)}-\\mathrm{H}\\left(\\boldsymbol{x}^{(0)}\\right)^{-1} \\nabla_{\\boldsymbol{x}} f\\left(\\boldsymbol{x}^{(0)}\\right) x∗=x(0)−H(x(0))−1∇xf(x(0))
牛顿法迭代地更新近似函数和跳到近似函数的最小点可以比梯度下降法更快地到达临界点。这在接近全局极小时是⼀个特别有用的性质，但是在鞍点附近是有害的。
对于我们的最小二乘法，我们可以计算得到，我们的海瑟矩阵其实就是
H
=
A
⊤
A
\\mathrm{H}=A^{\\top} A
H=A⊤A
进一步计算我们的最优解
x
∗
=
x
(
0
)
−
(
A
⊤
A
)
−
1
(
A
⊤
A
x
(
0
)
−
A
⊤
b
)
=
(
A
⊤
A
)
−
1
A
⊤
b
x^{*}=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
x∗=x(0)−(A⊤A)−1(A⊤Ax(0)−A⊤b)=(A⊤A)−1A⊤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
m
m个函数
g
(
i
)
g^{(i)}
g(i)和
n
n
n个函数
h
(
j
)
h^{(j)}
h(j) 描述
S
S
S,那么
S
S
S可以表示为
S
=
{
x
∣
∀
i
,
g
(
i
)
(
x
)