第五章_矩阵代数数值计算.ppt

第五章矩阵代数数值计算 一 矩阵的基本运算二 矩阵的三角分解三 矩阵的正交变换四 矩阵的谱分解五 IMSL中的线性系统 特征值分析模块 矩阵代数运算是统计模型的基础 统计模型的所有估计几乎都是用矩阵代数运算计算出结果 例如最小二乘估计 典型相关分析 因子分析以及各类回归分析 从计算的角度来说为使计算结果可靠 我们总是先对矩阵进行三角分解 然后进行各种计算例如 矩阵的逆 求解线性方程组以及对矩阵进行谱分解等 本章首先介绍矩阵的三角分解 然后引导学习者使用IMSL和SASD中的丰富矩阵的算法 将它们拼接起来就可以解决各种矩阵的计算 5 1引言 矩阵代数运算在数值计算中起着基础性的作用 无论我们建立了多么复杂的数学模型 最终我们总是要把它变为矩阵代数的形式 特别是统计模型 无论是多元线性回归 广义线性模型 多元统计分析无不与矩阵代数有着密切的联系 我们所研究的对象 即样本可以看成是一个矩阵 而统计上的协方差 实际上是该矩阵的转置与该矩阵相乘形成的新矩阵的元素 而回归的最小二乘估计的算法为 5 1 1 包含了矩阵的乘 矩阵的求逆及矩阵与向量的乘法等等 而特征值与特征向量在数理统计理都有明确的条件统计含义 因此我们将在这一章介绍矩阵运算的基本数值方法 以及如何调用IMSL和SASD中丰富的算法模块 5 2矩阵数值计算基础 对于一般的二维样本 我们总可以写成如下的矩阵形式 5 2 1 从计算的角度处理矩阵问题的一个最有效的方法是 将一个一般的矩阵分解为几个简单矩阵的乘积形式 其中最便于计算的是三角矩阵 以上三角矩阵为例 由其性质 两个上三角阵的和 积仍为上三角阵 上三角阵的特征值就是其对角线元素 系数矩阵为上三角阵的线性线性方程组是最容易求解的 上三角阵的逆阵仍然是上三角阵 因此处理矩阵计算问题的关键是将一般矩阵化为三角阵和对角阵的形式 然后进行计算 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分解是统计计算中最常用的分解方法之一 因为我们的协方差矩阵 相关矩阵都是使用这种分解方法 5 2 2矩阵的三角分解算法以上四种分解是类似的 使用待定系数法 1 以LR 分解为例 设 其中A为正定阵 并记分解已为以下形式 利用矩阵相对应元素相等的事实 我们立即得到 现在我们可以计算矩阵L的第一列 由第二行 第二列相等 以及用前面的计算结果我们有 矩阵R 的第二行 矩阵L的第二列 即有以下公式 从而我们可以推出一般的计算公式 2 Cholesky分解算法同样 利用待定系数法以及矩阵A的正定对称性 我们有 我们可以推导出Cholesky三角分解得算法 为保证除法运算时 我们由以下定理定理5 2 1当A为对称正定阵时 A的Cholesky分解必存在 并且当限定T的对角元素为正时 其分解是唯一的 有了矩阵三角 三角分解后 各种矩阵的求解就十分方便了 例如 求解线性方程组 对A作LR分解 有 则解方程变换为解 5 2 3矩阵的正交变换我们从另一个角度来考虑LR分解 由前面的结论我们有 此表达式可以了解为对A线性变换后变成了三角阵R 其中为变换阵 问题是我们能否用更为简单的一系列变换将A变为上三角阵 1 矩阵的正交 三角分解矩阵的正交分解可以写成以下形式 其中Q是正交矩阵 即 R是上三角矩阵 从而我们有 5 2 6 这种变换在矩阵的运算中是非常重要的 以下我们将分解为一系列较为简单的正交变换 2 Householder变换为产生尽量简单的正交变换 我们考虑以下形式的正交方阵 5 2 7 这里In是单位矩阵 u为n维向量 为正实数 具有这种形式的正交变换称为H型变换 我们可以通过以下步骤将矩阵A变换为上三角阵R 先用H型变换将A的第一列变量变为 再用H型变换将A的第二列变为 第i步有 为实现这一过程 我们先考虑以下简单问题 设 我们要求一个H型正交矩阵Hi 使得后n I个元素为0 其中为常数 为使后n i个元素为0 可以取 这里 称此为由向量定义的Householder变换 并有性质 1 2 Hi是正交的 3 3 Gives变换 Gives变换具有以下性质 第j列 第i列 Gives变换具有以下性质 1 是正交矩阵 2 用左乘 结果只改变A的第i行第j行元素 用左乘 结果只改变A的第i行第j行元素 3 对向量 这里 5 2 4矩阵的谱分解前面的方法是用正交变换方法将矩阵A变为三角阵 以下我们用同样的方法将A变换为对角矩阵 1 对称矩阵的谱分解 设是n阶方阵 以下分解式称为A的谱分解式 或称为特征值分解式 5 2 8 其中U为n阶正交方阵 D为对角阵 称为矩阵A的谱或称特征值 若记 记 则上式可以写成 5 2 9 如果A是实对称矩阵 则A的谱分解一定存在 2 矩阵谱分解的计算方法 5 2 9 可以改写如下 5 2 10 即A是经过正交变换后化为对角阵的 我们可以利用Householder和Givens方法的思路来构造这样的正交变换 具体来讲 我们可以将 5 2 8 式中的U分解为一系列简单的正交矩阵乘积的形式 具体算法为 即 以下介绍如何适当选取 使在k充分大时接近于一个对角阵 Jacobi旋转法 在Gives矩阵中取 其中待定 对于任意 可以验证是一个正交变换 与解析几何中的旋转变换类似 在n维空间中 若对其中的二维 i j 作旋转变换 称其为 i j 平面上的旋转矩阵 可以证明Jacobi算法必有为对角矩阵 在 5 2 9 如果取为的正交 三角变换 则 为著名的求特征值的QR算法 3 QR算法 a b 对做的正交分解 取 A为任意方阵 可以证明 基本收敛 于一个上三角矩阵 而对角线元素为其特征值 5 2 5矩阵的奇异值分解 设是任意非零矩阵 则 为A的奇异值分解 其中U和V分别为n和m阶正交方阵 D为n m阶对角阵 其非对角元素均为0 D的对角线元素称为矩阵A的奇异值 奇异值的分解与矩阵的谱分解方法类似 5 2 6矩阵的广义逆 矩阵的广义逆在统计计算中具有重要的作用 它是由矩阵的逆的概念进一步一般化而来 设为n m阶矩阵 G为n m矩阵 如果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为A的加号逆极为 5 2 7线性方程组的解我们可以用矩阵的广义逆这一有力的工具来解线性方程组 其中A为n m矩阵 b为n 1的向量 x为m 1的位知向量 若 5 2 11 有解 则称其为相容方程 5 2 11 如果存在使得 则称为方程 5 2 11 的最小二乘解 其中表示向量y的模 其定义为 以下不证明 给出相容方程的一般表达式 或 u任意 5 2 8矩阵的范数 模 矩阵的范数与条件数是矩阵代数运算的重要概念之一 范数为任意矩阵定义了一个函数 而条件数是将来进行计算时对计算精度的一种衡量 请见范数的定义 定义 设Rm为m维实数空间 是Rm到实数轴R1的一个映射 若满足 1 时成立 2 对任意常数 3 则称为向量x的范数 常用的范数有 1 2 3 矩阵范数的定义 定义 设是n m维实数空间 是到实数轴的一个映射 如果满足 2 对任意常数 有 3 则称为矩阵A的范数 特别 如果 1 4 对 则称为矩阵A的相容范数 矩阵A的常用范数对于任意矩阵定义 容易证明是矩阵A的相容范数 则常用的矩阵范数 1 2 3 这里我们主要讨论 并记以 定理 设A是n m矩阵 则有 这里为矩阵A的绝对值最大的特征值 奇异值 从而我们给矩阵定义了一个范数 5 2 9矩阵的条件数定义 称为的条件数 记为 条件数有以下性质 1 即矩阵的条件数为矩阵的最大奇异值比最小奇异值的绝对值 特别当A为对称阵时即为A的最大特征值和最小特征值之比的绝对值 2 当A为正交阵时等式成立 3 对任意非零常数有 4 若 则 5 矩阵条件数的重要意义在于 可以判别一个求逆矩阵的病态性 当系数矩阵的条件数非常大时 即存在接近与0的特征值 将导致解的误差急剧放大 或者说得出的解不可信 对于最小二乘估计 我们最关心的是对称矩阵的条件数是否正常 从而判别最小二乘估计是否可信 5 3IMSL库中的矩阵计算模块 IMSL库中的maths lib库中有非常丰富的矩阵计算模块包括 矩阵乘 矩阵求逆 广义逆矩阵与向量乘等 矩阵的谱分解 矩阵的奇异值分解 矩阵 向量的模计算 矩阵的条件数的计算 求解线性方程组 打开maths lib我们可以看到 基本矩阵向量运算模块 求解线性方程组 矩阵的特征值分析 5 3 1矩阵的基本运算模块 基本线性代数计算 矩阵的变换 矩阵相乘 矩阵与向量乘 IIntegerSRealCComplexDDoubleZDoublecomplexSDSingleandDoubleCZSingleanddoublecomplexDQDoubleandQuadrupleZQDoubleandquadruplecomplex 基本代数运算子程序名的第一个字符的含义 例如 SGEMMMatrix matrixmultiply general 表示单精度的矩阵与矩阵的乘法 例5 3 1计算一个矩阵A与其转置阵的乘积 这里 打开 注意B是输出变量 编制程序 在工程中建立数据文件 在程序中打开数据文件和结果文件 再调用矩阵乘子程序和矩阵打印程序 5 3 2解线性方程组系统 点击线性系统 可看到列表 解线性系统的子程序名的意义 解方法示意图 程序名的意义 LSARG 高精度求解 一般矩形矩阵 LFCRG 三角分解与条件数 一般矩形矩阵 LSACG 矩形复数阵 高精度求解 LSVRR 实系数矩阵 奇异值分解 例5 3 2 对系数矩阵为正定对称阵的线性方程组 先对系数矩阵进行Cholesky分解 再求解线性方程组解 我们使用LCHRG子程序功能 对对称正定矩阵 对称半正定矩阵进行Cholesky分解 现在我们来看LCHRG子程序 UsageCALLLCHRG N A LDA PIVOT IPVT FAC LDFAC 编程 先将对称正定矩阵建立一个数据文件CHE DAT 然后编程如下 计算结果为 5 3 3求矩阵的谱分解 求矩阵的谱分解 求矩阵的特征值和特征向量 常规线性特征值系统求解问题是这里为一个矩阵 而为特征值 为对应于的一个向量 广义线性特征值系统求解问题是 常规线性特征值系统求解程序以E开头 而广义线性特征值系统求解以GE开头 具体的程序分类见下列表格 例5 3 3计算以下矩阵的全部特征值和特征向量 UsageCALLEVCRG N A LDA EVAL EVEC LDEVEC ArgumentsN Orderofthematrix Input A Floating pointarraycontainingthematrix Input LDA LeadingdimensionofAexactlyasspecifiedinthedimensionstatementofthecallingprogram Input EVAL ComplexarrayofsizeNcontainingtheeigenvaluesofAindecreasingorderofmagnitude Output EVEC Complexarraycontainingthematrixofeigenvectors Output TheJ theigenvector correspondingtoEVAL J isstoredintheJ thcolumn EachvectorisnormalizedtohaveEuclideanlengthequaltothevalueone LDEVEC LeadingdimensionofEVECexactlyasspecifiedinthedimensionstatementofthecallingprogram Input 计算结果如下 EV