预备知识
Hesse 矩阵
函数
f
(
x
)
f(x)
f(x)为自变量为为向量的实值函数,其中
x
=
[
x
1
,
x
2
,
.
.
.
,
x
n
]
x = [x_1, x_2,...,x_n]
x=[x1,x2,...,xn],则Hesse矩阵的定义为:
H
(
f
)
=
[
∂
2
f
∂
x
1
2
∂
2
f
∂
x
1
∂
x
2
⋯
∂
2
f
∂
x
1
∂
x
n
∂
2
f
∂
x
2
∂
x
1
∂
2
f
∂
x
2
2
⋯
∂
2
f
∂
x
2
∂
x
n
⋮
⋮
⋱
⋮
∂
2
f
∂
x
n
∂
x
1
∂
2
f
∂
x
n
∂
x
2
⋯
∂
2
f
∂
x
n
2
]
\Large H(f)=\left[\begin{array}{cccc} \frac{\partial^{2} f}{\partial x_{1}^{2}} & \frac{\partial^{2} f}{\partial x_{1} \partial x_{2}} & \cdots & \frac{\partial^{2} f}{\partial x_{1} \partial x_{n}} \\ \frac{\partial^{2} f}{\partial x_{2} \partial x_{1}} & \frac{\partial^{2} f}{\partial x_{2}^{2}} & \cdots & \frac{\partial^{2} f}{\partial x_{2} \partial x_{n}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^{2} f}{\partial x_{n} \partial x_{1}} & \frac{\partial^{2} f}{\partial x_{n} \partial x_{2}} & \cdots & \frac{\partial^{2} f}{\partial x_{n}^{2}} \end{array}\right]
H(f)=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡∂x12∂2f∂x2∂x1∂2f⋮∂xn∂x1∂2f∂x1∂x2∂2f∂x22∂2f⋮∂xn∂x2∂2f⋯⋯⋱⋯∂x1∂xn∂2f∂x2∂xn∂2f⋮∂xn2∂2f⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤
问题
求解线性方程( 系数矩阵A对称且正定):
A
x
=
b
\Large Ax = b
Ax=b
显然,
x
=
A
−
1
b
x = A^{-1}b
x=A−1b。但是求矩阵的逆计算量太大,所以实际中使用迭代的方式求
x
x
x。
首先构造一个二次函数:
ϕ
(
x
)
=
1
2
x
T
A
x
−
b
T
x
(1)
\Large \phi(x) = \frac{1}{2}x^TAx-b^Tx \tag{1}
ϕ(x)=21xTAx−bTx(1)
对(1)求导,并令导数为0得:
∇
ϕ
(
x
)
=
A
x
−
b
=
0
(2)
\Large \nabla \phi(x)=A x-b=0 \tag{2}
∇ϕ(x)=Ax−b=0(2)
从(2)可以看出
ϕ
(
x
)
\phi(x)
ϕ(x)的导数为0,就是线性方程
A
x
=
b
Ax=b
Ax=b的解,即求
ϕ
(
x
)
\phi(x)
ϕ(x)的极小值点。现在问题转为求二次函数的极值问题。
二次函数求极值的方法(迭代方法)
迭代的方法都使用以下的原则:
最重要的就是要找到更新方向和步长。
1. 最速下降法(梯度下降法)(一阶优化方法)
梯度下降法选择函数值下降最快的方向,即负梯度方向。则有:
x
(
t
+
1
)
=
x
(
t
)
−
λ
(
t
)
∇
ϕ
(
x
)
\Large x^{(t+1)} = x^{(t)} - \lambda^{(t)}\nabla \phi(x)
x(t+1)=x(t)−λ(t)∇ϕ(x)
其中
λ
(
t
)
\lambda^{(t)}
λ(t)为步长,可以通过求
λ
(
t
)
=
a
r
g
m
i
n
f
(
x
(
t
)
−
λ
(
t
)
∇
ϕ
(
x
)
)
\lambda^{(t)} = argmin f(x^{(t)}-\lambda^{(t)}\nabla \phi(x) )
λ(t)=argminf(x(t)−λ(t)∇ϕ(x)),但由于计算代价太大,实际中设置步长为常数。
梯度下降法的收敛速度太慢,甚至可能比直接求矩阵的逆还慢。
2. 牛顿法(二阶优化方法)
给一个初始值
x
(
t
)
x^{(t)}
x(t),则
ϕ
(
x
)
\phi(x)
ϕ(x)在
x
(
t
)
x^{(t)}
x(t)处的二阶泰勒展开为:
f
(
x
)
=
ϕ
(
x
(
t
)
)
+
∇
ϕ
(
x
(
t
)
)
(
x
−
x
(
t
)
)
+
1
2
(
x
−
x
(
t
)
)
T
∇
2
ϕ
(
x
(
t
)
)
(
x
−
x
(
t
)
)
\Large f(x) =\phi(x^{(t)}) + \nabla \phi(x^{(t)})(x-x^{(t)}) + \frac{1}{2}(x-x^{(t)})^T\nabla^2 \phi(x^{(t)})(x-x^{(t)})
f(x)=ϕ(x(t))+∇ϕ(x(t))(x−x(t))+21(x−x(t))T∇2ϕ(x(t))(x−x(t))
显然
f
(
x
(
t
)
)
=
ϕ
(
x
(
t
)
)
f(x^{(t)}) = \phi(x^{(t)})
f(x(t))=ϕ(x(t)),
f
(
x
)
f(x)
f(x)可以作为
ϕ
(
x
)
\phi(x)
ϕ(x)在
x
(
t
)
x^{(t)}
x(t)的近似,则问题就变为求
f
(
x
)
f(x)
f(x)的极值, 对
f
(
x
)
f(x)
f(x)求导等于0得:
d
f
(
x
)
d
x
=
∇
ϕ
(
x
(
t
)
)
+
∇
2
ϕ
(
x
(
t
)
)
(
x
−
x
(
t
)
)
=
0
x
=
x
(
t
)
−
∇
2
ϕ
(
x
(
t
)
)
−
1
∇
ϕ
(
x
(
t
)
)
\Large \frac{df(x)}{dx} = \nabla \phi(x^{(t)}) + \nabla^2 \phi(x^{(t)})(x-x^{(t)}) = 0 \\ \Large x = x^{(t)} - \nabla^2 \phi(x^{(t)})^{-1} \nabla \phi(x^{(t)})
dxdf(x)=∇ϕ(x(t))+∇2ϕ(x(t))(x−x(t))=0x=x(t)−∇2ϕ(x(t))−1∇ϕ(x(t))
把求得得极值点
x
x
x作为
x
(
t
+
1
)
x^{(t+1)}
x(t+1),然后迭代求解。因此牛顿法,在迭代过程中需要求一阶导和二阶导,且
ϕ
(
x
)
\phi(x)
ϕ(x)的Hesse矩阵可逆。更新公式如下:
x
(
t
+
1
)
=
x
(
t
)
−
A
−
1
∇
ϕ
(
x
(
t
)
\Large x^{(t+1)} = x^{(t)} - A^{-1}\nabla \phi(x^{(t)}
x(t+1)=x(t)−A−1∇ϕ(x(t)
但是牛顿法仍需要计算矩阵的逆,因此不适用。
3. 共轭梯度法
如果要使用共轭梯度法求解线性方程,必须要求系数矩阵对称且正定。
共轭的定义:
对于一组向量
{
p
(
0
)
,
p
(
1
)
,
.
.
.
,
p
(
n
−
1
)
}
\{p^{(0)},p^{(1)},...,p^{(n-1)}\}
{p(0),p(1),...,p(n−1)},如果任意两个向量间(
i
≠
j
)
i \ne j)
i=j)满足:
(
p
(
i
)
)
T
A
p
(
j
)
=
0
(3)
\Large (p^{(i)})^T A p^{(j)} = 0 \tag{3}
(p(i))TAp(j)=0(3)
则称这组向量与对阵正定矩阵A共轭。
共轭梯度法是介于梯度下降法和牛顿法之间的一种方法。共轭梯度法初始选择负梯度方向进行更新,在后面的迭代中取负梯度方向和前一搜索方向的线性组合作为搜索方向。
对于优化问题(1)的二维情况,如下图:
其中
x
0
x
1
→
=
−
∇
ϕ
(
x
0
)
\overrightarrow{x_{0}x_{1}}=-\nabla \phi\left(x_{0}\right)
x0x1=−∇ϕ(x0)、
x
1
x
∗
→
=
−
A
−
1
∇
ϕ
(
x
1
)
\overrightarrow{x_{1} x^{*}}=-A^{-1} \nabla \phi\left(x_{1}\right)
x1x∗=−A−1∇ϕ(x1)。于是有:
x
0
x
1
→
T
A
x
1
x
∗
→
=
∇
ϕ
(
x
0
)
A
A
−
1
∇
ϕ
(
x
1
)
=
∇
ϕ
(
x
0
)
∇
ϕ
(
x
1
)
=
0
\Large \begin{aligned} \overrightarrow{x_{0} x_{1}}^{T} A \overrightarrow{x_{1} x^{*}} &=\nabla \phi\left(x_{0}\right) A A^{-1} \nabla \phi\left(x_{1}\right) \\ &=\nabla \phi\left(x_{0}\right) \nabla \phi\left(x_{1}\right) \\ &=0 \end{aligned}
x0x1TAx1x∗=∇ϕ(x0)AA−1∇ϕ(x1)=∇ϕ(x0)∇ϕ(x1)=0
这表明,两次迭代的方向是一组共轭向量。从二维推广到N维,只要找到一组共轭向量
{
p
(
0
)
,
p
(
1
)
,
.
.
.
,
p
(
n
−
1
)
}
\{p^{(0)},p^{(1)},...,p^{(n-1)}\}
{p(0),p(1),...,p(n−1)},然后依次沿着每个向量方向优化,最终在N次迭代以后就可以达到极小值。现在的问题是,如何找到一组共轭向量?
假设起始点为 x ( 0 ) x^{(0)} x(0)。 首先,选 p ( 0 ) = − ∇ ϕ ( x ( 0 ) ) p^{(0)}=- \nabla \phi(x^{(0)}) p(0)=−∇ϕ(x(0)),然后求后续的向量 p ( 1 ) , p ( 2 ) . . . p ( t ) , . . . , p ( n − 1 ) p^{(1)}, p^{(2)}...p^{(t)},...,p^{(n-1)} p(1),p(2)...p(t),...,p(n−1)。
当t=1时, 求
p
(
1
)
p^{(1)}
p(1),已知
p
(
0
)
p^{(0)}
p(0)和
∇
ϕ
(
x
(
1
)
)
\nabla\phi(x^{(1)})
∇ϕ(x(1))。因为这两个向量一定是线性无关的,所以可以在这两个向量构成的平面上寻找
p
(
1
)
p^{(1)}
p(1), 则
p
(
1
)
=
−
∇
ϕ
(
x
(
1
)
)
+
β
1
p
(
0
)
p^{(1)} = -\nabla\phi(x^{(1)}) + \beta_1 p^{(0)}
p(1)=−∇ϕ(x(1))+β1p(0)。将
p
(
0
)
,
p
(
1
)
p^{(0)}, p^{(1)}
p(0),p(1)带入(3)得:
β
1
=
(
p
(
0
)
)
T
A
∇
ϕ
(
x
(
1
)
)
(
p
(
0
)
)
T
A
p
(
0
)
\Large \beta_{1}=\frac{(p^{(0)})^{T} A \nabla \phi\left(x^{(1)}\right)}{(p^{(0)})^{T} A p^{(0)}}
β1=(p(0))TAp(0)(p(0))TA∇ϕ(x(1))
有了
β
1
\beta_1
β1,就可以得到
p
(
1
)
p^{(1)}
p(1)。依次类推,则有:
p
(
t
)
=
−
∇
ϕ
(
x
(
t
)
)
+
β
t
p
(
t
−
1
)
\Large p^{(t)} = -\nabla\phi(x^{(t)})+\beta_t p^{(t-1)}
p(t)=−∇ϕ(x(t))+βtp(t−1)
其中:
β
t
=
(
p
(
t
−
1
)
)
T
A
∇
ϕ
(
x
(
t
)
)
(
p
(
t
−
1
)
)
T
A
p
(
t
−
1
)
\Large \beta_{t}=\frac{(p^{(t-1)})^{T} A \nabla \phi\left(x^{(t)}\right)}{(p^{(t-1)})^{T} A p^{(t-1)}}
βt=(p(t−1))TAp(t−1)(p(t−1))TA∇ϕ(x(t))
现在知道了每一次更新的方向
p
(
t
)
p^{(t)}
p(t),就可以计算出
ϕ
(
x
(
t
)
)
\phi(x^{(t)})
ϕ(x(t))在该方向的最优步长
α
t
\alpha_t
αt。因为更新
x
(
t
)
x^{(t)}
x(t)的公式为:
x
(
t
+
1
)
=
x
(
t
)
+
α
t
p
(
t
)
\Large x^{(t+1)} = x^{(t)} + \alpha_t p^{(t)}
x(t+1)=x(t)+αtp(t)
则最优的
α
t
\alpha_t
αt,可以通过
ϕ
(
x
(
t
+
1
)
)
\phi(x^{(t+1)})
ϕ(x(t+1))对
α
t
\alpha_t
αt求导等于0计算,则有:
α
t
=
−
∇
ϕ
(
x
(
t
)
)
T
p
(
t
)
(
p
(
t
)
)
T
A
p
(
t
)
\Large \alpha_{t}=-\frac{\nabla \phi\left(x^{(t)}\right)^{T} p^{(t)}}{(p^{(t)})^{T} A p^{(t)}}
αt=−(p(t))TAp(t)∇ϕ(x(t))Tp(t)
综上,共轭梯度法的更新流程如下:
α
t
←
−
(
r
(
t
)
)
T
p
(
t
)
(
p
(
t
)
)
T
A
p
(
t
)
x
(
t
+
1
)
←
x
(
t
)
+
α
t
p
(
t
)
r
(
t
+
1
)
←
A
x
(
t
+
1
)
−
b
β
t
+
1
←
(
p
(
t
)
)
T
A
r
(
t
+
1
)
(
p
(
t
)
)
T
A
p
(
t
)
p
(
t
+
1
)
←
−
r
(
t
+
1
)
+
β
t
+
1
p
(
t
)
k
←
k
+
1
\Large \alpha_{t} \leftarrow -\frac{(r^{(t)})^T p^{(t)}}{(p^{(t)})^{T} A p^{(t)}} \\ \Large x^{(t+1)} \leftarrow x^{(t)} + \alpha_t p^{(t)} \\ \Large r^{(t+1)} \leftarrow Ax^{(t+1)}-b \\ \Large \beta_{t+1} \leftarrow \frac{(p^{(t)})^{T} A r^{(t+1)}}{(p^{(t)})^{T} A p^{(t)}} \\ \Large p^{(t+1)} \leftarrow -r^{(t+1)}+\beta_{t+1} p^{(t)} \\ \Large k \leftarrow k+1
αt←−(p(t))TAp(t)(r(t))Tp(t)x(t+1)←x(t)+αtp(t)r(t+1)←Ax(t+1)−bβt+1←(p(t))TAp(t)(p(t))TAr(t+1)p(t+1)←−r(t+1)+βt+1p(t)k←k+1
其中, r ( t ) = ∇ ϕ ( x ( t ) ) r^{(t)} = \nabla \phi(x^{(t)}) r(t)=∇ϕ(x(t)) , r ( 0 ) = A x ( 0 ) − b , p ( 0 ) = − r ( 0 ) r^{(0)} = Ax^{(0)} - b, p^{(0)} = -r^{(0)} r(0)=Ax(0)−b,p(0)=−r(0)。
共轭梯度法只需要一阶导数信息,就可以计算步长和更新方向。收敛速度快,占用空间低,适用于求解大规模的线性方程。