矩阵代数数值计算.ppt

《矩阵代数数值计算.ppt》由会员分享，可在线阅读，更多相关《矩阵代数数值计算.ppt（66页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第五章第五章矩阵代数数值计算矩阵代数数值计算一、一、矩阵的基本运算矩阵的基本运算二、二、矩阵的三角分解矩阵的三角分解三、三、矩阵的正交变换矩阵的正交变换四、四、矩阵的谱分解矩阵的谱分解五、五、IMSL中的线性系统、特征中的线性系统、特征值分析模块值分析模块矩矩阵阵代代数数运运算算是是统统计计模模型型的的基基础础，统统计计模模型型的的所所有有估估计计几几乎乎都都是是用用矩矩阵阵代代数数运运算算计计算算出出结结果果。例例如如最最小小二二乘乘估估计计、典典型型相相关关分分析析、因因子子分分析析以以及及各各类类回回归归分分析析。从从计计算算的的角角度度来来说说为为使使计计算算结结果果可可靠靠，我我们们

2、总总是是先先对对矩矩阵阵进进行行三三角角分分解解，然然后后进进行行各各种种计计算算例例如如，矩阵的逆、求解线性方程组以及对矩阵进行谱分解等。矩阵的逆、求解线性方程组以及对矩阵进行谱分解等。本本章章首首先先介介绍绍矩矩阵阵的的三三角角分分解解，然然后后引引导导学学习习者者使使用用IMSL和和SASD中中的的丰丰富富矩矩阵阵的的算算法法，将将它它们们拼拼接接起来就可以解决各种矩阵的计算。起来就可以解决各种矩阵的计算。5.1引言引言n矩矩阵阵代代数数运运算算在在数数值值计计算算中中起起着着基基础础性性的的作作用用，无无论论我我们们建建立立了了多多么么复复杂杂的的数数学学模模型型，最最终终我我们们总总

3、是是要要把把它它变变为为矩矩阵阵代代数数的的形形式式。特特别别是是统统计计模模型型，无无论论是是多多元元线线性性回回归归、广广义义线线性性模模型型、多多元元统统计计分分析析无无不不与与矩矩阵阵代代数数有有着着密密切切的的联联系系。我我们们所所研研究究的的对对象象，即即样样本本可可以以看看成成是是一一个个矩矩阵阵，而而统统计计上上的的协协方方差差，实实际际上上是是该该矩矩阵阵的的转转置置与与该该矩矩阵阵相相乘乘形形成成的的新矩阵的元素。而回归的最小二乘估计的算法为：新矩阵的元素。而回归的最小二乘估计的算法为：（5.1.1)包包含含了了矩矩阵阵的的乘乘、矩矩阵阵的的求求逆逆及及矩矩阵阵与与向向量量

4、的的乘乘法法等等等等。而而特特征征值值与与特特征征向向量量在在数数理理统统计计理理都都有有明明确确的的条条件件统统计计含含义义。因因此此我我们们将将在在这这一一章章介介绍绍矩矩阵阵运运算算的的基基本本数数值值方方法法，以以及及如如何何调用调用IMSL和和SASD中丰富的算法模块。中丰富的算法模块。5.2矩阵数值计算基础矩阵数值计算基础n对于一般的二维样本，我们总可以写成如下的矩阵形式。对于一般的二维样本，我们总可以写成如下的矩阵形式。（5.2.1)从计算的角度处理矩阵问题的一个最有效的方法是，将一个一从计算的角度处理矩阵问题的一个最有效的方法是，将一个一般的矩阵分解为几个简单矩阵的乘积形式。其

5、中最便于计算的是三般的矩阵分解为几个简单矩阵的乘积形式。其中最便于计算的是三角矩阵，以上三角矩阵为例，由其性质，两个上三角阵的和、积仍角矩阵，以上三角矩阵为例，由其性质，两个上三角阵的和、积仍为上三角阵，上三角阵的特征值就是其对角线元素。为上三角阵，上三角阵的特征值就是其对角线元素。系数矩阵为上三角阵的线性线性方程组是最容易求解的，系数矩阵为上三角阵的线性线性方程组是最容易求解的，上三角阵的逆阵仍然是上三角阵。因此处理矩阵计算问题的关上三角阵的逆阵仍然是上三角阵。因此处理矩阵计算问题的关键是将一般矩阵化为三角阵和对角阵的形式，然后进行计算。键是将一般矩阵化为三角阵和对角阵的形式，然后进行计算。

6、5.2.1矩阵的三角矩阵的三角三角分解三角分解（1）L*R分解分解（5.2.2)其中其中L单位下三角阵（主对角线元素为单位下三角阵（主对角线元素为1），），R为上三角阵。为上三角阵。（2）LDR*分解分解(5.2.3)其中其中L为单位下三角阵，为单位下三角阵，D为对角阵，为单位上三角阵。为对角阵，为单位上三角阵。（3）Crout分解分解（5.2.4)其中为单位下三角阵，为单位上三角阵。其中为单位下三角阵，为单位上三角阵。(4)Cholesky分解分解（5.2.5)其中其中A为正定对称阵，为正定对称阵，T 为上三角阵。为上三角阵。Cholesky Cholesky 分解是统计计分解是统计计算中最

7、常用的分解方法之一。因为我们的协方差矩阵、相关矩算中最常用的分解方法之一。因为我们的协方差矩阵、相关矩阵都是使用这种分解方法。阵都是使用这种分解方法。5.2.2矩阵的三角分解算法矩阵的三角分解算法以上四种分解是类似的，使用待定系数法。以上四种分解是类似的，使用待定系数法。(1)以以LR*分解为例，设分解为例，设其中其中A为正定阵，并记分解已为以下形式：为正定阵，并记分解已为以下形式：=利用矩阵相对应元素相等的事实，我们立即得到利用矩阵相对应元素相等的事实，我们立即得到现在我们可以计算矩阵现在我们可以计算矩阵L的第一列的第一列由第二行，第二列相等由第二行，第二列相等,以及用前面的计算结果我们有：

8、以及用前面的计算结果我们有：矩阵矩阵R*的第二行的第二行矩阵矩阵L的第二列的第二列即有以下公式：即有以下公式：从而我们可以推出一般的计算公式：从而我们可以推出一般的计算公式：（2）Cholesky分解算法分解算法同样，利用待定系数法以及矩阵同样，利用待定系数法以及矩阵A的正定对称性，我们有：的正定对称性，我们有：我们可以推导出我们可以推导出 Cholesky三角分解得算法三角分解得算法：为保证除法运算时为保证除法运算时 ，我们由以下定理，我们由以下定理定定理理 当当A A 为为对对称称正正定定阵阵时时，A A 的的Cholesky Cholesky 分分解解必必存存在在，并并且当限定且当限定T

9、 T 的对角元素为正时，其分解是唯一的。的对角元素为正时，其分解是唯一的。有有了了矩矩阵阵三三角角三三角角分分解解后后，各各种种矩矩阵阵的的求求解解就就十十分分方方便便了了。例如：求解线性方程组例如：求解线性方程组对对A作作LR分解，有分解，有，则解方程变换为解，则解方程变换为解5.2.3矩阵的正交变换矩阵的正交变换我们从另一个角度来考虑我们从另一个角度来考虑LR分解，由前面的结论我们有分解，由前面的结论我们有此此表表达达式式可可以以了了解解为为对对A线线性性变变换换后后变变成成了了三三角角阵阵R，其其中中为为变变换换阵阵。问问题题是是我我们们能能否否用用更更为为简简单单的的一一系系列列变变换

10、换将将A变变为为上上三三角阵。角阵。（1）矩阵的正交矩阵的正交三角分解三角分解矩阵的正交分解可以写成以下形式矩阵的正交分解可以写成以下形式其中其中 Q是正交矩阵，是正交矩阵，即即 ，R是上三角矩阵，从而我们有是上三角矩阵，从而我们有(5.2.6)这这种种变变换换在在矩矩阵阵的的运运算算中中是是非非常常重重要要的的。以以下下我我们们将将分分解解为一系列较为简单的正交变换。为一系列较为简单的正交变换。（2）Householder变换变换为产生尽量简单的正交变换，我们考虑以下形式的正交方阵为产生尽量简单的正交变换，我们考虑以下形式的正交方阵(5.2.7)这这里里In是是单单位位矩矩阵阵，u为为n维维

11、向向量量，为为正正实实数数。具具有有这这种种形形式式的的正正交交变变换换称称为为H型型变变换换，我我们们可可以以通通过过以以下下步步骤骤将将矩矩阵阵A变变换换为上三角阵为上三角阵R，先用，先用H型变换将型变换将A的第一列变量变为：的第一列变量变为：再用再用 H型变换将型变换将A的第二列变为的第二列变为：第第 i步有步有为实现这一过程，我们先考虑以下简单问题。设为实现这一过程，我们先考虑以下简单问题。设我们要求一个我们要求一个H 型正交矩阵型正交矩阵Hi，使得后，使得后 n-I个元素为个元素为 0，其中其中 为常数，为常数，为使后为使后n-i个元素为个元素为0，可以取，可以取这里这里称此称此为由

12、向量为由向量定义的定义的Householder变换，并有性质。变换，并有性质。1）2）Hi是正交的是正交的3）（3）Gives变换变换Gives变换具有以下性质：变换具有以下性质：第第j列列第第i列列Gives变换具有以下性质：变换具有以下性质：1）是正交矩阵是正交矩阵2）用用 左左 乘，结果只改变乘，结果只改变 A的第的第 i行第行第 j行元素。行元素。用用 左左 乘，结果只改变乘，结果只改变 A的第的第 i行第行第 j行元素。行元素。3）对向量对向量这里这里5.2.4矩阵的谱分解矩阵的谱分解前前面面的的方方法法是是用用正正交交变变换换方方法法将将矩矩阵阵A变变为为三三角角阵阵，以以下下我我

13、们们用同样的方法将用同样的方法将A变换为对角矩阵。变换为对角矩阵。（1）对称矩阵的谱分解）对称矩阵的谱分解设设是是n阶方阵，以下分解式称为阶方阵，以下分解式称为A的谱分解式，或称为的谱分解式，或称为特征值分解式：特征值分解式：（5.2.8)其中其中 U为为 n阶正交方阵，阶正交方阵，D为对角阵，为对角阵，称称为矩阵为矩阵 A的谱或称特征值，若记的谱或称特征值，若记记记，则上式可以写成，则上式可以写成(5.2.9)如果如果A是实对称矩阵，则是实对称矩阵，则A的谱分解一定存在。的谱分解一定存在。（2）矩阵谱分解的计算方法）矩阵谱分解的计算方法（5.2.9)可以改写如下：可以改写如下：（5.2.10

14、)即即 A是是 经经 过过 正正 交交 变变 换换 后后 化化 为为 对对 角角 阵阵 的的，我我 们们 可可 以以 利利 用用Householder和和Givens方方法法的的思思路路来来构构造造这这样样的的正正交交变变换换，具具体体来来讲讲，我我们们可可以以将将（）式式中中的的U分分解解为为一一系系列列简简单单的的正正交交矩矩阵阵乘积的形式，具体算法为：乘积的形式，具体算法为：即即以下介绍如何适当选取以下介绍如何适当选取，使，使在在k充分大时接近于充分大时接近于一个对角阵。一个对角阵。Jacobi旋转法旋转法在在Gives矩阵中取矩阵中取其其中中待待定定，对对于于任任意意，可可以以验验证证

15、是是一一个个正正交交变变换换，与与解解析析几几何何中中的的旋旋转转变变换换类类似似，在在n维维空空间间中中，若若对对其其中中的的二二维维(i,j)作旋转变换，称其为作旋转变换，称其为(i,j)平面上的旋转矩阵。平面上的旋转矩阵。可以证明可以证明 Jacobi算法必有算法必有 为对角矩阵。为对角矩阵。在（）如果取在（）如果取为为的的正交正交三角变换，则三角变换，则为著名的求特征值的为著名的求特征值的QR算法算法3）QR算法算法（a）（b）对）对做做的正交分解的正交分解取取A为任意方阵，可以证明为任意方阵，可以证明 “基本收敛基本收敛”于一个上三角矩阵，于一个上三角矩阵，而对角线元素为其特征值。而

16、对角线元素为其特征值。5.2.5矩阵的奇异值分解矩阵的奇异值分解设设是任意非零矩阵，则是任意非零矩阵，则为为A的奇异值分解，其中的奇异值分解，其中U和和V分别为分别为n和和m阶正交方阵，阶正交方阵，D为为nm阶对角阵，其非对角元素均为阶对角阵，其非对角元素均为0，D的对角线元素的对角线元素称为矩阵称为矩阵A的奇异值。奇异值的分解与矩阵的谱分解方法类的奇异值。奇异值的分解与矩阵的谱分解方法类似。似。5.2.6矩阵的广义逆矩阵的广义逆矩阵的广义逆在统计计算中具有重要的作用，它是由矩阵的逆矩阵的广义逆在统计计算中具有重要的作用，它是由矩阵的逆的概念进一步一般化而来。的概念进一步一般化而来。设设为为n

17、m阶矩阵，阶矩阵，G为为nm矩阵，如果矩阵，如果G满足：满足：（1）AGA=A（2）GAG=G（3）（）（AG）=AG（即（即AG为对称矩阵）为对称矩阵）（4）（）（GA）=GA（即（即GA为对称矩阵）为对称矩阵）满足这四个条件的某几个或全部，则称矩阵满足这四个条件的某几个或全部，则称矩阵G为矩阵为矩阵A的广义逆。的广义逆。定义：定义：1）满足上面第一条的矩阵）满足上面第一条的矩阵G称为称为A的减号逆，记为的减号逆，记为G=A-2）满足上面条件（）满足上面条件（1），（），（2）的矩阵）的矩阵G称为称为A的自反广义的自反广义逆记为逆记为3）满足上面所有的四个条件称）满足上面所有的四个条件称G为

18、为A的加号逆极为的加号逆极为5.2.7线性方程组的解线性方程组的解我们可以用矩阵的广义逆这一有力的工具来解线性方程组我们可以用矩阵的广义逆这一有力的工具来解线性方程组其中其中A 为为nm矩阵，矩阵，b 为为n1的向量，的向量，x 为为m1的位知的位知向量，若（向量，若（5.2.11）有解，则称其为相容方程。）有解，则称其为相容方程。（5.2.11)如果存在如果存在使得使得则称则称为方程（为方程（5.2.11)的最小二乘解，其中的最小二乘解，其中表示向量表示向量y的模，其定义为。的模，其定义为。以下不证明，给出相容方程的一般表达式以下不证明，给出相容方程的一般表达式或或u任意任意5.2.8矩阵的

19、范数（模）矩阵的范数（模）矩阵的范数与条件数是矩阵代数运算的重要概念之一，范数为矩阵的范数与条件数是矩阵代数运算的重要概念之一，范数为任意矩阵定义了一个函数，而条件数是将来进行计算时对计算任意矩阵定义了一个函数，而条件数是将来进行计算时对计算精度的一种衡量。请见范数的定义：精度的一种衡量。请见范数的定义：定义：设定义：设Rm为为m维实数空间，维实数空间，是是Rm到实数轴到实数轴R1的一的一个映射，若个映射，若满足满足（1）时成立时成立（2）对任意常数）对任意常数，（3）则称则称为向量为向量x的范数的范数常用的范数有：常用的范数有：1）2）3）矩阵范数的定义：矩阵范数的定义：定义：设定义：设是是

20、nm维实数空间，维实数空间，是是到实数到实数轴的一个映射，如果满足：轴的一个映射，如果满足：（2）对任意常数）对任意常数，有有（3）则称则称为矩阵为矩阵A的范数的范数，特别，如果，特别，如果（1）（4）对）对则称则称为矩阵为矩阵A的相容范数。的相容范数。矩阵矩阵A的常用范数的常用范数对于任意矩阵对于任意矩阵定义：定义：容易证明容易证明是矩阵是矩阵A的相容范数。则常用的矩阵范数：的相容范数。则常用的矩阵范数：（1）（2）（3）这里我们主要讨论这里我们主要讨论，并记以，并记以定理：设定理：设A是是nm矩阵，则有矩阵，则有这里这里为矩阵为矩阵A的绝对值最大的特征值（奇异值），的绝对值最大的特征值（奇

21、异值），从而我们给矩阵定义了一个范数。从而我们给矩阵定义了一个范数。5.2.9矩阵的条件数矩阵的条件数定义：称定义：称为为的条件数，记为：的条件数，记为：条件数有以下性质：条件数有以下性质：（1）即矩阵的条件数为矩阵即矩阵的条件数为矩阵的最大奇异值比最小奇异值的绝对值。特别当的最大奇异值比最小奇异值的绝对值。特别当A为对称阵时为对称阵时即为即为A的最大特征值和最小特征值之比的绝对值。的最大特征值和最小特征值之比的绝对值。（2），当当A为正交阵时等式成立。为正交阵时等式成立。（3）对任意非零常数有，）对任意非零常数有，（4）若）若，则，则（5）矩阵条件数的重要意义在于，可以判别一个求逆矩阵的病矩

22、阵条件数的重要意义在于，可以判别一个求逆矩阵的病态性。当系数矩阵的条件数非常大时，即存在接近与态性。当系数矩阵的条件数非常大时，即存在接近与0的特征的特征值，将导致解的误差急剧放大。或者说得出的解不可信。值，将导致解的误差急剧放大。或者说得出的解不可信。对于最小二乘估计：对于最小二乘估计：我们最关心的是对称矩阵我们最关心的是对称矩阵的条件数是否正常，从而判的条件数是否正常，从而判别最小二乘估计是否可信。别最小二乘估计是否可信。5.3IMSL库中的矩阵计算模块库中的矩阵计算模块nIMSL库中的库中的maths.lib库中有非常丰富的矩阵库中有非常丰富的矩阵计算模块包括：计算模块包括：n矩阵乘、矩

23、阵求逆、广义逆矩阵与向量乘等。矩阵乘、矩阵求逆、广义逆矩阵与向量乘等。n矩阵的谱分解、矩阵的奇异值分解。矩阵的谱分解、矩阵的奇异值分解。n矩阵、向量的模计算。矩阵、向量的模计算。n矩阵的条件数的计算。矩阵的条件数的计算。n求解线性方程组。求解线性方程组。打开打开maths.lib我们可以看到我们可以看到基本矩阵向量运算模块基本矩阵向量运算模块求解线性方程组求解线性方程组矩阵的特征值分析矩阵的特征值分析5.3.1矩阵的基本运算模块矩阵的基本运算模块基本线性代数计算基本线性代数计算矩阵的变换矩阵的变换矩阵相乘矩阵相乘矩阵与向量乘矩阵与向量乘IIntegerSRealCComplexDDoubleZ

24、Double complexSDSingle and DoubleCZ Single and double complexDQDouble and QuadrupleZQDouble and quadruple complex基本代数运算子程序名的第一个字符的含义：基本代数运算子程序名的第一个字符的含义：例如：例如：SGEMM Matrix-matrix multiply,general表示单精度的矩阵与矩阵的乘法表示单精度的矩阵与矩阵的乘法 计算一个矩阵计算一个矩阵A A 与其转置阵的乘积，这里与其转置阵的乘积，这里打开打开注意注意B是输出变量是输出变量编制程序：在工程中建立数据文件，在程序

25、中打开数据文件编制程序：在工程中建立数据文件，在程序中打开数据文件和结果文件，再调用矩阵乘子程序和矩阵打印程序。和结果文件，再调用矩阵乘子程序和矩阵打印程序。5.3.2解线性方程组系统解线性方程组系统点击线性系统，可看到列表点击线性系统，可看到列表解线性系统的子程序名的意义：解线性系统的子程序名的意义：解方法示意图解方法示意图程序名的意义：程序名的意义：LSARG高精度求解高精度求解一般矩形矩阵一般矩形矩阵LFCRG三角分解与条件数三角分解与条件数一般矩形矩阵一般矩形矩阵LSACG矩形复数阵矩形复数阵高精度求解高精度求解LSVRR实系数矩阵实系数矩阵奇异值分解奇异值分解例：例：对系数矩阵为正定

26、对称阵的线性方程组对系数矩阵为正定对称阵的线性方程组先对系数矩阵进行先对系数矩阵进行Cholesky分解，再求解线性方程组分解，再求解线性方程组解：我们使用解：我们使用LCHRG子程序子程序功能：对对称正定矩阵、对称半正定矩阵进行功能：对对称正定矩阵、对称半正定矩阵进行Cholesky分解。分解。现在我们来看现在我们来看LCHRG子程序。子程序。UsageCALL LCHRG(N,A,LDA,PIVOT,IPVT,FAC,LDFAC)编程，先将对称正定矩阵建立一个数据文件编程，先将对称正定矩阵建立一个数据文件CHE.DAT,然后编程如然后编程如下：下：计算结果为：计算结果为：5.3.3求矩阵的

27、谱分解求矩阵的谱分解求求矩矩阵阵的的谱谱分分解解，求求矩矩阵阵的的特特征征值值和和特特征征向向量量。常常规规线线性性特特征值系统求解问题是征值系统求解问题是这里为一个矩阵，而为特征值，为对应于的一个向量。这里为一个矩阵，而为特征值，为对应于的一个向量。广义线性特征值系统求解问题是广义线性特征值系统求解问题是常常规规线线性性特特征征值值系系统统求求解解程程序序以以E开开头头，而而广广义义线线性性特特征征值值系系统统求求解以解以GE开头，具体的程序分类见下列表格。开头，具体的程序分类见下列表格。计算以下矩阵的全部特征值和特征向量计算以下矩阵的全部特征值和特征向量nUsagenCALLEVCRG(N

28、,A,LDA,EVAL,EVEC,LDEVEC)nArgumentsnNOrderofthematrix.(Input)nAFloating-pointarraycontainingthematrix.(Input)nLDALeadingdimensionofAexactlyasspecifiedinthedimensionstatementofthecallingprogram.(Input)nEVALComplexarrayofsizeNcontainingtheeigenvaluesofAindecreasingorderofmagnitude.(Output)nEVECComplexa

29、rraycontainingthematrixofeigenvectors.(Output)TheJ-theigenvector,correspondingtoEVAL(J),isstoredintheJ-thcolumn.EachvectorisnormalizedtohaveEuclideanlengthequaltothevalueone.nLDEVECLeadingdimensionofEVECexactlyasspecifiedinthedimensionstatementofthecallingprogram.(Input)n计算结果如下：计算结果如下：nEVALn123n(2.0

30、00,4.000)(2.000,-4.000)(1.000,.000)nEVECn123n1(.3162,.3162)(.3162,-.3162)(.4082,.0000)n2(.0000,.6325)(.0000,-.6325)(.8165,.0000)n3(.6325,.0000)(.6325,.0000)(.4082,.0000)nPerformanceindex=.015n从而我们得到三个复数特征值，及其对应的三个特征向量。从而我们得到三个复数特征值，及其对应的三个特征向量。习题习题（1）编程求解矩阵编程求解矩阵A与其转置矩阵与其转置矩阵A的乘积的乘积AA（2）编程求下面对称阵的三角分解和编程求下面对称阵的三角分解和CHOLESKY分解分解（3）编程求解下面的方程组编程求解下面的方程组（4）求矩阵的特征值与特征向量求矩阵的特征值与特征向量