矩阵分解及无约束最优化方法.docx

矩阵分解及无约束最优化方法.docx

《矩阵分解及无约束最优化方法.docx》由会员分享，可在线阅读，更多相关《矩阵分解及无约束最优化方法.docx（28页珍藏版）》请在冰豆网上搜索。

矩阵分解及无约束最优化方法

矩阵分解及无约束最优化

方法的原理和应用简介

最优化方法课程实验报告

学院：

数学与统计学院

班级：

硕2041班

姓名：

王彭

学号：

3112054028

指导教师：

阮小娥

同纟组人：

陈莹钱东东

矩阵分解及无约束最优化方法

的原理和应用简介

摘要

应课程学习的需要，本文主要对矩阵分解中的LU分解、LDLT分解、乔列斯基分解，以及无约束最优化领域中的最速下降法、牛顿法、拟牛顿法的原理、步骤和算法进行了简要介绍，并对各种方法进行了Matlab编程实验，得到了较好的结果。

关键字：

LU分解，LDLT分，、乔列斯基分解，最速下降法，牛顿法，拟牛顿法，Matlab编程。

【目录】

摘要-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-

232DFP法-13-

2.3.3BFGS法-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-

1矩阵分解

1.1矩阵的LU分解

1.1.1定义

若n阶矩阵A的各阶顺序主子式

1.1.2矩阵的LU分解过程

矩阵A的LU分解计算过程如下:

U1j,j1,2,,n

n1

Unn

ann

lnkUknk1

Lyb,

Uxy

来实现，而这两个方程组的求解只须前代和回带即可。

1.1.3矩阵LU分解的应用

1.1.3.1对矩阵进行LU分解

（1）问题描述

对矩阵

4

2

1

5

8

7

2

10

A

4

8

3

6

6

8

4

9

进行LU分解。

⑵用Matlab程序实现矩阵的

LU分解

Matlab程序见附录3.1结果为：

L=

1.0000

0

0

D

2.000Q

1.0000

0

0

l_0C00

2-0000

UOOOO

0

1.5000

1,^667

L2500

1.0000

U=

4.0000

2.0000

1.0000

5.0000

0

3.0000

0

0

0

0

2.0000

1.0000

0

0

0

C.2500

1.1.3.2LU分解在解方程组Axb中的应用

（1）问题描述

解方程组

Axb，

其中

4

2

1

5

A

8

7

2

10

4

8

3

6

6

8

4

9

b

[3

5

9

10]T

⑵用Matlab求解方程组AxbMatlab程序见附录

结果为：

30/3

-1/3

29/3-34/3

经验证，此解正确。

1.2对称矩阵的ldlt分解

1.2.1定义

设n阶矩阵A的各阶顺序主子式均不等于零，则存在惟一的单位下三角矩阵L，对角矩阵D和单位上三角矩阵MT，使得

ALDMT.

特别地，当A是对称矩阵时，ML,即矩阵A可以唯一地分解为

ALDLt，

其中L是单位下三角矩阵，D是对角阵。

1.2.2对称矩阵的LDLt分解过程

对称矩阵A的LDLT分解计算过程如下：

d1a11,

li1ai1/d1,i2,3,,n

forj2,3,,n1

j1

djajj1jkdk

k1

forij1,j2,,n

j1

1ij（aij1ikdk1jk）/dj

k1

n1

dnann1nkdk

k1

对矩阵A进行了LDLt分解后，求解线性方程组Axb，可以通过依次求解

F列三个方程组

Ly

b,

Dz

y,

Ltx

z

来实现。

1.2.3对称矩阵的LDLt分解应用

（1）问题描述

对矩阵

5

24

A2

12

4

25

进行ldlt分解。

⑵LDLt分解的Matlab程序实现

Matlab程序见附录3.2

实验结果为：

L=

1-0000

0

0

0.400Q

l.CQOO

0

-0-9000

-2.0000

1-0000

D=

5*0000

0

0

0

0.2000

0

0

0

1,0000

1.3对称正定矩阵的

ggt分解

1.3.1定义

设A是n阶对称正定矩阵，则存在一个可逆的下三角阵G，使得

AGGt.

当限定G的对角元为正时，这种分解时惟一的，称该分解为矩阵A的GGt分解

或乔列斯基（Cholesky）分解。

132对称正定矩阵的乔列斯基分解过程

对称正定矩阵A的乔列斯基分解计算过程如下：

J

g11aii

gnn

gi1ai1/gi1,i2,3,,n

for

j

2,3,

id

n1

gjj

（ajj

j1

2、1/2

gjk）

k1

for

ij

1,j

j1

2,,n

gj

（aij

gikgjk）/gjj

k1

（ann

n1

1/2

将正定矩阵A进行乔列斯基分解之后，求解线性方程组Axb，可以通过依次求解一下两个三角形方程组：

Gyb,

Gtxy

来实现。

133对称矩阵的乔列斯基分解应用

（1）问题描述对矩阵

425

进行乔列斯基分解。

⑵乔列斯基分解的Matlab程序实现

Matlab程序见附录3.3。

实验结果为：

2.236100

0.S9440.44720

-I.70S90*3705-0*759311.G&65

2无约束最优化方法

2.1最速下降法

2.1.1最速下降法的原理

最速下降法的搜索法向是目标函数的负梯度方向，最速下降法从目标函数的

负梯度方向一直前进，直到到达目标函数的最低点。

已知目标函数在X（k）点的梯度为：

当求目标函数的最小点时，由于函数沿负梯度方向下降最快，故在X（k）点的探索

显然，s（k）为单位向量。

这样第k1次迭代计算所得的新点为

负梯度仅给出了最优化方向，而没有给出步长的大小，所以可能有各种各样的最速下降的过程，它们依赖于

（k）

的大小。

（k）

步长有两种取法:

种方法是任意给定一个初始步长，使满足条件：

f（X（k）（k）s（k））f（x（k））

另外一种方法是沿负梯度方向做一维探索，以求解一维最优化问题的最优步长，即对目标函数极小，以得到最优步长：

（k）

。

minf（x（k）S（k））f（X（k）（k）S（k））

以此最优步长作为由x（k）点出发沿该点的负梯度方向探索的步长

这种方法的迭代计算的收敛性，可用以下三式中的任一式或二式作为准则来进行判断：

2.1.2最速下降法的步骤

用最速下降法求无约束多维极值问题minf（x）,xRn的算法步骤如下:

（1）取初始点x（0），精度0，令k0

（2）计算搜索方向v（k）f（x（k）），其中f（x（k））表示函数f（x）在点x（k）处的梯度；

（3）若IV"l，则停止计算；否则，从x（k）出发，沿v（k）进行一维搜索，割法等算法，当然也可以用MATLAB勺fminbnd函数;

（k1）（k）（k）

（4）令xx*,kk1，转步骤

（2）

2.1.3最速下降法的应用

2.1.3.1问题描述

22

用最速下降法求函数f（t,s）（t4）（s2）1极小值，取初始点取

x0（1,3）.

2.1.3.2用Matlab实现最速下降法求解函数极值

Matlab程序见附录

运行结果为：

4-0000-2.OOCO

mf=

1

2.2牛顿法

221牛顿法的原理

牛顿法是一种收敛很快的方法，其基本思路是利用二次曲线来逐点近似原目标函数，以二次曲线的极小值点来近似原目标函数的极小值点并逐渐逼近改点。

一维目标函数f（x）在x（k）点逼近用的二次曲线（即泰勒二次多项式）为

（（k）（k）（k）（k）1（k）（k）2

（x）f（x）f（x）（xx）-f（x）（xx）

2

此二次函数的极小点可由（x（k））0求得。

对于n维问题，n为目标函数f（X）在X（k）点逼近用的二次曲线为：

（X（k））f（x（k））f（X（k））.[XX（k）]^[XX（k）]T.2f（X（k））.[XX（k）]

2

令式中的Hessian2f（X（k））H（X（k）），则上式可改写为：

（X（k））f（x（k））f（X（k））.[XX（k）]-[XX（k）]T.H（X（k））.[XX（k）]

2

f（X）

当（X）0时可求得二次曲线（X）的极值点，且当且仅当改点处的

Hessian矩阵为正定时有极小值点。

由上式得：

（X）f（X（k））H（X（k））[XX（k）]

令（X）0，贝Uf（X（k））H（X（k））[XX（k）]0

1

若H（X（k））为可逆矩阵，将上式等号两边左乘H（X（k）），则得

H（X（k））1f（X（k））In[XX（k）]0

整理后得

1

XX（k）H（X（k））f（X（k））

当目标函数f（X）是二次函数时，牛顿法变得极为简单、有效，这时H（x（k））是一个常数矩阵，式

（X（k））f（x（k））f（x（k））.[XX（k）]2[Xx（k）]T.H（x（k））.[xX（k）]变成精

f（X）

1

确表达式，而利用式Xx（k）H（X（k））f（X（k））作一次迭代计算所得的X就

是最优点X*。

在一般情况下f（X）不一定是二次函数，则不能一步就能求出极小值，即极小值点不在H（X（k））f（X（k））方向上，但由于在X（k）点附近函数

（X）与f（X）是近似的，所以这个方向可以作为近似方向，可以用式

Xx（k）H（X（k））f（X（k））求出点X作为一个逼近点x（k1）。

这时式

1

Xx（k）H（X（k））f（X（k））可改成牛顿法的一般迭代公式：

X（k1）x（k）H（x（k））1f（x（k））

1

式中H（X（k））f（X（k））称为牛顿方向，通过这种迭代，逐次向极小值点X*逼近。

牛顿法无一维探索而直接代入公式计算，当初始点选得合适且当f（X）为二

次函数时收敛很快。

即使目标函数f（X）不是二次函数，当初始点选得好时，例如满足初始误差x（0）X*1时，也会很快收敛。

但是这种方法对初始点的选择

要求较高，如不满足x（0）X*1就不能保证有比较好的收敛性。

初始点的选择

有时甚至会影响到是否收敛。

如果选择不当使初始点离较大点近时，计算结果就

可能收敛于极大点。

初始点选择不当有时也会导致收敛到鞍点或不收敛，基于这

种原因，对古典的牛顿法要做些修改，于是便出现了修正牛顿法。

其修正方法是：

x（k）求x（k1）时不是直接用原来的迭代公式，而且沿着X（k）点处的牛顿方向进行一维搜索，将该方向上的最优点最为X（k1}0这样就会避免收敛于极大点或鞍点。

1

于是式X（k1）x（k）H（x（k））f（x（k））改写成：

X（k1）X（k）（k）H（X（k））1f（X（k））

S（k）H（X（k））1f（X（k））

式中的探索步长（k）为

minf（X（k）S（k））f（X（k）（k）S（k））

这种修正牛顿法虽然计算工作量多一些，但是具有收敛快的优点，并且，即使初始点选择不当，用这种探索方法也会成功。

2.2.2牛顿法的步骤

1

牛顿法是基于多元函数的泰勒展开而来的，它将H（X（k））f（X（k））作为

探索方向，因此它的迭代公式可直接写出来：

X（k1）X（k）H（x（k））1f（x（k））

牛顿法算法步骤如下所示：

（1）给定初始点x（0），及精度0，令k0；

（2）若|f（X（k））|，停止，极小点为x（k），否则转步骤（3）；

11

（3）计算2f（X（k）），令s（k）H（X（k））f（X（k））；

（4）令x（k1）x（k）s（k），kk1，转步骤

（2）。

2.2.3牛顿法的应用

2.2.3.1问题描述

220

用牛顿法求函数f（t，s）（t4）（s2）1极小值，取初始点取X（1,3）.

223.2用Matlab实现牛顿法求解函数极值

Matlab程序见附录3.5。

运行结果为：

1

2.3拟牛顿法

231拟牛顿法的原理

牛顿法的收敛速度虽然较快，但要求海森矩阵要可逆，要计算二阶导数和逆矩阵，就加大了就算机计算量。

为了克服牛顿法的缺点，同时保持较快收敛速度的优点，就产生了拟牛顿法。

拟牛顿法是牛顿法的直接推广，通过在试探点附近的二次逼近引进牛顿条件来确定线搜索方向，它主要有DFP和BFGS两种形式,

拟牛顿法的一般步骤为：

（1）给定初始点x（0），初始对称正定矩阵H。

，gog（x（0））及精度0;

（2）计算搜索方向p（k）Hkgk；

（3）作直线搜索x（k1）F（x（k）,p（k）），计算fkif（x（k1））,gkig（x（k1）），

Skx（k1）x（k）,ykgkigk

（4）判断终止准则是否满足；

（5）

令HkEk置kk1，转步骤

（2）；

2.3.2DFP法

2.3.2.1算法原理

DFP算法中的校正公式为:

为了保证Hk的正定性，在下面的算法中迭代一定次数后，重置初始点和迭代矩阵再进行迭代

2.3.1.2算法步骤

（1）给定初始点X（0），初始矩阵H。

In及精度0；

（2）若fx（0），停止，极小点为x（0）;否则转步骤（3）;

（3）取p（0）Hofx（0），且令k0；

（4）用一维搜索法求tk，使得tkf（X（k）tkp（k））minf（X（k）tp（k）），令

0

x（kdx（k）tp（k），转步骤（5）；

zIIfx（k1），停止，极小值点为x（k1）;否则转步骤（6）；

（5）11

（6）若k1n，令x（0）x（n），转步骤（3）；否则转步骤（7）；

Hk1

H

x

（k1）x（k）

zv

（k1）x

（k）x

T

k

（k

x

1）

（k）T

Xf

（k1）x

f

x

（k）

Hk

f

x

（k1）

fx（k）

（k

fx

1）

f

（k）T

xHk

f

（kx

1）

f

x（k）THk

f

（k1）

x

（k）

fx

（7）

令

取p（k）

Hk1

x（k1），置kk1，转步骤（4）

2.3.3BFGS法

2.3.3.1算法原理

BFGS算法中的校正公式为

1S（k）y（k）THHy（k）G（k）T

S（k）ty（k）SyHkHkyS

为了保证Hk的正定性，在下面算法步骤中迭代一定次数后，重置初始点和迭代矩阵再进行迭代。

2.3.3.2算法步骤

（1）给定初始点x（0），初始矩阵H°In及精度0；

（2）若fx（0），停止，极小点为x（0）;否则转步骤（3）；

（4）用一维搜索法求tk，使得tkf（X（k）tkp（k））minf（X（k）tp（k）），令

（5）

（4）。

X0（1,3）

2.342用Matlab实现拟牛顿法求解函数极值

Matlab程序见附录3.6。

运行结果为：

1.0000

-2.0000

L0000

3总结

通过本次实验，笔者对无约束最优化问题的原理及应用有了较深入的理解,对以后的进一步学习和科研工作有一定的指导价值。

但由于时间原因，本次实验与笔者的初始计划还是有一定的差距，例如笔者本打算对最速下降法、牛顿法、拟牛顿法等几种方法的异同从原理上进行分析，并给出实验验证，但由于时间原因，在本次报告中并未完成此部分内容，笔者希望在以后可以继续补上。

4附录

4.1矩阵LU分解的matlab程序：

（1）LU分解函数

function[LU]=LU1（A）

n=length（A（1,:

））;

L=zeros（n）;

U=zeros（n）;

fori=1:

n

L（i,i）=1;

U（1,i）=A（1,i）;

L（i,1）=A（i,1）/U（1,1）;

end

fori=2:

n

U（i,i）=A（i,i）-L（i,1:

i-1）*U（1:

i-1,i）;

forj=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=[4215;

87210;

4836;

6849];

[LU]=LU1（A）

⑶求解方程组Axb的Matlab主程序:

clc;

clear;

A=[4215;

87210;

4836;

6849];

[LU]=LU1（A）;

b=[35910]';

symsx1x2x3x4y1y2y3y4x;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=[y11y22y33y44]';

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=[x11x22x33x44]'

Ax=A*x

4.2对称矩阵的LDLT分解

4.2.1LDLT分解函数

function[LD]=MyLDLT（A）

n=length（A（1,:

））;

L=zeros（n）;

D=zeros（n）;

D（1,1）=A（1,1）;

fori=1:

n

L（i,i）=1;

end

fori=2:

n

L（i,1）=A（i,1）/D（1,1）;

end

summ=0;

forj=2:

n-1

fork=1:

j-1summ=summ+L（j,kF2*D（k,k）;

end

D（j,j）=A（j,j）-summ;

summm=0;

fori=j+1:

n

fork=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;

fori=1:

n-1

summ=summ+L（n,iF2*D（i,i）;end

D（n,n）=A（n,n）-summ;

4.2.2LDLT分解函数主程序

clc;

clear;

A=[52-4;21-2;-4-25]

[LD]=MyLDLT（A）

4.3正定举证的乔列斯基分解

4.3.1乔列斯基分解函数functionG=MyChol（A）n=length（A（1,:

））;G=zeros（n）;

G（1,1）=sqrt（A（1,1））;fori=2:

n

G（i,1）=A（i,1）/G（1,1）;

end

forj=2:

n-1

Ga,j）=（A（j,j）-G（j,1：

j-1）即（1/2）;

fori=j+1:

n

G（i,j）=（A（i,j）-G（i,1:

j-1FG（j,1：

j-1））A（1/2）;end

end

G（n,n）=（A（n,n）-G（n,1:

n-1）*G（n,1:

n-1）'F（1/2）;

4.3.2乔列斯基分解主程序

clc;

clear;

A=[52-4;21-2;-4-25]

G=MyChol（A）

4.4最速下降法

4.4.1最速下降法函数

function[x,minf]=minFD（f,x0,var,eps）%目标函数:

f;

%初始点:

x0;%自变量向量:

var;%精度:

eps;%目标函数取最小值时的自变量值:

x;

%目标函数的最小值:

minf;

formatlong;

ifnargin==3eps=1.0e-6;

end

symsl;

tol=1;

whiletol>eps

gradf=-jacobian（f,var）;%负梯度方向