矩阵分解及无约束最优化 方法的原理和应用简介
——最优化方法课程实验报告
学 院：数学与统计学院 班 级：硕2041班 姓 名：王彭 学 号：指导教师：阮小娥 同 组 人：陈莹 钱东东
矩阵分解及无约束最优化方法的原理和应用简介
矩阵分解及无约束最优化方法
的原理和应用简介
摘要
应课程学习的需要，本文主要对矩阵分解中的LU分解、LDLT分解、乔列斯基分解，以及无约束最优化领域中的最速下降法、牛顿法、拟牛顿法的原理、步骤和算法进行了简要介绍，并对各种方法进行了Matlab编程实验，得到了较好的结果。
关键字：LU分解，LDLT分，、乔列斯基分解，最速下降法，牛顿法，拟牛顿法，Matlab编程。
- 1 -
《最优化方法》课程实验报告
【目录】
摘要 ............................................................................................................................................. - 1 - 1 矩阵分解 ................................................................................................................................. - 3 -
1.1 矩阵的LU分解 ........................................................................................................... - 3 -
1.1.1 定义 .................................................................................................................... - 3 - 1.1.2 矩阵的LU分解过程 ........................................................................................ - 3 - 1.1.3 矩阵LU分解的应用 ........................................................................................ - 4 - 1.2 对称矩阵的LDL分解 ................................................................................................ - 5 -
1.2.1 定义 .................................................................................................................... - 5 - 1.2.2对称矩阵的LDL分解过程 ............................................................................... - 5 - 1.2.3对称矩阵的LDL分解应用 ............................................................................... - 6 - 1.3 对称正定矩阵的GG分解 ........................................................................................... - 6 -
1.3.1 定义 .................................................................................................................... - 6 - 1.3.2 对称正定矩阵的乔列斯基分解过程 ................................................................ - 7 - 1.3.3对称矩阵的乔列斯基分解应用 ......................................................................... - 7 -
2 无约束最优化方法.................................................................................................................. - 8 -
2.1 最速下降法 ................................................................................................................... - 8 -
2.1.1 最速下降法的原理 ............................................................................................ - 8 - 2.1.2 最速下降法的步骤 ............................................................................................ - 9 - 2.1.3 最速下降法的应用 ............................................................................................ - 9 - 2.2 牛顿法 ......................................................................................................................... - 10 -
2.2.1 牛顿法的原理 .................................................................................................. - 10 - 2.2.2 牛顿法的步骤 .................................................................................................. - 12 - 2.2.3 牛顿法的应用 .................................................................................................. - 12 - 2.3 拟牛顿法 ..................................................................................................................... - 13 -
2.3.1 拟牛顿法的原理 .............................................................................................. - 13 - 2.3.2 DFP法 .............................................................................................................. - 13 - 2.3.3 BFGS法 ............................................................................................................ - 14 - 2.3.4 拟牛顿法的应用 .............................................................................................. - 15 -
3 总结 ....................................................................................................................................... - 15 - 4 附录 ....................................................................................................................................... - 16 -
4.1 矩阵LU分解的matlab程序： ................................................................................ - 16 - 4.2 对称矩阵的LDL分解 .............................................................................................. - 17 - 4.3 正定举证的乔列斯基分解 ......................................................................................... - 18 - 4.4 最速下降法 ................................................................................................................. - 18 - 4.5 牛顿法 ......................................................................................................................... - 19 - 4.6 拟牛顿法 ..................................................................................................................... - 20 -
- 2 -
矩阵分解及无约束最优化方法的原理和应用简介
1 矩阵分解
1.1 矩阵的LU分解
1.1.1 定义
若n阶矩阵A的各阶顺序主子式
a11a21 ai1
a12a22 ai2
a1i a2i
≠0, aii
i=1,2, ,n，
则存在惟一的单位下三角矩阵L和可逆的上三角阵U，满足
A=LU，
称该式为矩阵A的LU分解。
1.1.2 矩阵的LU分解过程
矩阵A的LU分解计算过程如下：
u1j=a1j,li1=ai1/u11,
j=1,2, ,ni=2,3, ,n
⎡fori=2,3, ,n-1
i-1⎢
⎢uii=aii-∑likuki
k=1⎢
⎢⎡
⎢⎢forj=i+1,i+2, ,n⎢⎢i-1
⎢⎢u=a-lu
∑ijijikkj
⎢⎢k=1
i-1⎢⎢
⎢⎢lji=(aji-∑ljkuki)/uii⎢⎢k=1⎣⎢⎣
unn=ann-∑lnkukn
k=1n-1
若对矩阵A进行了LU分解，求解线性方程组Ax=b，可以通过依次求解以下两个三角方程组：
Ly=b,
Ux=y
来实现，而这两个方程组的求解只须前代和回带即可。
- 3 -
《最优化方法》课程实验报告
1.1.3 矩阵LU分解的应用
1.1.3.1 对矩阵进行LU分解 (1)问题描述
对矩阵
⎡4⎢8A=⎢
⎢427815⎤210⎥⎥ 36⎥
⎢⎣6
8
49⎥⎦
进行LU分解。
(2)用Matlab程序实现矩阵的LU分解
Matlab程序见附录3.1。 结果为：
1.1.3.2 LU分解在解方程组Ax=b中的应用 (1)问题描述
解方程组
Ax=b，
其中
⎡⎢4215⎤A=⎢
87210⎥⎢⎢4836⎥⎥
； ⎣6
8
49⎥⎦
b=[35910]T. (2)用Matlab求解方程组Ax=b
Matlab程序见附录 结果为：
- 4 -
矩阵分解及无约束最优化方法的原理和应用简介
经验证，此解正确。
1.2 对称矩阵的LDLT分解
1.2.1 定义
设n阶矩阵A的各阶顺序主子式均不等于零，则存在惟一的单位下三角矩阵
L，对角矩阵D和单位上三角矩阵MT，使得
A=LDMT.
特别地，当A是对称矩阵时，M=L，即矩阵A可以唯一地分解为
A=LDLT，
其中L是单位下三角矩阵，D是对角阵。
1.2.2对称矩阵的LDLT分解过程
对称矩阵A的LDLT分解计算过程如下：
d1=a11,
li1=ai1/d1,i=2,3, ,n⎡
⎢⎢
⎢forj=2,3, ,n-1
j-1
⎢
dj=ajj-∑l2⎢jkdk
k=1⎢
⎢⎡fori=j+1,j+2, ,n
j-1⎢⎢
⎢⎢lij=(aij-∑likdkljk)/dj
k=1⎣⎣
2
dn=ann-∑lnkdk
k=1n-1
对矩阵A进行了LDLT分解后，求解线性方程组Ax=b，可以通过依次求解下列三个方程组
- 5 -
《最优化方法》课程实验报告
Ly=b,Dz=y, LTx=z
来实现。
1.2.3对称矩阵的LDLT分解应用
(1)问题描述
对矩阵
2-4⎤⎡5
⎥ A=⎢21-2⎢⎥
⎢⎣-4-25⎥⎦
进行LDLT分解。
(2) LDLT分解的Matlab程序实现 Matlab程序见附录3.2。
实验结果为：
1.3 对称正定矩阵的GGT分解
1.3.1 定义
设A是n阶对称正定矩阵，则存在一个可逆的下三角阵G，使得
A=GGT.
当限定G的对角元为正时，这种分解时惟一的，称该分解为矩阵A的GGT分解或乔列斯基(Cholesky)分解。
- 6 -
矩阵分解及无约束最优化方法的原理和应用简介
1.3.2 对称正定矩阵的乔列斯基分解过程
对称正定矩阵A的乔列斯基分解计算过程如下：
g11=a11
gi1=ai1/g11,i=2,3, ,n⎡⎢⎢
⎢forj=2,3, ,n-1
j-1
⎢21/2 ⎢gjj=(ajj-∑gjk)
k=1⎢
⎢⎡fori=j+1,j+2, ,n
j-1⎢⎢
⎢⎢gij=(aij-∑gikgjk)/gjj
k=1⎣⎣
21/2
gnn=(ann-∑gnk)
k=1n-1
将正定矩阵A进行乔列斯基分解之后，求解线性方程组Ax=b，可以通过
依次求解一下两个三角形方程组：
Gy=b,GTx=y
来实现。
1.3.3对称矩阵的乔列斯基分解应用
(1)问题描述
对矩阵
2-4⎤⎡5
⎥ A=⎢21-2⎢⎥
⎢⎣-4-25⎥⎦
进行乔列斯基分解。
(2)乔列斯基分解的Matlab程序实现 Matlab程序见附录3.3。 实验结果为：
- 7 -
《最优化方法》课程实验报告
2 无约束最优化方法
2.1 最速下降法
2.1.1 最速下降法的原理
最速下降法的搜索法向是目标函数的负梯度方向，最速下降法从目标函数的负梯度方向一直前进，直到到达目标函数的最低点。
已知目标函数在X
(k)
点的梯度为：
∂f(X(k))∂x2
∂f(X(k))⎤
⎥...
∂xn⎥
⎦
(k)
⎡∂f(X(k))
∇f(X(k))=⎢
∂x1⎢⎣
T
当求目标函数的最小点时，由于函数沿负梯度方向下降最快，故在X方向应取该点的负梯度方向，即
点的探索
S
(k)
=-
∇f(X(k))∇fX(k)
(k)
S显然，为单位向量。这样第k+1次迭代计算所得的新点为
X(k+1)=X(k)+α(k)S(k)=X(k)-
α(k)∇f(X(k))
∇fX(k)
负梯度仅给出了最优化方向，而没有给出步长的大小，所以可能有各种各样的最速下降的过程，它们依赖于
α(k)
∇fX(k)的大小。
步长α
(k)
有两种取法：
一种方法是任意给定一个初始步长，使满足条件：
f(X(k)+α(k)S(k))
另外一种方法是沿负梯度方向做一维探索，以求解一维最优化问题的最优步长α，即对目标函数极小，以得到最优步长：
minf(X(k)+αS(k))=f(X(k)+α(k)S(k))
α>0
(k)
(k)
以此最优步长作为由X点出发沿该点的负梯度方向探索的步长α
。
- 8 -
矩阵分解及无约束最优化方法的原理和应用简介
这种方法的迭代计算的收敛性，可用以下三式中的任一式或二式作为准则来进行判断：
⎧∇f(X(k))≤ε1⎪
⎪fX(k)-fX(k-1)
)()⎪(≤ε2⎨(k)fX⎪
⎪(k)
(k-1)
≤ε3⎪X-X
⎩
2.1.2 最速下降法的步骤
n
minf(x),x∈R用最速下降法求无约束多维极值问题的算法步骤如下： (0)
（1）取初始点x，精度ε>0，令k=0
（2）计算搜索方向v的梯度；
（3）若
v(k)≤ε
(k)
=-∇f(x(k))，其中∇f(x(k))表示函数f(x)在点x(k)处
，则停止计算；否则，从x
(k)
出发，沿v
(k)
进行一维搜索，
f(x(k)+λkv(k))=minf(x(k)+λv(k))λλ≥0即求k，使得。此处的一维搜索可以用黄金分
割法等算法，当然也可以用MATLAB的fminbnd函数；
(k+1)(k)(k)x=x+λv,k=k+1，转步骤（2）k（4）令。
2.1.3 最速下降法的应用
2.1.3.1 问题描述
22
f(t,s)=(t-4)+(s+2)+1极小值，取初始点取用最速下降法求函数
x0=(1,-3).
2.1.3.2 用Matlab实现最速下降法求解函数极值
Matlab程序见附录 运行结果为：
- 9 -
《最优化方法》课程实验报告
2.2 牛顿法
2.2.1 牛顿法的原理
牛顿法是一种收敛很快的方法，其基本思路是利用二次曲线来逐点近似原目标函数，以二次曲线的极小值点来近似原目标函数的极小值点并逐渐逼近改点。
一维目标函数f(x) 在x(k)点逼近用的二次曲线（即泰勒二次多项式）为
ϕ(x(k))=f(x(k))+f'(x(k))(x-x(k))+
此二次函数的极小点可由ϕ'(x(k))=0求得。
1
f''(x(k))(x-x(k))2 2
对于n维问题，n为目标函数f(X)在X(k)点逼近用的二次曲线为：
(k)(k)(k)T2(k)(k)
⎤ϕ(X(k))=f(x(k))+⎡∇f(X).[X-X]+[X-X].∇f(X).[X-X] ⎣⎦
1
2
令式中的Hessian∇2f(X(k))=H(X(k))，则上式可改写为：
(k)(k)(k)T(k)(k)
⎤ϕ(X(k))=f(x(k))+⎡∇f(X).[X-X]+[X-X].H(X).[X-X]⎣⎦
1
2
≈f(X)
当∇ϕ(X)=0时可求得二次曲线ϕ(X)的极值点，且当且仅当改点处的
Hessian矩阵为正定时有极小值点。
由上式得：
∇ϕ(X)=∇f(X(k))+H(X(k))[X-X(k)]
令∇ϕ(X)=0，则∇f(X(k))+H(X(k))[X-X(k)]=0
(k)
H(X)⎤若H(X(k))为可逆矩阵，将上式等号两边左乘⎡⎣⎦，则得
(k)(k)(k)
⎡⎤H(X)∇f(X)+I[X-X]=0 n⎣⎦
-1
-1
整理后得
- 10 -
矩阵分解及无约束最优化方法的原理和应用简介
(k)(k)
⎤X=X(k)-⎡H(X)∇f(X) ⎣⎦
-1
当目标函数f(X)是二次函数时，牛顿法变得极为简单、有效，这时H(X(k))是一个常数矩阵，式
(k)(k)(k)T(k)(k)
⎤ϕ(X(k))=f(x(k))+⎡∇f(X).[X-X]+[X-X].H(X).[X-X]⎣⎦
1
2
变成精
≈f(X)
(k)(k)
⎤H(X)∇f(X)作一次迭代计算所得的X就确表达式，而利用式X=X(k)-⎡⎣⎦
-1
是最优点X*。在一般情况下f(X)不一定是二次函数，则不能一步就能求出极小
(k)(k)(k)
⎤XH(X)∇f(X)值，即极小值点不在-⎡方向上，但由于在点附近函数⎣⎦
-1
ϕ(X)与f(X)是近似的，所以这个方向可以作为近似方向，可以用式
(k+1)(k)(k)
⎤XXX=X(k)-⎡H(X)∇f(X)求出点作为一个逼近点。这时式⎣⎦
-1
(k)k()⎤X=X(k)-⎡H(X)∇f(X)可改成牛顿法的一般迭代公式： ⎣⎦
-1
X
-1
(k+1)
=X
(k)
(k)(k)
-⎡⎣H(X)⎤⎦∇f(X)
-1
*(k)(k)
⎤XH(X)∇f(X)式中-⎡称为牛顿方向，通过这种迭代，逐次向极小值点逼⎣⎦
近。
牛顿法无一维探索而直接代入公式计算，当初始点选得合适且当f(X)为二次函数时收敛很快。即使目标函数f(X)不是二次函数，当初始点选得好时，例如满足初始误差X(0)-X*
可能收敛于极大点。初始点选择不当有时也会导致收敛到鞍点或不收敛，基于这种原因，对古典的牛顿法要做些修改，于是便出现了修正牛顿法。其修正方法是：
X(k)求X(k+1)时不是直接用原来的迭代公式，而且沿着X(k)点处的牛顿方向进行
一维搜索，将该方向上的最优点最为X(k+1)。这样就会避免收敛于极大点或鞍点。于是式X
(k+1)
=X
(k)
(k)(k)
-⎡⎣H(X)⎤⎦∇f(X)改写成：
-1
- 11 -
《最优化方法》课程实验报告
(k)(k)
⎤X(k+1)=X(k)-α(k)⎡H(X)∇f(X) ⎣⎦
-1
令
(k)(k)
⎤S(k)=-⎡H(X)∇f(X) ⎣⎦
-1
则
X(k+1)=X(k)+α(k)S(k)
式中的探索步长α(k)为
minf(X(k)+αS(k))=f(X(k)+α(k)S(k))
α≥0
这种修正牛顿法虽然计算工作量多一些，但是具有收敛快的优点，并且，即
使初始点选择不当，用这种探索方法也会成功。
2.2.2 牛顿法的步骤
(k)(k)
⎤H(X)∇f(X)作为牛顿法是基于多元函数的泰勒展开而来的，它将-⎡⎣⎦
-1
探索方向，因此它的迭代公式可直接写出来：
X
(k+1)
=X
(k)
(k)(k)
-⎡⎣H(X)⎤⎦∇f(X)
-1
牛顿法算法步骤如下所示：
（1）给定初始点x(0)，及精度ε>0，令k=0；
（2）若∇f(X(k))≤ε，停止，极小点为x(k)，否则转步骤（3）；
2(k)(k)(k)(k)
（3）计算⎡⎣∇f(X)⎤⎦，令s=-⎡⎣H(X)⎤⎦∇f(X)；
-1
-1
（4）令x(k+1)=x(k)+s(k)，k=k+1，转步骤（2）。
2.2.3 牛顿法的应用
2.2.3.1 问题描述
220
用牛顿法求函数f(t,s)=(t-4)+(s+2)+1极小值，取初始点取x=(1,-3).
2.2.3.2 用Matlab实现牛顿法求解函数极值
Matlab程序见附录3.5。 运行结果为：
- 12 -
矩阵分解及无约束最优化方法的原理和应用简介
2.3 拟牛顿法
2.3.1 拟牛顿法的原理
牛顿法的收敛速度虽然较快，但要求海森矩阵要可逆，要计算二阶导数和逆矩阵，就加大了就算机计算量。为了克服牛顿法的缺点，同时保持较快收敛速度的优点，就产生了拟牛顿法。拟牛顿法是牛顿法的直接推广，通过在试探点附近的二次逼近引进牛顿条件来确定线搜索方向，它主要有DFP和BFGS两种形式，拟牛顿法的一般步骤为：
（1）给定初始点x(0)，初始对称正定矩阵H0，g0=g(x(0))及精度ε>0； （2）计算搜索方向p(k)=-Hkgk；
（3）作直线搜索x(k+1)=F(x(k),p(k))，计算fk+1=f(x(k+1)),gk+1=g(x(k+1))，
Sk=x(k+1)-x(k),yk=gk+1-gk
（4）判断终止准则是否满足；
（5）令Hk+1=Hk+Ek置k=k+1，转步骤（2）；
不同的拟牛顿法对应不同的Ek，主要介绍DFP和BFGS两种拟牛顿法。
2.3.2 DFP法
2.3.2.1算法原理
DFP算法中的校正公式为：
Hk+1=Hk+
SkSkTSkTyk
-
HkykyTHk
k
yTHkykk
为了保证Hk的正定性，在下面的算法中迭代一定次数后，重置初始点和迭代矩阵再进行迭代。 2.3.1.2算法步骤
- 13 -
《最优化方法》课程实验报告
（1）给定初始点x(0)，初始矩阵H0=In及精度ε>0；
（2）若∇f(x(0))≤ε，停止，极小点为x(0)；否则转步骤（3）； （3）取p(0)=-H0∇f(x(0))，且令k=0；
(（4）用一维搜索法求tk，使得tkf(X(k)+tkp(k))=minfX
α≥0
()
，转步骤（5）； x(k+1)=x(k)+tpk
(k)
+tp
k()
)，令
（5）
∇f(x(k+1))≤ε，停止，极小值点为x(k+1)；否则转步骤（6）；
（6）若k+1=n，令x(0)=x(n)，转步骤（3）；否则转步骤（7）；
Hk+1=Hk+
(x(k+1)-x(k))(x(k+1)-x(k))(x
(k+1)
T
)∇f(x)-∇f(x)（7）令
H(∇f(x)-∇f(x))(∇f(x)-∇f(x))H-
∇f(x)-∇f(x)H∇f(x)-∇f(x)-x
(k)T
(k+1)
(k)
(k+1)
(k)
(k+1)
(k)
T
k
(k+1)
(k)
T
(k+1)
(k)
k
，
k
取p(k)=-Hk+1∇f(x(k+1))，置k=k+1，转步骤（4）。
2.3.3 BFGS法
2.3.3.1算法原理
BFGS算法中的校正公式为
TT
S(k)(S(k))⎡(y(k))Hky(k)⎤
⎢1+⎥Hk+1=Hk+TT(k)(k)⎥Sy(S(k))y(k)⎢()⎣⎦
-
1
(S(k))y(k)
T
⎡S(k)(y(k))TH+Hy(k)(S(k))T⎤
kk⎢⎥⎣⎦
为了保证Hk的正定性，在下面算法步骤中迭代一定次数后，重置初始点和迭代矩阵再进行迭代。
2.3.3.2算法步骤
（1）给定初始点x(0)，初始矩阵H0=In及精度ε>0；
（2）若∇f(x(0))≤ε，停止，极小点为x(0)；否则转步骤（3）； （3）取p(0)=-H0∇f(x(0))，且令k=0；
- 14 -
矩阵分解及无约束最优化方法的原理和应用简介
(（4）用一维搜索法求tk，使得tkf(X(k)+tkp(k))=minfX
α≥0
()
，转步骤（5）； x(k+1)=x(k)+tpk
(k)
+tp
k()
)，令
（5）
∇f(x(k+1))≤ε，停止，极小值点为x(k+1)；否则转步骤（6）；
（6）若k+1=n，令x(0)=x(n)，转步骤（3）；否则转步骤（7）； （7）令
TT
S(k)(S(k))⎡(y(k))Hky(k)⎤
⎢1+⎥Hk+1=Hk+TT(k)(k)⎥Sy(S(k))y(k)⎢()⎣⎦
-
1
(S(k))y(k)
T
⎡S(k)(y(k))TH+Hy(k)(S(k))T⎤
kk⎢⎥⎣⎦
其中：
S(k)=x(k+1)-x(k)y
(k)
=∇f(x
(k+1)
)-∇f(x)
(k)
+)
，取p(k)=-Hk+1∇f(x(k1)，置k=k+1，转步骤
（4）。
2.3.4 拟牛顿法的应用
2.3.4.1 问题描述
22
用拟牛顿法求函数f(t,s)=(t-4)+(s+2)+1极小值，取初始点取
x0=(1,-3).
2.3.4.2 用Matlab实现拟牛顿法求解函数极值
Matlab程序见附录3.6
。 运行结果为：
3 总结
通过本次实验，笔者对无约束最优化问题的原理及应用有了较深入的理解，
- 15 -
《最优化方法》课程实验报告
对以后的进一步学习和科研工作有一定的指导价值。但由于时间原因，本次实验与笔者的初始计划还是有一定的差距，例如笔者本打算对最速下降法、牛顿法、拟牛顿法等几种方法的异同从原理上进行分析，并给出实验验证，但由于时间原因，在本次报告中并未完成此部分内容，笔者希望在以后可以继续补上。
4 附录
4.1 矩阵LU分解的matlab程序：
(1) LU分解函数
function [L U]=LU1(A) n=length(A(1,:)); L=zeros(n); U=zeros(n); for i=1:n L(i,i)=1; U(1,i)=A(1,i); L(i,1)=A(i,1)/U(1,1); end for i=2:n
U(i,i)=A(i,i)-L(i,1:i-1)*U(1:i-1,i); for j=i+1:n
U(i,j)=A(i,j)-L(i,1:i-1)*U(1:i-1,j);
L(j,i)=(A(j,i)-L(j,1:i-1)*U(1:i-1,i))/U(i,i); end end
(2) 对矩阵进行LU分解主程序：
clc; clear; A=[4 2 1 5; 8 7 2 10; 4 8 3 6; 6 8 4 9]; [L U]=LU1(A)
(3) 求解方程组Ax b的Matlab主程序：
clc; clear;
A=[4 2 1 5; 8 7 2 10; 4 8 3 6; 6 8 4 9]; [L U]=LU1(A); b=[3 5 9 10]';
- 16 -
矩阵分解及无约束最优化方法的原理和应用简介
syms x1 x2 x3 x4 y1 y2 y3 y4 x; y11=solve(L(1,1)*y1-b(1),y1);
y22=solve(L(2,1)*y11+L(2,2)*y2-b(2),y2);
y33=solve(L(3,1)*y11+L(3,2)*y22+L(3,3)*y3-b(3),y3);
y44=solve(L(4,1)*y11+L(4,2)*y22+L(4,3)*y33+L(4,4)*y4-b(4),y4); y=[y11 y22 y33 y44]';
x44=solve(U(4,4)*x4-y(4),x4);
x33=solve(U(3,3)*x3+U(3,4)*x44-y(3),x3);
x22=solve(U(2,2)*x2+U(2,3)*x33+U(2,4)*x44-y(2),x2);
x11=solve(U(1,1)*x1+U(1,2)*x22+U(1,3)*x33+U(1,4)*x44-y(1),x1); x=[x11 x22 x33 x44]' Ax=A*x
4.2 对称矩阵的LDLT分解
4.2.1 LDL分解函数
function [L D]=MyLDLT(A) n=length(A(1,:)); L=zeros(n); D=zeros(n); D(1,1)=A(1,1); for i=1:n L(i,i)=1; end for i=2:n
L(i,1)=A(i,1)/D(1,1); end summ=0;
for j=2:n-1 for k=1:j-1
summ=summ+L(j,k)^2*D(k,k); end
D(j,j)=A(j,j)-summ; summm=0; for i=j+1:n for k=1:j-1
summm=summm+L(i,k)*D(k,k)*L(j,k); end
L(i,j)=(A(i,j)-summm)/D(j,j); end end summ=0;
- 17 -
T
《最优化方法》课程实验报告
for i=1:n-1
summ=summ+L(n,i)^2*D(i,i); end
D(n,n)=A(n,n)-summ;
4.2.2 LDLT分解函数主程序
clc; clear;
A=[5 2 -4;2 1 -2;-4 -2 5] [L D]=MyLDLT(A)
4.3 正定举证的乔列斯基分解
4.3.1 乔列斯基分解函数
function G=MyChol(A) n=length(A(1,:)); G=zeros(n);
G(1,1)=sqrt(A(1,1)); for i=2:n
G(i,1)=A(i,1)/G(1,1); end
for j=2:n-1
G(j,j)=(A(j,j)-G(j,1:j-1)^2)^(1/2); for i=j+1:n
G(i,j)=(A(i,j)-G(i,1:j-1)^G(j,1:j-1))^(1/2); end end
G(n,n)=(A(n,n)-G(n,1:n-1)*G(n,1:n-1)')^(1/2);
4.3.2 乔列斯基分解主程序
clc; clear;
A=[5 2 -4;2 1 -2;-4 -2 5] G=MyChol(A)
4.4 最速下降法
4.4.1 最速下降法函数
function [x,minf]=minFD(f,x0,var,eps) %目标函数:f; %初始点:x0; %自变量向量:var; %精度:eps;
%目标函数取最小值时的自变量值:x; %目标函数的最小值:minf;
- 18 -
矩阵分解及无约束最优化方法的原理和应用简介
format long; if nargin==3 eps=1.0e-6; end syms l; tol=1; while tol>eps
gradf=-jacobian(f,var);%负梯度方向 v=Funval(gradf,var,x0); tol=norm(v); y=x0+l*v;
yf=Funval(f,var,y); [a,b]=minJT(yf,0,0.1);
xm=minHJ(yf,a,b);%用黄金分割法进行一维搜索 x1=x0+xm*v; x0=x1; end x=x1;
minf=Funval(f,var,x); format short;
4.4.2 最速下降法求解极值问题主函数
clc; clear; syms t s;
f=(t-4)^2+(s+2)^2+1;
[x,mf]=minFD(f,[1 -3],[t s])
4.5 牛顿法
4.5.1 牛顿法函数
function [x,minf]=minNT(f,x0,var,eps) %目标函数:f; %初始点:x0; %自变量向量:var;
%目标函数取最小值时的自变量值:x; %目标函数的最小值:minf; format long; if nargin==3 eps=1.0e-6; end tol=1;
x0=transpose(x0); while tol>eps
gradf=jacobian(f,var); %梯度方向
- 19 -
《最优化方法》课程实验报告
jacf=jacobian(gradf,var); %雅克比矩阵
v=Funval(gradf,var,x0);
tol=norm(v);
pv=Funval(jacf,var,x0);
p=-inv(pv)*transpose(v);
x1=x0+p;
x0=x1;
end
x=x1;
minf=Funval(f,var,x);
format short;
4.5.2 牛顿法求解极值问题主函数
clc;
clear;
syms t s;
f=(t-4)^2+(s+2)^2+1;
[x,mf]=minNT(f,[1000 -300000000],[t s])
4.6 拟牛顿法
4.6.1 DFP拟牛顿法函数
function [x,minf]=minDFP(f,x0,var,eps)
%目标函数:f;
%初始点:x0;
%自变量向量:var;
%目标函数取最小值时的自变量值:x;
%目标函数的最小值:minf;
format long;
if nargin==3
eps=1.0e-6;
end
x0=transpose(x0);
n=length(var);
syms l;
H=eye(n,n);
gradf=jacobian(f,var);
v0=Funval(gradf,var,x0);
p=-H*transpose(v0);
k=0;
while 1
v=Funval(gradf,var,x0);
tol=norm(v);
if tol
x=x0;
- 20 -
矩阵分解及无约束最优化方法的原理和应用简介
break;
end
y=x0+l*p;
yf=Funval(f,var,y);
[a,b]=minJT(yf,0,0.1);
xm=minHJ(yf,a,b);
x1=x0+xm*p;
vk=Funval(gradf,var,x1);
tol=norm(vk);
if tol
x=x1;
break;
end
if k+1==n
x0=x1;
continue;
else
dx=x1-x0;
dgf=vk-v;
dgf=transpose(dgf);
dxT=transpose(dx);
dgfT=transpose(dgf);
mdx=dx*dxT;
mdgf=dgf*dgfT;
fz=H*(dgf*(dgfT*H));
H=H+mdx/(dxT*dgf)-inv(dgfT*(H*dgf))*fz;
p=-H*transpose(vk);
k=k+1;
x0=x1;
end
end
minf=Funval(f,var,x);
format short;
4.6.2 BFGS拟牛顿法函数
function [x,minf]=minBFGS(f,x0,var,eps)
%目标函数:f;
%初始点:x0;
%自变量向量:var;
%目标函数取最小值时的自变量值:x;
%目标函数的最小值:minf;
format long;
if nargin==3
eps=1.0e-6;
end
- 21 -
《最优化方法》课程实验报告
x0=transpose(x0);
n=length(var);
syms l;
H=eye(n,n);
gradf=jacobian(f,var);
v0=Funval(gradf,var,x0);
p=-H*transpose(v0);
k=0;
while 1
v=Funval(gradf,var,x0);
tol=norm(v);
if tol
x=x0;
break;
end
y=x0+l*p;
yf=Funval(f,var,y);
[a,b]=minJT(yf,0,0.1);
xm=minHJ(yf,a,b);
x1=x0+xm*p;
vk=Funval(gradf,var,x1);
tol=norm(vk);
if tol
x=x1;
break;
end
if k+1==n
x0=x1;
continue;
else
dx=x1-x0;
dgf=vk-v;
dgf=transpose(dgf);
dxT=transpose(dx);
dgfT=transpose(dgf);
mdx=dx*dxT;
mdgf=dgf*dgfT;
H=H+(1+dgfT*(H*dgf)/(dxT*dgf))*mdx/(dxT*dgf)-(dx*dgfT*H+H*dgf*dxT)/(dxT*dgf);
p=-H*transpose(vk);
k=k+1;
x0=x1;
end
- 22 -
矩阵分解及无约束最优化方法的原理和应用简介
end
minf=Funval(f,var,x);
format short;
4.6.3 拟牛顿法求解极值问题主函数
clc;
clear;
syms t s;
f=(t-4)^2+(s+2)^2+1;
[x,mf]=minDFP(f,[1 -3],[t s])
- 23 -
矩阵分解及无约束最优化 方法的原理和应用简介
——最优化方法课程实验报告
学 院：数学与统计学院 班 级：硕2041班 姓 名：王彭 学 号：指导教师：阮小娥 同 组 人：陈莹 钱东东
矩阵分解及无约束最优化方法的原理和应用简介
矩阵分解及无约束最优化方法
的原理和应用简介
摘要
应课程学习的需要，本文主要对矩阵分解中的LU分解、LDLT分解、乔列斯基分解，以及无约束最优化领域中的最速下降法、牛顿法、拟牛顿法的原理、步骤和算法进行了简要介绍，并对各种方法进行了Matlab编程实验，得到了较好的结果。
关键字：LU分解，LDLT分，、乔列斯基分解，最速下降法，牛顿法，拟牛顿法，Matlab编程。
- 1 -
《最优化方法》课程实验报告
【目录】
摘要 ............................................................................................................................................. - 1 - 1 矩阵分解 ................................................................................................................................. - 3 -
1.1 矩阵的LU分解 ........................................................................................................... - 3 -
1.1.1 定义 .................................................................................................................... - 3 - 1.1.2 矩阵的LU分解过程 ........................................................................................ - 3 - 1.1.3 矩阵LU分解的应用 ........................................................................................ - 4 - 1.2 对称矩阵的LDL分解 ................................................................................................ - 5 -
1.2.1 定义 .................................................................................................................... - 5 - 1.2.2对称矩阵的LDL分解过程 ............................................................................... - 5 - 1.2.3对称矩阵的LDL分解应用 ............................................................................... - 6 - 1.3 对称正定矩阵的GG分解 ........................................................................................... - 6 -
1.3.1 定义 .................................................................................................................... - 6 - 1.3.2 对称正定矩阵的乔列斯基分解过程 ................................................................ - 7 - 1.3.3对称矩阵的乔列斯基分解应用 ......................................................................... - 7 -
2 无约束最优化方法.................................................................................................................. - 8 -
2.1 最速下降法 ................................................................................................................... - 8 -
2.1.1 最速下降法的原理 ............................................................................................ - 8 - 2.1.2 最速下降法的步骤 ............................................................................................ - 9 - 2.1.3 最速下降法的应用 ............................................................................................ - 9 - 2.2 牛顿法 ......................................................................................................................... - 10 -
2.2.1 牛顿法的原理 .................................................................................................. - 10 - 2.2.2 牛顿法的步骤 .................................................................................................. - 12 - 2.2.3 牛顿法的应用 .................................................................................................. - 12 - 2.3 拟牛顿法 ..................................................................................................................... - 13 -
2.3.1 拟牛顿法的原理 .............................................................................................. - 13 - 2.3.2 DFP法 .............................................................................................................. - 13 - 2.3.3 BFGS法 ............................................................................................................ - 14 - 2.3.4 拟牛顿法的应用 .............................................................................................. - 15 -
3 总结 ....................................................................................................................................... - 15 - 4 附录 ....................................................................................................................................... - 16 -
4.1 矩阵LU分解的matlab程序： ................................................................................ - 16 - 4.2 对称矩阵的LDL分解 .............................................................................................. - 17 - 4.3 正定举证的乔列斯基分解 ......................................................................................... - 18 - 4.4 最速下降法 ................................................................................................................. - 18 - 4.5 牛顿法 ......................................................................................................................... - 19 - 4.6 拟牛顿法 ..................................................................................................................... - 20 -
- 2 -
矩阵分解及无约束最优化方法的原理和应用简介
1 矩阵分解
1.1 矩阵的LU分解
1.1.1 定义
若n阶矩阵A的各阶顺序主子式
a11a21 ai1
a12a22 ai2
a1i a2i
≠0, aii
i=1,2, ,n，
则存在惟一的单位下三角矩阵L和可逆的上三角阵U，满足
A=LU，
称该式为矩阵A的LU分解。
1.1.2 矩阵的LU分解过程
矩阵A的LU分解计算过程如下：
u1j=a1j,li1=ai1/u11,
j=1,2, ,ni=2,3, ,n
⎡fori=2,3, ,n-1
i-1⎢
⎢uii=aii-∑likuki
k=1⎢
⎢⎡
⎢⎢forj=i+1,i+2, ,n⎢⎢i-1
⎢⎢u=a-lu
∑ijijikkj
⎢⎢k=1
i-1⎢⎢
⎢⎢lji=(aji-∑ljkuki)/uii⎢⎢k=1⎣⎢⎣
unn=ann-∑lnkukn
k=1n-1
若对矩阵A进行了LU分解，求解线性方程组Ax=b，可以通过依次求解以下两个三角方程组：
Ly=b,
Ux=y
来实现，而这两个方程组的求解只须前代和回带即可。
- 3 -
《最优化方法》课程实验报告
1.1.3 矩阵LU分解的应用
1.1.3.1 对矩阵进行LU分解 (1)问题描述
对矩阵
⎡4⎢8A=⎢
⎢427815⎤210⎥⎥ 36⎥
⎢⎣6
8
49⎥⎦
进行LU分解。
(2)用Matlab程序实现矩阵的LU分解
Matlab程序见附录3.1。 结果为：
1.1.3.2 LU分解在解方程组Ax=b中的应用 (1)问题描述
解方程组
Ax=b，
其中
⎡⎢4215⎤A=⎢
87210⎥⎢⎢4836⎥⎥
； ⎣6
8
49⎥⎦
b=[35910]T. (2)用Matlab求解方程组Ax=b
Matlab程序见附录 结果为：
- 4 -
矩阵分解及无约束最优化方法的原理和应用简介
经验证，此解正确。
1.2 对称矩阵的LDLT分解
1.2.1 定义
设n阶矩阵A的各阶顺序主子式均不等于零，则存在惟一的单位下三角矩阵
L，对角矩阵D和单位上三角矩阵MT，使得
A=LDMT.
特别地，当A是对称矩阵时，M=L，即矩阵A可以唯一地分解为
A=LDLT，
其中L是单位下三角矩阵，D是对角阵。
1.2.2对称矩阵的LDLT分解过程
对称矩阵A的LDLT分解计算过程如下：
d1=a11,
li1=ai1/d1,i=2,3, ,n⎡
⎢⎢
⎢forj=2,3, ,n-1
j-1
⎢
dj=ajj-∑l2⎢jkdk
k=1⎢
⎢⎡fori=j+1,j+2, ,n
j-1⎢⎢
⎢⎢lij=(aij-∑likdkljk)/dj
k=1⎣⎣
2
dn=ann-∑lnkdk
k=1n-1
对矩阵A进行了LDLT分解后，求解线性方程组Ax=b，可以通过依次求解下列三个方程组
- 5 -
《最优化方法》课程实验报告
Ly=b,Dz=y, LTx=z
来实现。
1.2.3对称矩阵的LDLT分解应用
(1)问题描述
对矩阵
2-4⎤⎡5
⎥ A=⎢21-2⎢⎥
⎢⎣-4-25⎥⎦
进行LDLT分解。
(2) LDLT分解的Matlab程序实现 Matlab程序见附录3.2。
实验结果为：
1.3 对称正定矩阵的GGT分解
1.3.1 定义
设A是n阶对称正定矩阵，则存在一个可逆的下三角阵G，使得
A=GGT.
当限定G的对角元为正时，这种分解时惟一的，称该分解为矩阵A的GGT分解或乔列斯基(Cholesky)分解。
- 6 -
矩阵分解及无约束最优化方法的原理和应用简介
1.3.2 对称正定矩阵的乔列斯基分解过程
对称正定矩阵A的乔列斯基分解计算过程如下：
g11=a11
gi1=ai1/g11,i=2,3, ,n⎡⎢⎢
⎢forj=2,3, ,n-1
j-1
⎢21/2 ⎢gjj=(ajj-∑gjk)
k=1⎢
⎢⎡fori=j+1,j+2, ,n
j-1⎢⎢
⎢⎢gij=(aij-∑gikgjk)/gjj
k=1⎣⎣
21/2
gnn=(ann-∑gnk)
k=1n-1
将正定矩阵A进行乔列斯基分解之后，求解线性方程组Ax=b，可以通过
依次求解一下两个三角形方程组：
Gy=b,GTx=y
来实现。
1.3.3对称矩阵的乔列斯基分解应用
(1)问题描述
对矩阵
2-4⎤⎡5
⎥ A=⎢21-2⎢⎥
⎢⎣-4-25⎥⎦
进行乔列斯基分解。
(2)乔列斯基分解的Matlab程序实现 Matlab程序见附录3.3。 实验结果为：
- 7 -
《最优化方法》课程实验报告
2 无约束最优化方法
2.1 最速下降法
2.1.1 最速下降法的原理
最速下降法的搜索法向是目标函数的负梯度方向，最速下降法从目标函数的负梯度方向一直前进，直到到达目标函数的最低点。
已知目标函数在X
(k)
点的梯度为：
∂f(X(k))∂x2
∂f(X(k))⎤
⎥...
∂xn⎥
⎦
(k)
⎡∂f(X(k))
∇f(X(k))=⎢
∂x1⎢⎣
T
当求目标函数的最小点时，由于函数沿负梯度方向下降最快，故在X方向应取该点的负梯度方向，即
点的探索
S
(k)
=-
∇f(X(k))∇fX(k)
(k)
S显然，为单位向量。这样第k+1次迭代计算所得的新点为
X(k+1)=X(k)+α(k)S(k)=X(k)-
α(k)∇f(X(k))
∇fX(k)
负梯度仅给出了最优化方向，而没有给出步长的大小，所以可能有各种各样的最速下降的过程，它们依赖于
α(k)
∇fX(k)的大小。
步长α
(k)
有两种取法：
一种方法是任意给定一个初始步长，使满足条件：
f(X(k)+α(k)S(k))
另外一种方法是沿负梯度方向做一维探索，以求解一维最优化问题的最优步长α，即对目标函数极小，以得到最优步长：
minf(X(k)+αS(k))=f(X(k)+α(k)S(k))
α>0
(k)
(k)
以此最优步长作为由X点出发沿该点的负梯度方向探索的步长α
。
- 8 -
矩阵分解及无约束最优化方法的原理和应用简介
这种方法的迭代计算的收敛性，可用以下三式中的任一式或二式作为准则来进行判断：
⎧∇f(X(k))≤ε1⎪
⎪fX(k)-fX(k-1)
)()⎪(≤ε2⎨(k)fX⎪
⎪(k)
(k-1)
≤ε3⎪X-X
⎩
2.1.2 最速下降法的步骤
n
minf(x),x∈R用最速下降法求无约束多维极值问题的算法步骤如下： (0)
（1）取初始点x，精度ε>0，令k=0
（2）计算搜索方向v的梯度；
（3）若
v(k)≤ε
(k)
=-∇f(x(k))，其中∇f(x(k))表示函数f(x)在点x(k)处
，则停止计算；否则，从x
(k)
出发，沿v
(k)
进行一维搜索，
f(x(k)+λkv(k))=minf(x(k)+λv(k))λλ≥0即求k，使得。此处的一维搜索可以用黄金分
割法等算法，当然也可以用MATLAB的fminbnd函数；
(k+1)(k)(k)x=x+λv,k=k+1，转步骤（2）k（4）令。
2.1.3 最速下降法的应用
2.1.3.1 问题描述
22
f(t,s)=(t-4)+(s+2)+1极小值，取初始点取用最速下降法求函数
x0=(1,-3).
2.1.3.2 用Matlab实现最速下降法求解函数极值
Matlab程序见附录 运行结果为：
- 9 -
《最优化方法》课程实验报告
2.2 牛顿法
2.2.1 牛顿法的原理
牛顿法是一种收敛很快的方法，其基本思路是利用二次曲线来逐点近似原目标函数，以二次曲线的极小值点来近似原目标函数的极小值点并逐渐逼近改点。
一维目标函数f(x) 在x(k)点逼近用的二次曲线（即泰勒二次多项式）为
ϕ(x(k))=f(x(k))+f'(x(k))(x-x(k))+
此二次函数的极小点可由ϕ'(x(k))=0求得。
1
f''(x(k))(x-x(k))2 2
对于n维问题，n为目标函数f(X)在X(k)点逼近用的二次曲线为：
(k)(k)(k)T2(k)(k)
⎤ϕ(X(k))=f(x(k))+⎡∇f(X).[X-X]+[X-X].∇f(X).[X-X] ⎣⎦
1
2
令式中的Hessian∇2f(X(k))=H(X(k))，则上式可改写为：
(k)(k)(k)T(k)(k)
⎤ϕ(X(k))=f(x(k))+⎡∇f(X).[X-X]+[X-X].H(X).[X-X]⎣⎦
1
2
≈f(X)
当∇ϕ(X)=0时可求得二次曲线ϕ(X)的极值点，且当且仅当改点处的
Hessian矩阵为正定时有极小值点。
由上式得：
∇ϕ(X)=∇f(X(k))+H(X(k))[X-X(k)]
令∇ϕ(X)=0，则∇f(X(k))+H(X(k))[X-X(k)]=0
(k)
H(X)⎤若H(X(k))为可逆矩阵，将上式等号两边左乘⎡⎣⎦，则得
(k)(k)(k)
⎡⎤H(X)∇f(X)+I[X-X]=0 n⎣⎦
-1
-1
整理后得
- 10 -
矩阵分解及无约束最优化方法的原理和应用简介
(k)(k)
⎤X=X(k)-⎡H(X)∇f(X) ⎣⎦
-1
当目标函数f(X)是二次函数时，牛顿法变得极为简单、有效，这时H(X(k))是一个常数矩阵，式
(k)(k)(k)T(k)(k)
⎤ϕ(X(k))=f(x(k))+⎡∇f(X).[X-X]+[X-X].H(X).[X-X]⎣⎦
1
2
变成精
≈f(X)
(k)(k)
⎤H(X)∇f(X)作一次迭代计算所得的X就确表达式，而利用式X=X(k)-⎡⎣⎦
-1
是最优点X*。在一般情况下f(X)不一定是二次函数，则不能一步就能求出极小
(k)(k)(k)
⎤XH(X)∇f(X)值，即极小值点不在-⎡方向上，但由于在点附近函数⎣⎦
-1
ϕ(X)与f(X)是近似的，所以这个方向可以作为近似方向，可以用式
(k+1)(k)(k)
⎤XXX=X(k)-⎡H(X)∇f(X)求出点作为一个逼近点。这时式⎣⎦
-1
(k)k()⎤X=X(k)-⎡H(X)∇f(X)可改成牛顿法的一般迭代公式： ⎣⎦
-1
X
-1
(k+1)
=X
(k)
(k)(k)
-⎡⎣H(X)⎤⎦∇f(X)
-1
*(k)(k)
⎤XH(X)∇f(X)式中-⎡称为牛顿方向，通过这种迭代，逐次向极小值点逼⎣⎦
近。
牛顿法无一维探索而直接代入公式计算，当初始点选得合适且当f(X)为二次函数时收敛很快。即使目标函数f(X)不是二次函数，当初始点选得好时，例如满足初始误差X(0)-X*
可能收敛于极大点。初始点选择不当有时也会导致收敛到鞍点或不收敛，基于这种原因，对古典的牛顿法要做些修改，于是便出现了修正牛顿法。其修正方法是：
X(k)求X(k+1)时不是直接用原来的迭代公式，而且沿着X(k)点处的牛顿方向进行
一维搜索，将该方向上的最优点最为X(k+1)。这样就会避免收敛于极大点或鞍点。于是式X
(k+1)
=X
(k)
(k)(k)
-⎡⎣H(X)⎤⎦∇f(X)改写成：
-1
- 11 -
《最优化方法》课程实验报告
(k)(k)
⎤X(k+1)=X(k)-α(k)⎡H(X)∇f(X) ⎣⎦
-1
令
(k)(k)
⎤S(k)=-⎡H(X)∇f(X) ⎣⎦
-1
则
X(k+1)=X(k)+α(k)S(k)
式中的探索步长α(k)为
minf(X(k)+αS(k))=f(X(k)+α(k)S(k))
α≥0
这种修正牛顿法虽然计算工作量多一些，但是具有收敛快的优点，并且，即
使初始点选择不当，用这种探索方法也会成功。
2.2.2 牛顿法的步骤
(k)(k)
⎤H(X)∇f(X)作为牛顿法是基于多元函数的泰勒展开而来的，它将-⎡⎣⎦
-1
探索方向，因此它的迭代公式可直接写出来：
X
(k+1)
=X
(k)
(k)(k)
-⎡⎣H(X)⎤⎦∇f(X)
-1
牛顿法算法步骤如下所示：
（1）给定初始点x(0)，及精度ε>0，令k=0；
（2）若∇f(X(k))≤ε，停止，极小点为x(k)，否则转步骤（3）；
2(k)(k)(k)(k)
（3）计算⎡⎣∇f(X)⎤⎦，令s=-⎡⎣H(X)⎤⎦∇f(X)；
-1
-1
（4）令x(k+1)=x(k)+s(k)，k=k+1，转步骤（2）。
2.2.3 牛顿法的应用
2.2.3.1 问题描述
220
用牛顿法求函数f(t,s)=(t-4)+(s+2)+1极小值，取初始点取x=(1,-3).
2.2.3.2 用Matlab实现牛顿法求解函数极值
Matlab程序见附录3.5。 运行结果为：
- 12 -
矩阵分解及无约束最优化方法的原理和应用简介
2.3 拟牛顿法
2.3.1 拟牛顿法的原理
牛顿法的收敛速度虽然较快，但要求海森矩阵要可逆，要计算二阶导数和逆矩阵，就加大了就算机计算量。为了克服牛顿法的缺点，同时保持较快收敛速度的优点，就产生了拟牛顿法。拟牛顿法是牛顿法的直接推广，通过在试探点附近的二次逼近引进牛顿条件来确定线搜索方向，它主要有DFP和BFGS两种形式，拟牛顿法的一般步骤为：
（1）给定初始点x(0)，初始对称正定矩阵H0，g0=g(x(0))及精度ε>0； （2）计算搜索方向p(k)=-Hkgk；
（3）作直线搜索x(k+1)=F(x(k),p(k))，计算fk+1=f(x(k+1)),gk+1=g(x(k+1))，
Sk=x(k+1)-x(k),yk=gk+1-gk
（4）判断终止准则是否满足；
（5）令Hk+1=Hk+Ek置k=k+1，转步骤（2）；
不同的拟牛顿法对应不同的Ek，主要介绍DFP和BFGS两种拟牛顿法。
2.3.2 DFP法
2.3.2.1算法原理
DFP算法中的校正公式为：
Hk+1=Hk+
SkSkTSkTyk
-
HkykyTHk
k
yTHkykk
为了保证Hk的正定性，在下面的算法中迭代一定次数后，重置初始点和迭代矩阵再进行迭代。 2.3.1.2算法步骤
- 13 -
《最优化方法》课程实验报告
（1）给定初始点x(0)，初始矩阵H0=In及精度ε>0；
（2）若∇f(x(0))≤ε，停止，极小点为x(0)；否则转步骤（3）； （3）取p(0)=-H0∇f(x(0))，且令k=0；
(（4）用一维搜索法求tk，使得tkf(X(k)+tkp(k))=minfX
α≥0
()
，转步骤（5）； x(k+1)=x(k)+tpk
(k)
+tp
k()
)，令
（5）
∇f(x(k+1))≤ε，停止，极小值点为x(k+1)；否则转步骤（6）；
（6）若k+1=n，令x(0)=x(n)，转步骤（3）；否则转步骤（7）；
Hk+1=Hk+
(x(k+1)-x(k))(x(k+1)-x(k))(x
(k+1)
T
)∇f(x)-∇f(x)（7）令
H(∇f(x)-∇f(x))(∇f(x)-∇f(x))H-
∇f(x)-∇f(x)H∇f(x)-∇f(x)-x
(k)T
(k+1)
(k)
(k+1)
(k)
(k+1)
(k)
T
k
(k+1)
(k)
T
(k+1)
(k)
k
，
k
取p(k)=-Hk+1∇f(x(k+1))，置k=k+1，转步骤（4）。
2.3.3 BFGS法
2.3.3.1算法原理
BFGS算法中的校正公式为
TT
S(k)(S(k))⎡(y(k))Hky(k)⎤
⎢1+⎥Hk+1=Hk+TT(k)(k)⎥Sy(S(k))y(k)⎢()⎣⎦
-
1
(S(k))y(k)
T
⎡S(k)(y(k))TH+Hy(k)(S(k))T⎤
kk⎢⎥⎣⎦
为了保证Hk的正定性，在下面算法步骤中迭代一定次数后，重置初始点和迭代矩阵再进行迭代。
2.3.3.2算法步骤
（1）给定初始点x(0)，初始矩阵H0=In及精度ε>0；
（2）若∇f(x(0))≤ε，停止，极小点为x(0)；否则转步骤（3）； （3）取p(0)=-H0∇f(x(0))，且令k=0；
- 14 -
矩阵分解及无约束最优化方法的原理和应用简介
(（4）用一维搜索法求tk，使得tkf(X(k)+tkp(k))=minfX
α≥0
()
，转步骤（5）； x(k+1)=x(k)+tpk
(k)
+tp
k()
)，令
（5）
∇f(x(k+1))≤ε，停止，极小值点为x(k+1)；否则转步骤（6）；
（6）若k+1=n，令x(0)=x(n)，转步骤（3）；否则转步骤（7）； （7）令
TT
S(k)(S(k))⎡(y(k))Hky(k)⎤
⎢1+⎥Hk+1=Hk+TT(k)(k)⎥Sy(S(k))y(k)⎢()⎣⎦
-
1
(S(k))y(k)
T
⎡S(k)(y(k))TH+Hy(k)(S(k))T⎤
kk⎢⎥⎣⎦
其中：
S(k)=x(k+1)-x(k)y
(k)
=∇f(x
(k+1)
)-∇f(x)
(k)
+)
，取p(k)=-Hk+1∇f(x(k1)，置k=k+1，转步骤
（4）。
2.3.4 拟牛顿法的应用
2.3.4.1 问题描述
22
用拟牛顿法求函数f(t,s)=(t-4)+(s+2)+1极小值，取初始点取
x0=(1,-3).
2.3.4.2 用Matlab实现拟牛顿法求解函数极值
Matlab程序见附录3.6
。 运行结果为：
3 总结
通过本次实验，笔者对无约束最优化问题的原理及应用有了较深入的理解，
- 15 -
《最优化方法》课程实验报告
对以后的进一步学习和科研工作有一定的指导价值。但由于时间原因，本次实验与笔者的初始计划还是有一定的差距，例如笔者本打算对最速下降法、牛顿法、拟牛顿法等几种方法的异同从原理上进行分析，并给出实验验证，但由于时间原因，在本次报告中并未完成此部分内容，笔者希望在以后可以继续补上。
4 附录
4.1 矩阵LU分解的matlab程序：
(1) LU分解函数
function [L U]=LU1(A) n=length(A(1,:)); L=zeros(n); U=zeros(n); for i=1:n L(i,i)=1; U(1,i)=A(1,i); L(i,1)=A(i,1)/U(1,1); end for i=2:n
U(i,i)=A(i,i)-L(i,1:i-1)*U(1:i-1,i); for j=i+1:n
U(i,j)=A(i,j)-L(i,1:i-1)*U(1:i-1,j);
L(j,i)=(A(j,i)-L(j,1:i-1)*U(1:i-1,i))/U(i,i); end end
(2) 对矩阵进行LU分解主程序：
clc; clear; A=[4 2 1 5; 8 7 2 10; 4 8 3 6; 6 8 4 9]; [L U]=LU1(A)
(3) 求解方程组Ax b的Matlab主程序：
clc; clear;
A=[4 2 1 5; 8 7 2 10; 4 8 3 6; 6 8 4 9]; [L U]=LU1(A); b=[3 5 9 10]';
- 16 -
矩阵分解及无约束最优化方法的原理和应用简介
syms x1 x2 x3 x4 y1 y2 y3 y4 x; y11=solve(L(1,1)*y1-b(1),y1);
y22=solve(L(2,1)*y11+L(2,2)*y2-b(2),y2);
y33=solve(L(3,1)*y11+L(3,2)*y22+L(3,3)*y3-b(3),y3);
y44=solve(L(4,1)*y11+L(4,2)*y22+L(4,3)*y33+L(4,4)*y4-b(4),y4); y=[y11 y22 y33 y44]';
x44=solve(U(4,4)*x4-y(4),x4);
x33=solve(U(3,3)*x3+U(3,4)*x44-y(3),x3);
x22=solve(U(2,2)*x2+U(2,3)*x33+U(2,4)*x44-y(2),x2);
x11=solve(U(1,1)*x1+U(1,2)*x22+U(1,3)*x33+U(1,4)*x44-y(1),x1); x=[x11 x22 x33 x44]' Ax=A*x
4.2 对称矩阵的LDLT分解
4.2.1 LDL分解函数
function [L D]=MyLDLT(A) n=length(A(1,:)); L=zeros(n); D=zeros(n); D(1,1)=A(1,1); for i=1:n L(i,i)=1; end for i=2:n
L(i,1)=A(i,1)/D(1,1); end summ=0;
for j=2:n-1 for k=1:j-1
summ=summ+L(j,k)^2*D(k,k); end
D(j,j)=A(j,j)-summ; summm=0; for i=j+1:n for k=1:j-1
summm=summm+L(i,k)*D(k,k)*L(j,k); end
L(i,j)=(A(i,j)-summm)/D(j,j); end end summ=0;
- 17 -
T
《最优化方法》课程实验报告
for i=1:n-1
summ=summ+L(n,i)^2*D(i,i); end
D(n,n)=A(n,n)-summ;
4.2.2 LDLT分解函数主程序
clc; clear;
A=[5 2 -4;2 1 -2;-4 -2 5] [L D]=MyLDLT(A)
4.3 正定举证的乔列斯基分解
4.3.1 乔列斯基分解函数
function G=MyChol(A) n=length(A(1,:)); G=zeros(n);
G(1,1)=sqrt(A(1,1)); for i=2:n
G(i,1)=A(i,1)/G(1,1); end
for j=2:n-1
G(j,j)=(A(j,j)-G(j,1:j-1)^2)^(1/2); for i=j+1:n
G(i,j)=(A(i,j)-G(i,1:j-1)^G(j,1:j-1))^(1/2); end end
G(n,n)=(A(n,n)-G(n,1:n-1)*G(n,1:n-1)')^(1/2);
4.3.2 乔列斯基分解主程序
clc; clear;
A=[5 2 -4;2 1 -2;-4 -2 5] G=MyChol(A)
4.4 最速下降法
4.4.1 最速下降法函数
function [x,minf]=minFD(f,x0,var,eps) %目标函数:f; %初始点:x0; %自变量向量:var; %精度:eps;
%目标函数取最小值时的自变量值:x; %目标函数的最小值:minf;
- 18 -
矩阵分解及无约束最优化方法的原理和应用简介
format long; if nargin==3 eps=1.0e-6; end syms l; tol=1; while tol>eps
gradf=-jacobian(f,var);%负梯度方向 v=Funval(gradf,var,x0); tol=norm(v); y=x0+l*v;
yf=Funval(f,var,y); [a,b]=minJT(yf,0,0.1);
xm=minHJ(yf,a,b);%用黄金分割法进行一维搜索 x1=x0+xm*v; x0=x1; end x=x1;
minf=Funval(f,var,x); format short;
4.4.2 最速下降法求解极值问题主函数
clc; clear; syms t s;
f=(t-4)^2+(s+2)^2+1;
[x,mf]=minFD(f,[1 -3],[t s])
4.5 牛顿法
4.5.1 牛顿法函数
function [x,minf]=minNT(f,x0,var,eps) %目标函数:f; %初始点:x0; %自变量向量:var;
%目标函数取最小值时的自变量值:x; %目标函数的最小值:minf; format long; if nargin==3 eps=1.0e-6; end tol=1;
x0=transpose(x0); while tol>eps
gradf=jacobian(f,var); %梯度方向
- 19 -
《最优化方法》课程实验报告
jacf=jacobian(gradf,var); %雅克比矩阵
v=Funval(gradf,var,x0);
tol=norm(v);
pv=Funval(jacf,var,x0);
p=-inv(pv)*transpose(v);
x1=x0+p;
x0=x1;
end
x=x1;
minf=Funval(f,var,x);
format short;
4.5.2 牛顿法求解极值问题主函数
clc;
clear;
syms t s;
f=(t-4)^2+(s+2)^2+1;
[x,mf]=minNT(f,[1000 -300000000],[t s])
4.6 拟牛顿法
4.6.1 DFP拟牛顿法函数
function [x,minf]=minDFP(f,x0,var,eps)
%目标函数:f;
%初始点:x0;
%自变量向量:var;
%目标函数取最小值时的自变量值:x;
%目标函数的最小值:minf;
format long;
if nargin==3
eps=1.0e-6;
end
x0=transpose(x0);
n=length(var);
syms l;
H=eye(n,n);
gradf=jacobian(f,var);
v0=Funval(gradf,var,x0);
p=-H*transpose(v0);
k=0;
while 1
v=Funval(gradf,var,x0);
tol=norm(v);
if tol
x=x0;
- 20 -
矩阵分解及无约束最优化方法的原理和应用简介
break;
end
y=x0+l*p;
yf=Funval(f,var,y);
[a,b]=minJT(yf,0,0.1);
xm=minHJ(yf,a,b);
x1=x0+xm*p;
vk=Funval(gradf,var,x1);
tol=norm(vk);
if tol
x=x1;
break;
end
if k+1==n
x0=x1;
continue;
else
dx=x1-x0;
dgf=vk-v;
dgf=transpose(dgf);
dxT=transpose(dx);
dgfT=transpose(dgf);
mdx=dx*dxT;
mdgf=dgf*dgfT;
fz=H*(dgf*(dgfT*H));
H=H+mdx/(dxT*dgf)-inv(dgfT*(H*dgf))*fz;
p=-H*transpose(vk);
k=k+1;
x0=x1;
end
end
minf=Funval(f,var,x);
format short;
4.6.2 BFGS拟牛顿法函数
function [x,minf]=minBFGS(f,x0,var,eps)
%目标函数:f;
%初始点:x0;
%自变量向量:var;
%目标函数取最小值时的自变量值:x;
%目标函数的最小值:minf;
format long;
if nargin==3
eps=1.0e-6;
end
- 21 -
《最优化方法》课程实验报告
x0=transpose(x0);
n=length(var);
syms l;
H=eye(n,n);
gradf=jacobian(f,var);
v0=Funval(gradf,var,x0);
p=-H*transpose(v0);
k=0;
while 1
v=Funval(gradf,var,x0);
tol=norm(v);
if tol
x=x0;
break;
end
y=x0+l*p;
yf=Funval(f,var,y);
[a,b]=minJT(yf,0,0.1);
xm=minHJ(yf,a,b);
x1=x0+xm*p;
vk=Funval(gradf,var,x1);
tol=norm(vk);
if tol
x=x1;
break;
end
if k+1==n
x0=x1;
continue;
else
dx=x1-x0;
dgf=vk-v;
dgf=transpose(dgf);
dxT=transpose(dx);
dgfT=transpose(dgf);
mdx=dx*dxT;
mdgf=dgf*dgfT;
H=H+(1+dgfT*(H*dgf)/(dxT*dgf))*mdx/(dxT*dgf)-(dx*dgfT*H+H*dgf*dxT)/(dxT*dgf);
p=-H*transpose(vk);
k=k+1;
x0=x1;
end
- 22 -
矩阵分解及无约束最优化方法的原理和应用简介
end
minf=Funval(f,var,x);
format short;
4.6.3 拟牛顿法求解极值问题主函数
clc;
clear;
syms t s;
f=(t-4)^2+(s+2)^2+1;
[x,mf]=minDFP(f,[1 -3],[t s])
- 23 -