多元回归中的几个常用的矩阵式

合集下载

多元线性回归实例分析

多元线性回归实例分析

SPSS--回归-多元线性回归模型案例解析!(一)多元线性回归,主要是研究一个因变量与多个自变量之间的相关关系,跟一元回归原理差不多,区别在于影响因素(自变量)更多些而已,例如:一元线性回归方程为:毫无疑问,多元线性回归方程应该为:上图中的 x1, x2, xp分别代表“自变量”Xp截止,代表有P个自变量,如果有“N组样本,那么这个多元线性回归,将会组成一个矩阵,如下图所示:那么,多元线性回归方程矩阵形式为:其中:代表随机误差,其中随机误差分为:可解释的误差和不可解释的误差,随机误差必须满足以下四个条件,多元线性方程才有意义(一元线性方程也一样)1:服成正太分布,即指:随机误差必须是服成正太分别的随机变量。

2:无偏性假设,即指:期望值为03:同共方差性假设,即指,所有的随机误差变量方差都相等4:独立性假设,即指:所有的随机误差变量都相互独立,可以用协方差解释。

今天跟大家一起讨论一下,SPSS---多元线性回归的具体操作过程,下面以教程教程数据为例,分析汽车特征与汽车销售量之间的关系。

通过分析汽车特征跟汽车销售量的关系,建立拟合多元线性回归模型。

数据如下图所示:点击“分析”——回归——线性——进入如下图所示的界面:将“销售量”作为“因变量”拖入因变量框内,将“车长,车宽,耗油率,车净重等10个自变量拖入自变量框内,如上图所示,在“方法”旁边,选择“逐步”,当然,你也可以选择其它的方式,如果你选择“进入”默认的方式,在分析结果中,将会得到如下图所示的结果:(所有的自变量,都会强行进入)如果你选择“逐步”这个方法,将会得到如下图所示的结果:(将会根据预先设定的“F统计量的概率值进行筛选,最先进入回归方程的“自变量”应该是跟“因变量”关系最为密切,贡献最大的,如下图可以看出,车的价格和车轴跟因变量关系最为密切,符合判断条件的概率值必须小于0.05,当概率值大于等于0.1时将会被剔除)“选择变量(E)" 框内,我并没有输入数据,如果你需要对某个“自变量”进行条件筛选,可以将那个自变量,移入“选择变量框”内,有一个前提就是:该变量从未在另一个目标列表中出现!,再点击“规则”设定相应的“筛选条件”即可,如下图所示:点击“统计量”弹出如下所示的框,如下所示:在“回归系数”下面勾选“估计,在右侧勾选”模型拟合度“ 和”共线性诊断“ 两个选项,再勾选“个案诊断”再点击“离群值”一般默认值为“3”,(设定异常值的依据,只有当残差超过3倍标准差的观测才会被当做异常值)点击继续。

多元线性回归

多元线性回归

36
目录 上页 下页 返回 结束
§5.4 回归方程的显著性检验
2019/11/5
中国人民大学六西格玛质量管理研究中心
37
目录 上页 下页 返回 结束
§5.4 回归方程的显著性检验
2019/11/5
中国人民大学六西格玛质量管理研究中心
38
目录 上页 下页 返回 结束
§5.4 回归方程的显著性检验
2019/11/5
中国人民大学六西格玛质量管理研究中心
16
目录 上页 下页 返回 结束
§5.2 多元回归参数的估计
2019/11/5
中国人民大学六西格玛质量管理研究中心
17
目录 上页 下页 返回 结束
§5.2 多元回归参数的估计
2019/11/5
中国人民大学六西格玛质量管理研究中心
18
目录 上页 下页 返回 结束
§5.4 回归方程的显著性检验
在一元线性回归中,回归系数显著性的t检验与回归方 程显著性的F检验是等价的,而在多元线性回归中,这 两种检验不同。
2019/11/5
中国人民大学六西格玛质量管理研究中心
43
目录 上页 下页 返回 结束
§5.4 回归方程的显著性检验
2019/11/5
中国人民大学六西格玛质量管理研究中心
27
目录 上页 下页 返回 结束
§5.3 参数估计量的性质
2019/11/5
中国人民大学六西格玛质量管理研究中心
28
目录 上页 下页 返回 结束
§5.3 参数估计量的性质
性质4 Gauss-Markov定理
2019/11/5
中国人民大学六西格玛质量管理研究中心
29

多元回归模型与回归方程

多元回归模型与回归方程
1. 用样本统计量 ˆ0 , ˆ1 , ˆ2 , , ˆ p 估计回归方 程中的 参数 0 , 1 , 2 , , p 时得到的方程
2. 由最小二乘法求得 3. 一般形式为
yˆ ˆ0 ˆ1x1 ˆ2x2 ˆp xp
▪ ˆ0 , ˆ1 , ˆ2 , , ˆ p是 0 , 1 , 2 , , p
p-1找出临界值F 4. 作出决策:若F>F ,拒绝H0
回归系数检验和推断
回归系数的检验
1. 线性关系检验通过后,对各个回归系数有选 择地进行一次或多次检验
2. 究竟要对哪几个回归系数进行检验,通常需 要在建立模型之前作出决定
3. 对回归系数检验的个数进行限制,以避免犯 过多的第一类错误(弃真错误)
估计值
▪ yˆ 是 y 的估计值
参数的最小二乘估计
普通最小二乘估计
普通最小二乘估计
随机抽取被解释变量和解释变量的 n 组样本观测值:
(Yi , X ji ), i 1,2, , n, j 0,1,2, k
如果模型的参数估计值已经得到,则有:
Yˆi ˆ0 ˆ1 X 1i ˆ2 X 2i ˆki X Ki
rPerformance,Runtime = -0.98841
= Error = Collinearity
标量符号 2、随机误差项具有零均值、同方差及不序列相关
E(i ) 0
i 1,2, , n
Var ( i
)
E
(
2 i
)
2
i 1,2, , n
Cov(i , j ) E(i j ) 0
1. 回归平方和占总平方和的比例 2. 计算公式为
3. 因变量取值的变差中,能被估计的多元回 归方程所解释的比例

多元的线性回归

多元的线性回归

多元线性回归模型一、多元线性回归模型的一般形式设随机变量y 与一般变量p x x x ,,,21 的线性回归模型为:εββββ+++++=p p x x x y 22110写成矩阵形式为:εβ+=X y 其中:⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=n y y y y 21 ⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=np n n p p x x x x x x x x x X 212222********* ⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=p ββββ 10 ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=n εεεε 21二、多元线性回归模型的基本假定1、解释变量p x x x ,,,21 是确定性变量,不是随机变量,且要求n p X rank <+=1)(。

这里的n p X rank <+=1)(表明设计矩阵X 中自变量列之间不相关,样本容量的个数应大于解释变量的个数,X 是一满秩矩阵。

2、随机误差项具有0均值和等方差,即:⎪⎩⎪⎨⎧⎩⎨⎧=≠====),,2,1,(,,0,),cov(,,2,1,0)(2n j i j i j i n i E j i i σεεε 0)(=i E ε,即假设观测值没有系统误差,随机误差i ε的平均值为0,随机误差iε的协方差为0表明随机误差项在不同的样本点之间是不相关的(在正态假定下即为独立),不存在序列相关,并且具有相同的精度。

3、正态分布的假定条件为:⎩⎨⎧=相互独立n i ni N εεεσε ,,,,2,1),,0(~212,矩阵表示:),0(~2n I N σε,由该假定和多元正态分布的性质可知,随机变量y 服从n 维正态分布,回归模型的期望向量为:βX y E =)(;n I y 2)var(σ= 因此有),(~2n I X N y σβ 三、多元线性回归方程的解释对于一般情况含有p 个自变量的回归方程p p x x x y E ββββ++++= 22110)(的解释,每个回归系数i β表示在回归方程中其他自变量保持不变的情况下,自变量i x 每增加一个单位时因变量y 的平均增加程度。

多元线性回归模型的矩阵表示课件

多元线性回归模型的矩阵表示课件
根据上述公式计算决定系数,需要先根据回归
直线计算 Yi的理论值,然后计算回归残差序列,
再结合样本数据进行计算。
25
第四节 统计推断和预测
一、参数估计量的标准化 二、统计推断和检验 三、预测
26
一、参数估计量的标准化
在满足模型假设的情况下,多元线性回归模型 参数的最小二乘估计量是线性无偏估计。
Y1 0 1 X 11 K X K1 1
Yn 0 1 X 1n K X K n
Y1
Y
Yn
X i1
X i
X i n
1
l
1
0
K
1
n
1 X11 X K1
X l, X1,, X K
1 X1n X Kn
Y 0 1 X 1 2 X 2 K X K X
S.E. of regression 0.007246 Akaike info criterion -6.849241
Sum squared resid 0.000683 Schwarz criterion -6.704381
Log likelihood 57.79393 F-statistic
(1)、变量Y和X1,X K 之间存在多元线性随
机函数关系 Y 0 1X1 K X K ;
(2)、Ei 0 对任意 i 都成立;
(3)、Vari 2 ,与 i 无关;
(4)、误差项不相关,当 i j 时,E i j 0
(5)、解释变量都是确定性的而非随机变量, 且解释变量之间不存在线性关系;
bk k
seˆ(bk )
= bk
seˆ(bk )
t / 2(n-K-1)
如果t 统计量数值不满足上述不等式,意味着 可以拒绝原假设,不能认为第k个解释变量是 不重要的,称模型的第k个解释变量通过了显

第六章_多元回归分析的矩阵运算

第六章_多元回归分析的矩阵运算

第六章_多元回归分析的矩阵运算多元回归分析是统计学中重要的分析方法之一,用于研究多个自变量对一个因变量的影响关系。

在进行多元回归分析时,矩阵运算是一个重要的工具,可以帮助我们简化计算过程,提高效率。

本文将介绍多元回归分析中的矩阵运算。

多元回归模型可以表示为:Y=Xβ+ε其中,Y是因变量的观测值向量,X是自变量的观测值矩阵,β是自变量的系数向量,ε是误差项的观测值向量。

我们将自变量的观测值矩阵X进行标准化处理,使得X的每一列均值为0,标准差为1,即mean(X) = 0,std(X) = 1、这样做的目的是消除自变量之间的量纲差异,方便进行比较。

在进行多元回归分析时,我们需要使用最小二乘法来估计模型的参数β。

最小二乘法的估计公式为:β=(X'X)^(-1)X'Y其中,X'表示X的转置,^(-1)表示矩阵的逆运算。

矩阵的转置运算可以通过将矩阵的行转换为列,列转换为行来实现。

例如,矩阵X的转置X'的第i行第j列元素等于X的第j行第i列元素,可表示为X'ij = Xji。

矩阵的逆运算表示将矩阵转换为与其相乘后得到单位矩阵的矩阵。

例如,矩阵A的逆矩阵A^(-1)满足A^(-1)*A=I,其中I为单位矩阵。

在进行最小二乘法估计时,我们需要计算矩阵X'X的逆矩阵。

若X'X为可逆矩阵,则矩阵X'X的逆矩阵可以写为(X'X)^(-1) = 1/,X'X, *adj(X'X),其中,X'X,表示矩阵X'X的行列式,adj(X'X)为X'X的伴随矩阵。

矩阵的行列式表示矩阵的性质,可以通过计算矩阵的特征值(即矩阵的特征多项式的根)来得到。

例如,矩阵A的行列式,A,可以通过计算A的特征值λ1,λ2,…,λn的乘积来得到,即,A,=λ1*λ2*…*λn。

矩阵的伴随矩阵可以通过矩阵的代数余子式来计算。

矩阵A的第i行第j列元素的代数余子式Aij表示在A中去掉第i行第j列后,剩余矩阵的行列式。

多元回归分析原理及例子

多元回归分析原理及例子

多元回归分析原理回归分析是一种处理变量的统计相关关系的一种数理统计方法。

回归分析的基本思想是: 虽然自变量和因变量之间没有严格的、确定性的函数关系, 但可以设法找出最能代表它们之间关系的数学表达形式。

回归分析主要解决以下几个方面的问题:(1) 确定几个特定的变量之间是否存在相关关系, 如果存在的话, 找出它们之间合适的数学表达式; (2) 根据一个或几个变量的值, 预测或控制另一个变量的取值, 并且可以知道这种预测或控制能达到什么样的精确度;(3) 进行因素分析。

例如在对于共同影响一个变量的许多变量(因素)之间, 找出哪些是重要因素, 哪些是次要因素, 这些因素之间又有什么关系等等。

回归分析有很广泛的应用, 例如实验数据的一般处理, 经验公式的求得, 因素分析, 产品质量的控制, 气象及地震预报, 自动控制中数学模型的制定等等。

多元回归分析是研究多个变量之间关系的回归分析方法, 按因变量和自变量的数量对应关系可划分为一个因变量对多个自变量的回归分析(简称为“一对多”回归分析)及多个因变量对多个自变量的回归分析(简称为“多对多”回归分析), 按回归模型类型可划分为线性回归分析和非线性回归分析。

本“多元回归分析原理”是针对均匀设计3.00软件的使用而编制的, 它不是多元回归分析的全面内容, 欲了解多元回归分析的其他内容请参阅回归分析方面的书籍。

本部分内容分七个部分, §1~§4介绍“一对多”线性回归分析, 包括数学模型、回归系数估计、回归方程及回归系数的显著性检验、逐步回归分析方法。

“一对多”线性回归分析是多元回归分析的基础, “多对多”回归分析的内容与“一对多”的相应内容类似, §5介绍“多对多”线性回归的数学模型,§6介绍“多对多”回归的双重筛选逐步回归法。

§7简要介绍非线性回归分析。

§1 一对多线性回归分析的数学模型§2 回归系数的最小二乘估计§3 回归方程及回归系数的显著性检验§4 逐步回归分析§5 多对多线性回归数学模型§6 双重筛选逐步回归§7 非线性回归模型§1 一对多线性回归分析的数学模型设随机变量与个自变量存在线性关系:, (1.1)(1.1)式称为回归方程, 式中为回归系数,为随机误差。

多元回归分析中常用的矩阵算法

多元回归分析中常用的矩阵算法

毕业设计(论文)任务书课题名称多元回归分析中常用的矩阵算法学院专业班级姓名学号毕业设计(论文)的主要内容及要求:(1)首先要了解多元回归分析中的矩阵算法的研究背景。

(2)查阅相关文献(至少4-5篇),并查阅1-2篇外文文献。

(3)熟悉相关矩阵算法;掌握多元回归分析的基本理论知识;(4)完成各种矩阵算法的程序编写,并将其运用于多元回归分析。

(5)通过实例验证算法的准确性,然后进行修改优化。

(6)整理相关资料,完成毕业论文的写作。

(7)对论文进行全面修改、完善,准备论文答辩。

指导教师签字:摘要在多元回归分析的计算中,观测数据一般用矩阵表示,对数据的分析转化为对数据矩阵的分析计算问题.如线性方程组的求解,矩阵的分解,矩阵的变换,特征值和特征向量的计算等.这些常见的矩阵计算问题也是多元回归分析中经常遇到的问题.本文主要介绍了多元回归分析中常用的矩阵分解及其算法,其中包括三角分解,正交三角-分解,正交分解. 然后针对每一种分解我们讨论了它们的一些常用算法,并在计算机上通过Matlab软件编程实现这些算法,最后再介绍了这些矩阵算法在多元回归分析中的应用.本文给出的算法是多元回归分析计算的基础,对应用多元回归分析解决实际问题具有很重要的意义.关键词:矩阵分解;矩阵变换;算法;回归分析AbstractIn the calculation of multiple regression analysis, the observed data generally represented by matrix, the analysis of datas often transform into the analysis of the matrix. Such as the solution of linear equations, matrix decomposition, matrix transform, the computation of eigenvalues and eigenvectors. These common matrix computation problems are often encountered in the multivariate regression analysis of the problem.This paper mainly introduces the commonly used matrix decomposition and its algorithm in the multiple regression analysis,including triangular decomposition, QR decomposition, orthogonal decomposition. Then for each decomposition, we discuss some algorithms and realize the algorithm by Matlab software programming in the computer, and introduce the application of the algorithm of matrix in the multivariate regression analysis.The presented algorithm in this paper is the base of the analysis of multiple regression on the calculation, it has the very vital significance for using multiple regression analysis to solve practical problems.Keywords: Matrix decomposition; Matrix transformation; The algorithm; Regression analysis目录摘要 (I)Abstract (II)第一章引言 (1)1.1本文的研究背景 (1)1.2本文的主要工作 (1)第二章矩阵的三角分解及其算法 (2)2.1矩阵的LR分解及其算法 (2)2.2正定阵的Cholesky分解及其算法 (6)第三章矩阵的正交-三角分解及其算法 (10)3.1 Householder变换 (10)3.2 Givens变换 (17)3.3 Gram-Schmidt正交化及其修正算法 (20)第四章矩阵的正交分解及其算法 (24)4.1对称阵的谱分解及Jacobi算法 (24)4.2矩阵的奇异值分解及其算法 (28)第五章矩阵算法在多元回归分析中的应用 (31)5.1多元线性回归模型的参数估计与假设检验 (31)5.2基于Cholesky分解的回归算法 (33)5.3基于Householder变换的回归算法 (35)5.4谱分解在岭回归估计中的应用 (37)5.5总结 (41)附录 (43)参考文献 (57)致谢 (58)第一章引言1.1本文的研究背景数理统计方法是以概率论为基础,通过样本来了解和推断总体统计特性的科学方法,内容及其丰富.随着计算机使用的日益广泛,为了更好地应用数理统计方法来解决实际问题,从事统计工作或实际工作的人们都很关心如何应用计算机来更快完成各种统计数据的分析处理工作.故而出现了“统计计算”(Statistical Computation)这个方向.统计计算是数理统计、计算数学、和计算机科学三者的结合,它是一门综合性学科.在科学研究和生产实际的各个领域中,普遍地存在着大量数据的分析处理工作.如何应用数理统计学中的回归分析、多元分析、时间序列分析等统计方法来解决实际问题,以及如何解决在应用中的计算问题,对实际工作者来说是极需解决的问题.而矩阵算法对于实际工作者应用计算机实现数理统计中的计算问题具有很关键的位置.我们知道在多元统计分析的计算中观测数据一般用矩阵表示,对数据的分析转化为对数据矩阵的分析计算问题.如线性方程组的求解,矩阵的分解,奇异值的分解,特征值和特征向量的计算,广义特征值的计算,广义逆矩阵的计算及扫描变换等.这些常见的矩阵计算问题都是统计计算中经常遇到的问题.本文研究的目的是力求把统计计算中常用的矩阵算法的思想、步骤及其在计算机上的实现结合起来,为实际工作者应用多元统计分析解决实际问题提供便利.1.2本文的主要工作本文第2-4章主要讲解矩阵的三大分解,包括三角分解,正交三角-分解,正交分解.其中三角分解中主要讲解了Doolilttle、LDU以及正定阵的cholesky分解;正交-三角分解主要讲解了Housholder变换、Givens变换、GS正交化以及修正的MGS正交化;正交分解则主要讲解了对称阵的jacobi变换,以及一般矩阵的奇异值分解.本文的最后一章首先介绍了多元回归分析的一些基本理论知识,然后介绍了矩阵三大分解在多元回归分析中的运用.最后在本文的附录中我们给出了矩阵三大分解的中常用矩阵算法的程序以及在多元回归分析中的运用程序.第二章 矩阵的三角分解及其算法我们知道用有回代的消去变换法求解线性方程组b Ax =的过程,实质上就是化系数矩阵A 为上三角矩阵的过程;系数矩阵为上三角形矩阵的线性方程组用回代法很容易求得方程组的解.将一个矩阵分解为三角形矩阵或其他简单形式的矩阵,是矩阵计算的一种基本的方法.本章讨论将一个矩阵分解为两个三角形矩阵乘积的方法.2.1矩阵的LR 分解及其算法(一)矩阵的LR 分解(Doolittle 分解) 1.LR 分解的存在唯一性用高斯(Gauss )消去变换求解n 阶线性方程组b Ax =时,记增广矩阵为()b A ,化A 为上三角形矩阵的过程为:R a 00a a 0a a a A a a 0a a 0a a a A A )n (nn )1(n 2)1(22n 11211)1(nn )1(2n )1(n 2)1(22n 112111∆=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=→→⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=→ )()(n ,其中R 为上三角矩阵.初等变换的过程可通过初等变换阵i P 的运算来表示.记()n 211e ,,e ,t P =,其中n),2,1i ( )001,0,,0(e ,)a a ,,a a ,1( t T i i T111n 11211 ,,,,,个元素第==--= 则(). a A P A)1(ij011==)()(一般地记 ()1)n ,1,2,(i e ,,e ,t ,e ,,e P n 1i i 1i 1i -==+- ,其中i P 的第i 列T)1i (ii)1i (ni )1i (ii )1i (i),1i (i a a ,,a a ,1,0,,0t ⎪⎪⎭⎫⎝⎛--=----+ 是由)(1i A -的第i 列元素定义的向量(前1i -个元素为0),i P 是单位下三角矩阵,则A P P P P R 122n 1n --=,LR R P P P A 11n 1211∆----==这里11n 1211P P P L ----= 是单位下三角矩阵.因此对n 阶方阵进行1n -次高斯消去变换,就实现对A 的三角分解:LR A =.我们称矩阵A 分解为单位下三角矩阵和上三角矩阵的乘积的分解法为矩阵的LR 分解或Doolittle 分解.线性方程组b Ax =的解法可化为:⎩⎨⎧==⇔=⇔=(2)y Rx(1) b Ly b LRx b Ax其中(1)可求出b L y 1-=,然后由(2)用逐步回代法求出.x 因此,求解线性方程组的高斯消去法实质上就是系数矩阵的LR 分解.由上面的介绍我们大概了解了矩阵的LR 分解,那么我们不禁要问矩阵A 的LR 分解是否一定存在,若存在是否唯一?对此,我们先给出一个例子来进行讨论,如下:例2.1 设.3210A ⎥⎦⎤⎢⎣⎡= 显然02A ≠-=,方程组b Ax =有唯一解,但A 不存在LR 分解.事实上,使得等式⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡=3210r 0r r 1l 01LR 22121121 成立的L 与R 不存在.如果把A 的第一、第二行交换位置,得⎥⎦⎤⎢⎣⎡=1032A 1, 则1A 存在LR 分解.103210011032⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡ A 不存在LR 分解的原因是0a 11=(即A 的一阶顺序主子式为0).定理2.1 设A 为n 阶方阵,A 的k 阶顺序主子式记为.k ,,2,1k ,,2,1A a a a a d kk1k k111k ⎪⎪⎭⎫⎝⎛==∆则A 的LR 分解存在且唯一的充要条件是().1n ,1,2,k , 0d k -=≠对更一般的矩阵,如A 是退化方阵,或者A 是m n ⨯矩阵,是否仍有LR 分解?定理2.2 设A 为m n ⨯矩阵,()m n,min r r(A)≤= ,如果 , 0d k ≠,1k (=r),2, ,则A 存在LR 分解:LR A =(但不一定唯一).在矩阵A 中,当1r m n +==时,A 的LR 分解是存在唯一的;当A 是非奇异矩阵但不满足各阶顺序主子式不等于0,这时先对A 作行变换,然后进行三角分解.定理2.3 设A 为非奇异的n 阶方阵,则存在行置换阵P ,使LR PA =. 此定理对应于选主元的高斯消去变换,1L -实际上就是一系列初等变换阵的乘积 .注:定理2.1、2.2以及2.3的证明请参考数值计算方法的有关书籍,如参考文献[4]、[5].2. LR 分解的算法算法2.1.1(Doolittle 分解的算法) 已知矩阵()nn ija A ⨯=,设A 存在LR 分解,即⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡nn n 222n 112112n 1n 21mn m2m12n 22211n 1211r 00r r 0r r r 1l l 01l 001a a a a a a a a a.r r l r l r l r l r l r r l r r l r l r r r nn n 22n n 11n 222n 121n 111n n2n 1212212211121n11211⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡++++++=(2.1)利用(2.1)式两边矩阵对应元素相等计算出L 与R 的元素,具体步骤如下: (1)对比矩阵的第一行元素相等,求得:).n ,,2,1j ( a r j 1j 1 ==(2)对比矩阵的第一列元素相等,求得:n),,2,1(i r /a l 111i 1i ==(3)对比第二行、第二列元素相等,求得:⎩⎨⎧=-===).n ,,4,3i ( r /)r l a (l ).n ,,3,2j ( r l a r 22121i 2i 2i ij 21j 2j 2 — (4)以此类推,得到这一过程的递推计算公式如下:对n ,,1,2k =⎪⎩⎪⎨⎧+=-=+==∑∑-=-= ).n ,,1i j ( r /r l a l ).n ,,1i ,i j ( r l a r 1k 1t kk tk it ik ik 1k 1t tj kt kj kj )(— (2.2) 注:类似于Doolittle 分解,对于三角分解LR A =,当L 为下三角矩阵,R 为单位上三角矩阵时,该分解称为Grout 分解.Grout 分解的算法与Doolittle 分解类似,这里不再从复累述.下面我们介绍矩阵的另一种三角分解. (二)矩阵的LDU 分解由LDU R LDD LR A 1∆-===,其中()nn 2211r ,,r ,r diag D =,R D U 1-=为单位上三角阵,记⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=100u 10u u 1 U 2n n 112 则 ().n ,1,i j n;,1,2,i , r /r u ii ij ij +===例2.2 设⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=064010032A ,求A 的LDU 分解.解 首先由(2.2)式可得A 的唯一LR 分解式为.000010032102010001LR A ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡==则 ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡==00001003/21000010002102010001L D U A ,注:此例告诉我们只要知道了矩阵A 的LR 分解,相应的LDU 分解也就随之求出来了.另外,此例中A 是退化矩阵,但A 的LR 分解存在且唯一,故而A 的LDU 分解也是存在且唯一.2.2正定阵的Cholesky 分解及其算法以上介绍了一般矩阵A 的三角分解,当()ij a A =为正定阵时,A 的三角分解具有特殊的形式:TGG A =(其中)(ij g G =为下三角矩阵).正定阵A 的这种形式的分解称为Cholesky 分解.在多元统计分析中,涉及到的矩阵如协差阵、相关阵等一般都是正定阵.Cholesky 分解在统计计算中更是一类重要的、常用的矩阵计算. (一)Cholesky 分解的存在唯一性定理2.4 设A 为n 阶正定阵,则A 的Cholesky 分解TGG A =必存在;当G 的对角元素均取正时,分解式是唯一的.证明 A 是正定阵,所以A 的k 阶顺序主子式().n ,1,k , 0d k =≠由定理2.1知,A 存在唯一的LR 分解式,因此其LDU 分解式也存在唯一,即LDU A =.又因为A 是对称的,所以A A T=,即LDU DL U L D U T T T T T ==.由分解式唯一可知L U T=,于是()T T2/12/1T 1/21/2T GG LD LDL D LD LDL A ∆====其中当()nn 22112/1r ,,r ,r diag D=时,分解式T GG A =是唯一的.2/1LD G =为下三角矩阵. [证毕] (二)Cholesky 分解的算法以下给出正定阵A 的Cholesky 分解的两种算法. 算法2.2.1(直接递推算法) 设()ij a A =已知,TGG A =可写为:⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡nn 2n 221n 2111nn 2n 1n 222111nn 2n 1n n 222211n 1211g 00g g 0g g g g g g 0g g 00g a a a a a a a a a (1.3)(1)对比(1.3)式两边第一行元素对应相等,可得:()⎪⎩⎪⎨⎧===n ,2,3,j g /a ga g 11j 1j11111 (2)对比(1.3)式两边第二行元素对应相等,可得:()⎪⎩⎪⎨⎧=-=-=n ,2,3,j g /)g g a (gg a g 221j 21j 2j22212222 例2.以此类推,可得出计算G 的直接递推公式为(对n ,,1,2i =):()()⎪⎪⎪⎩⎪⎪⎪⎨⎧-==+=⎪⎭⎫ ⎝⎛-=-=∑∑-=-=1i ,1,2,j 0g n ,1,i j g /g g a g g a g ji ii 1i 1k jk ik ij ji 1i 1k 2ik ii ii (1.4)算法2.2.2(顺序Cholesky 分解算法或称平方根分解算法)设TGG A =,记()n g ,,g ,g G 21 =,其中() g ,,g ,0,,0g Tni ii i , =()n ,1,2,i =,则∑===n1i Ti i Tg g GG A (1.5)该算法是根据i g 的特点(前1i -个元素为0)及关系式(1.5)依次求出n 21g ,g ,g 的算法.记. a A A (0)ij (0))(∆== 具体步骤如下:例2.令)((1)ij T nn T 22T11(0)(1)a g g g g g g A A∆=++=-= . 由于n 2g ,g 的第一个元素均为0,所以)(1A的第一行和第一列元素全为0,即⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡---------)1(nn )1(2n )1(n 2)1(222n1)0(nn 211n )0(n2111n )0(n1121)0(2n 221)0(221121)0(211n 11(0)1n 2111)0(12211(0)11a a 0a a 0000g a g g a g g a g g a g a g g a g g a g g a g an 对比上式两端,求得1g 及(1)ij a :⎪⎩⎪⎨⎧=-=-===n),,2j (i, a /a a a g g a a n),1,2,(j a /a g )0(11(0)1j (0)1i )0(ij 1j i1)0(ij (1)ij )0(11(0)1j j1 (2)令)((2)ij T nn T 44T 33T22(1)(2)a g g g g g g g g A A∆=+++=-= ,类似的,由)(2A 的前两行和前两列元素均为0,可求得2g 及(2)ij a :⎪⎩⎪⎨⎧=-===n),,3j (i, a /a a a a n),2,3,(j a /a g )1(22(1)2j (1)2i )1(ij (2)ij )1(22(1)2j j2 例2.依次类推,1k -步后得出1k 21g ,,g ,g - 及1)(k A -,求k g 及(k)A 的计算公式为(对n ,1,2,k =):⎪⎩⎪⎨⎧+=-=+==------n),,1k j (i, a /a a a a n),1,k k,(j a /a g )1k (kk 1)(k kj 1)(k ki )1k (ij (k)ij )1k (kk 1)(k kj jk (1.6) (4)令()n g ,,g ,g G 21 =,则TGG A =就是A 的Cholesky 分解式.注:观察(1.6)式可知,用算法2.2.2求A 的Cholesky 分解式的过程实际上就是一种用初等变换化A 为特殊上三角矩阵的过程.例2.3 用平方根分解算法求矩阵⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=382910295861064A 的Cholesky 分解式. 解 对3,2,1k =利用递推公式(1.6)即可求出321g ,g ,g ,从而得到A 的Cholesky 分解式.下面我们来用初等变换法求A 的Cholesky 分解式.>==+⨯-+⨯-⨯⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡==TT 131211(0)(2,3,5)(4,6,10)4g r r 410 r r 46 r 41 382910295861064A A 其中,,27,0)14,49,0(49gr r 4914r 491 1314014490532A TT 23221>==+⨯-⨯⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=),(其中,)(>==⨯⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=TT 332)3,0,0()9,0,0(9gr 91900270532A 其中)(⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=300270532A )3( 记)3(T A G =,则T GG A = 由例2.3也可以看出用递推公式(1.6)求解矩阵A 的Cholesky 分解式实质上就是一种用初等变换化A 为特殊上三角矩阵的算法.第三章 矩阵的正交-三角分解及其算法矩阵A 的LR 分解,实质上是对A 作初等变换,化A 为上三角矩阵的过程.其实,我们还可以对A 作正交变换,即存在正交阵),,(n ,21i P i=,使得 R A P P P P 121-n n = (R 为上三角形矩阵), 则QR A =,其中112n P P P Q -=)( 是正交阵.我们称矩阵A 分解为正交阵Q 和上三角矩阵R 的乘积的分解法为矩阵的正交-三角分解,或简称QR 分解.通过正交变换化矩阵A 为上三角矩阵的这种变换思想在统计计算中是非常重要的.下面我们将介绍几种简单有效的正交变换方法.3.1 Householder 变换 (一)r Householde 矩阵定义3.1 设e 是欧式空间nR 中的单位向量,形如Tee 2I H -=的n 阶矩阵称为r Householde 矩阵(或r Householde 变换).性质 设H 为n n R ⨯中的一个r Householde 矩阵,则(1)H H T= (对称性); (2)T 1H H=- (正交性);(3)I H 2= (对合性).将Tee 2I H -=代入以上三式即可证得结论.为了进一步了解r Householde 矩阵,我们有必要探讨一下r Householde 矩阵的几何意义.当2n =时,x 与Hx 的关系如图3-1所示:如果把L 看成一面镜子,e 看成这面镜子的单位法向量,则Hx 即为x 在镜面L 下所成的像.因此,r Householde 矩阵的几何意义可描述为:把e 看成是nR中的一个1n -维超平面的单位法向量,对任意nR x ∈,Hx y =即为x 在1n -维超平面下所成的镜像.根据其几何意义,r Householde 矩阵也称为反射(镜像)矩阵或反射(镜像)变换.如果有nR u ∈,将u 看成是nR 中一个1n -维超平面的法向量,那么r Householde 矩阵又可写为T 12T 2Tuu I u u u u 2I ee 2I H -∆-=∙-=-=β (3.1)其中.u 21221=-β(3.1)式告诉我们,只要给定一个超平面的法方向,我们就可以通过这个法方向去定义一个r Householde 矩阵.再结合图3-1可以看出,法方向e 与y x -是共线的.所以,我们有以下定理.定理3.1 设nR 中有非零向量y x ≠且22y x =,则必存在rHouseholde 矩阵H ,使得y Hx =.证明 令2y x y x e --=,e 为与y x -共线的单位向量.设T2ee I H -=,则)y x ()y x ()x y x x )(y x (2x )y x ()y x (x )y x )(y x (2x x )2ee I (Hx T T T TT T-----=-----=-=由于22y x=,即y y x x T T =,所以x )y x 2(x y y x y 2x x )y x ()y x (T T T T T T -=+-=--从而,y.)y x (x Hx =--= [证毕]下面我们将r Householde 矩阵运用于矩阵的正交-三角化. 化A 为上三角矩阵的第i 步要求把)(1i A -的第i 列化为T(i)ii (i)1i ,0),0,a ,(a 且)(1i A-的前1i -列不变.因此,我们的想法是:若给定一个n 维向量T n 21)x ,,x ,(x x =,是否存在r Householde 矩阵i H ,使得y.,0),0,x ,x ,,(x x H T(1)i(1)1i (1)1i ∆-==事实上,由定理3.1可知,如果满足22x y=(记为条件1),那么我们便可构造一个r Householde 矩阵i H ,使得y x H i =(只需令y x u -=即可).但由于化A 为上三角矩阵的过程中,第i 步我们不仅要求将第i 列化为T (i)ii (i)1i ,0),0,a ,(a ,而且必须保证)(1i A -的前1i -列不变,因此构造的i H 应该形如⎥⎦⎤⎢⎣⎡*=-00I H 1i i (记为条件2),这样才能保证)(1i i A H -的前1i -列不变. 根据条件1和条件2,我们可以令(),0),0s x sign ,x ,,(x y i i 1i 1 ⋅-=-, 则()),x ,,x ,s x sign x ,0,(0,y x u n 1i i i i +⋅+=-= .u u I 0I uu I H T i i 1i 1i n 1T1ii ⎥⎦⎤⎢⎣⎡-=-=-+---ββi 其中().)x ,,x ,s x sign x (u ,s x s u 21 ,x s Tn 1i i i i i i i 2i 22i ni j 2ji +=⋅+=+===∑β这样选择的y ,即满足了条件1(22x y=),又满足了条件2.而对于y 的第i 个元素的符号为什么选择“()i x sign -”,那是因为当计算u 的第i 个元素()i i i s x signx ⋅+时,以“()i x sign ”作为i s 的系数可以避免i x 与i s 作减法运算从而防止计算结果损失有效数字.(二)利用r H o u s e h o l d e 变换的QR 分解设A 为n 阶非奇异矩阵,记A A (0)=,对n ,1,2k =计算)(1k k (k)AH A -=,其中k H 是由)(1k A-的第k 列向量定义的r Householde 矩阵,则R HA A H H H A 012n (n)===)(即QR A =,其中T1H H Q ==-为正交阵.对k 用数学归纳法,利用i H 的性质容易证明上述结论.以上结论说明非奇异矩阵存在QR 分解,并且通过r Householde 变换能够得出QR 分解式.当A 为m n ⨯矩阵(n m ≤),且m r(A)=(即A 为列满秩矩阵)时,类似地,对m ,1,2,k =计算)(1k k (k )AH A-=,⎥⎦⎤⎢⎣⎡===0R R HA A1(m),从而QR A =,1R 为m 阶上三角矩阵,T 1H H Q ==-为n 阶正交阵.若记[]21Q Q Q =(1Q 为m n ⨯正交阵,2Q 为m)n (n -⨯正交阵),则矩阵A 的QR 分解式还可以写为:[]11121R Q 0R Q ,Q QR A =⎥⎦⎤⎢⎣⎡==定理3.2 设A 为m n ⨯列满秩矩阵(n m ≤),则A 可以分解为:QR A =,其中Q 为m n ⨯正交阵,R 为m 阶上三角矩阵.如果规定R 的对角元素取正时,分解式是唯一的.证明 由m r(A)=知,A A T为m 阶正定阵,利用定理 2.4知,A A T存在Cholesky 分解:R R A A T T =(R 为m 阶非奇异上三角矩阵).记1AR Q -=,则I AR A R Q Q 1T 1T T ==--)(,即Q 为m n ⨯正交阵,且QR A =;即A 存在QR 分解.当规定R 的对角元素取正时,A A T 的Cholesky分解唯一,从而A 的QR 分解也唯一. [证毕]设A 为m n ⨯矩阵,n}.min{m,r r(A)≤=对A 作QR 分解的思路如下:(1)对A 作列变换,使A 的前r 列线性无关.即存在列置换阵P ,使[]10A A AP =,其中0A 为r n ⨯列满秩矩阵,而1A 可由0A 的r 个列向量线性表出,即存在)r m (r -⨯矩阵B ,使得B A A 01=,此时[]. B A A AP 00 =(2)对0A 作QR 分解,由定理3.2知,存在H 使得⎥⎦⎤⎢⎣⎡=0R HA 0,R 为r 阶非奇异上三角矩阵.故.00RB R B)HA (HA HAP 00⎥⎦⎤⎢⎣⎡== 在实际计算时,列置换不必在对A 作r Householde 变换之前执行,这一步可以顺便在执行r Householde 变换的过程中实现.下面我们根据这一思路来导出求解R 的具体步骤如下:记A A0=)(,对m ,1,2,k =计算(不妨设n m ≤):(2)由)(1k A-的第k 列计算k H 阵.()().u u I 0I H )a ,,a ,h (u s a sign a h ,s a s u 21 ,a s T k k 1k 1k n 1k k T)1k (nk)1k (k 1,k k k k )1k (kk)1k (kk k k )1k (kk 2k 22k k nkj 2)1k (jkk ⎥⎦⎤⎢⎣⎡-==⋅+=+===-+----+---=-∑ββ, 当出现0s k =时,则0a a a 1)(k nk1)(k k 1,k 1)(k kk====--+- ,此时)(1k A -的第k 列可由前1k -列线性表出.把第k 列移到最后一列,即令)m ,k (I A A 1k 1k )()(--=,其中)m ,k (I 表示交换m 阶单位矩阵m I 的第k 列与第m 列所得到的初等置换阵.最后对新的)(1k A-的第k 列(即原)(1k A-的第1k +列)重新计算k H 阵.(3)计算.Ak )(⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡-=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡-==--∆---------+---)()()(k 22)1k (121k )1k (22T k k 1k )1k (22)1k (121k )1k (22)1k (121k T k k 1k 1k n 1k 1k k k A 0A R A u u A 0A R A 0A R u u I 00I A H A ββ (3.2)由(3.2)式可知,我们只需计算 .A u u A A )1k (22T k k 1k )1k (22)k (22----=β (3.3) 为此,我们记)a ,,a (A )k (n )k (k )k (22 =,)a ,,a (A )1k (n )1k (k )1k (22---= .第一步:计算)k (k a ,由(3.3)式知)1k (kT k k k)1k (k)k (k a u u 1a a ---=β 由于()()()())s a s i gn (a s a sign s a sign a s a s s a sign a h k )1k (kk )1k (kk k )1k (kk k)1k (kk )1k (kk k)1k (kk 2k k )1k (kk )1k (kk k k⋅+⋅⋅⋅+=+⋅+=--------β 即()k)1k (kkk ks a sign h -=β 所以()()()()()k )1k (kk2k )1k (kk k )1k (kkk)1k (kk )1k (kkn k j 2)1k (jk )1k (kk k )1k (kk k )1k (kk )1k (kkn1k j 2)1k (jk )1k (kk k kk)1k (kk)1k (kT k kk)1k (kk)k (kks a sign )s a s a sign (s a sign a))a (a s a sign (s a sign a ))a (a (h h aau h aa⋅-=+⋅⋅-=+⋅⋅-=+-=-=-----=-----+=-----∑∑ββn),1,k (i 0)s s a (aa))a (a (h a aau a aa2k k )1k (kk k)1k (ik)1k (ik1k t 2)1k (ik )1k (kk k k)1k (ik)1k (ik)1k (kT k k)1k (ik)1k (ik)k (ik+==+-=+-=-=---+=-------∑βββn第二步:计算)k (j a ,由(3.3)式知)1k (jT k k k)1k (j)k (j a u u 1a a ---=β,)n ,,1k j ( +=()()()())1k (kk)1k (jT )1k (kk)1k (kk)1k (tjkt )1k (tk )1k (tj k t )1k (tk )1k (kj k )1k (kk k)1k (kk )1k (kj )1k (tj 1k t )1k (tk )1k (kj k kk)1k (kj)1k (jTk kk)1k (kj)k (kja aasasign a a )a a a s a sign (s a sign a)a a a (h h aau h aa-----=--=------+=-----=⋅-=+⋅⋅-=+-=-=∑∑∑nn nββ()()()n),1,k j (i, )a a (h a a )s a sign a a a (h a a)a a a s a (sign h s a sign a aa u a a a)k (kj )1k (kj k)1(k ik )1k (kj k)1k (kk )1k (tjk t )1k (tk)1k (kj k )1(k ik )1k (kj)1k (tj k t )1k (tk )1k (kj k )1k (kk k k )1k (kk )1(k ik )1k (ij)1k (jT k k)1(k ik)1k (ij)k (ij+=--=⋅+-=+⋅⋅⋅⋅-=-=-----=-----=---------∑∑nn β通过以上两步,我们便得到计算)k (22A 中各元素的计算公式:()() n),1,k j (i, )a a (h a aa n),1,k (j aaaa n),1,k (i 0a s a sign a )k (kj )1k (kj k)1(k ik )1k (kj)k (ij )1k (kk)1k (jT)1k (k)k (kj )k (ik k )1k (kk)k (kk ⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧+=--=+==+==⋅-=------- (3.4)经过m 次r Householde 变换后,得R A(m)=(m n ⨯上三角矩阵),令T 12m )H H H (Q =,则QR.A =算法3.1.1(用r Householde 变换化A 为上三角矩阵) (1)输入矩阵A ,并置A A 0=)(,对m ,1,2,k =,重复(2)~(5);(2)由)(1k A-的第k 列计算()∑=-=nkj 2)1k (jk k a s ,并置m k I P =;(3)判断0s k =是否成立:若0s k =成立,则令)m ,k (I A A1k 1k )()(--=,置)m ,k (I P P k k =,转第(2)步;若0s k =不成立,则转下一步(4)(4)由)(1k A-的第k 列计算k H 阵().u u I 0I H )a ,,a ,h (u s a sign a h ,s a s T k k 1k 1k n 1k k T)1k (nk)1k (k 1,k k k k )1k (kk )1k (kk k k )1k (kk2k k ⎥⎦⎤⎢⎣⎡-==⋅+=+=-+----+---ββ, (5)计算)(k A,首先令)()(1k k A A-=,再进行以下赋值()() m ,1,k j n ,1,k i )a a (h a aa m),1,k (j aaaa n),1,k (i 0a s a sign a )k (kj )1k (kj k )1(k ik )1k (ij)k (ij )k (kk)1k (jT)1k (k)k (kj)k (ik k )1k (kk)k (kk ⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧⎪⎪⎭⎫ ⎝⎛+=+=--=+==+==⋅-=------(6)重复(2)~(5),直至循环结束,令(m)A R =,令T12m )H H H (Q =,置m 21P P P P =,则QR.AP =用r Householde 变换化A 为上三角矩阵的算法是一种稳定有效的算法.3.2 Givens 变换对矩阵A 作正交变换化A 为上三角矩阵的另一种方法——Givens 变换,这是Givens 1954年提出的方法. (一)Givens 旋转变换阵首先考虑平面上的二维向量T 21)x ,x (x =,用正交阵⎥⎦⎤⎢⎣⎡-=θθθθcos sin sin cos G 左乘x ,则使得向量x 顺时针旋转θ角后变成Gx ,当取αθ=(见图3-2-1)时,⎥⎦⎤⎢⎣⎡=0s Gx .其中α满足:)x x s ( d s x sin c s x cos 222121+=⎪⎪⎩⎪⎪⎨⎧====∆∆αα 这样的矩阵G 称为平面旋转变换阵,它是正交阵,且作用在向量x 上,可把向量x 简化.这种形式的矩阵推广到n 维空间nR 中,有以下定义.定义 3.2 设向量)w ,,(w w n 1 =,令2j2i w w s +=(不妨设j i <),s w c i=,sw d j =,则称n 阶矩阵 行第行第i k 1c d 1d c 1G ij ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-= 为Givens 变换阵很明显,Givens 阵ij G 是由向量w 的第j i,两个元素定义的,它与单位阵nI 只在第j i,行和j i,列相应位置的四个元素上有差别.性质 设ij G 是由向量w 定义的Givens 变换阵(j i <),则有以下性质: (1)ij G 为正交阵;(2)设Tn 21ij )u ,,,u (u w G =,则s u i =,0u j =,)j ,i k ( w u k k ≠=; (3)用ij G 左乘任一矩阵A ,A G ij 只改变A 的第j i,行;用ij G 右乘任一矩阵A ,ij AG 只改变A 的第j i,列根据ij G 的特点及矩阵乘法可立即得到上述结论. 为把向量w 化为后面i n -个元素为0,记w w 0=)(,i).n 1,2,(k w G w )1k ()k i (,i (k)-==-+其中k)(i i,G +是由1)(k w-定义的Givens 变换阵.利用性质(2)可知:T 1i 1)0()1i (,i i)(n i in )1i n (in i)(n )0,,0,s ,w ,,w (w G G G w G w -+----====, 其中.w s nij 2j ∑==(二)用Givens 变换化A 为上三角矩阵 设A 为m n ⨯列满秩矩阵,记).(a A A (0)ij 0==)( 则用Givens 变换化A 为上三角矩阵的步骤如下: (1)令)()()()()1(ij010121n ,1n 1(1)a A G A G G G A∆-=== ,其中1G 是由)(0A 得第1个列向量定义的一系列Givens 变换阵的乘积.则⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=)1(nn )1(2n )1(n 2)1(22)1(n 1)1(1211a a 0a a 0a a s A )(,其中. )(a s n 1k 2(0)k11∑==(4)一般地,令)()()()()k (ij 1k k 1k )1k (,k 1n ,k kn (k)a AG A G G G A ∆--+-=== , 其中k G 是由1)(k A-的第k 个列向量定义的一系列Givens 变换阵的乘积.则⎥⎦⎤⎢⎣⎡=(k)22k 12k(k)A 0A RA )(, 其中)m ,,2,1k (R k =是k 阶上三角矩阵.故⎥⎦⎤⎢⎣⎡====-0R R GA A G G G G Am 0121m m (m))( ,即m n n n m n R Q A ⨯⨯⨯=,其中()()()T1)(m m,mn T32n 2T21n 1TG G G G G G G Q +== 是正交阵.算法3.2.1(用Givens 变换化A 为上三角矩阵)设A 为m n ⨯矩阵,假设已进行到第1k -步正交变换,即用正交阵1k 21G ,,G ,G - 依次左乘(0)A ,得)a (A )1k (ij 1)(k --=,第k 步拟把1)(k A -中的第k 列对角元)(1k kk a -以下的元素消为0,故第k 步对n ,1,k i +=执行以下步骤: (1)令n k I G =,并计算2)1k (ik21)(k kk )a ()(a s --+=; (2)若0s =,则取0d ,1c ==;若0s ≠,计算sa c 1)(k kk -=,s ad 1)(k ik -=,并令行第行第i k 1c d 1d c 1G ki ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-= ,k ki k G G G =; (3)令0a ,s a (k)ik )k (kk ==;(4)对m ,1,k j +=;计算⎪⎩⎪⎨⎧+-=+=----.ca da a da ca a 1)(k ij )1k (kj )k (ij 1)(k ij )1k (kj (k)kj 用Givens 变换化A 为上三角矩阵的变换是一个正交变换,它常用于多元统计分析中;另外在增加或删去一个观测数据后的回归分析计算中,Givens 变换还有特殊的作用.3.3 Gram-Schmidt 正交化及其修正算法(一)Schmidt -Gram 正交化方法这是高等代数中我们很熟悉的正交化方法.设m 21a ,,a ,a 是n 维空间中m 个线性无关的向量,Schmidt -Gram 方法是求单位向量m 21e ,,e ,e 使得: (1)由m 21a ,,a ,a 生成的线性空间与由m 21e ,,e ,e 生成的线性空间相等,即)e ,,e ,e (L )a ,,a ,L(a m 21m 21 =; (2)m 21e ,,e ,e 两两正交. 具体步骤如下:令 11a =β,单位化得111e ββ=;11222e )e ,a (a -=β,单位化得222e ββ=;∑-=-=1m 1k k k m m m e )e ,a (a β,单位化得m m m e ββ=;则m 21e ,,e ,e 满足以上要求,且⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=m 2m 21m 121m 21m 2100)e ,a (0)e ,a ()e ,a ()e ,,e ,e ()a ,,a ,(a βββ即.R Q A m m m n m n ⨯⨯⨯=(二)Schmidt -Gram 算法(GS 算法)记)a ,,a ,a (A m 21 =,)q ,,q ,(q Q m 21 =,⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=mm m 222m 11211r 000r r 0r r r R ,假设m r(A)=,则由QR A =得:)q r q r ,,q r q r ,q r ()a ,,a ,a (m mm 1m 1222112111m 21+++=即⎪⎪⎩⎪⎪⎨⎧++=+==mmm 1m 1m 22211221111q r q r a q r q r a q r a (3.5)注意到m 21q ,,q ,q 是两两正交的单位向量,由(3.5)式第一个方程可求出111a r =及1q ;然后用T 1q 左乘(3.5)式第二个方程,可求出2T112a q r =及22r 和2q ; 依此进行下去最后可求出mm 1m r ,,r 和.q m 这一算法通常称为Schmidt -Gram 算法,简称GS 算法.算法3.3.1(GS 算法)已知)a ,,a ,a (A m 21 =且m.r(A)= (1)111a r =,1111r /a q =;(2)对m.,2,3,k = 利用(3.5)式第k 个方程计算⎪⎪⎪⎩⎪⎪⎪⎨⎧-=-=-==∑∑-=-=.r /)q r a (q q r a r )1k ,,2,1i ( a q r kk 1k 1i i ik k k 1k 1i iik k kk k T i ik 以上算法是将矩阵A 的m 个线性无关的列向量正交化为m 21q ,,q ,q ,故称为正交化方法.分析GS 算法将发现,此算法的计算是先求出(3.5)式第一个方程右边的11r 和1q ;然后求出12r ,22r 和2q ; 直到求出R 和Q 为止.由于Q 各列j q 是作为A 中各列的线性组合,因此一般说来,用GS 算法求解列正交阵Q 时舍入误差较大,以致影响m 21q ,,q ,q 的正交性.特别是将算法用于求解最小二乘解时,常出现不稳定现象.针对GS 算法的缺点,下面给出修正的Schmidt -Gram 算法. (三)修正的Schmidt -Gram 正交化算法(MGS 算法)MGS 算法也是利用方程组(3.5),只是它求解的顺序与GS 算法不一样,MGS 算法求解的顺序为:先求解(3.5)式右边的第一列,直到求出R 和Q 为止.记)m ,,2,1i ( a a i (1)i ==,由(3.5)式的第一个方程首先求得(1)111a r =,11(1)11r /a q =,用T1q 左乘(3.5)式的第m ~2个方程,利用m 21q ,,q ,q 正交性可得).m ,,3,2j ( a q r (1)j T 1j 1 ==这样(3.5)式中右端第一列相应的m)1,2,(j r j 1 =和1q 便全部求出来了.为了求(3.5)式中右端第二列,令m),2,3,(j q r a a 1j 1(1)j(2)j =-=, 则(3.5)式变为⎪⎪⎩⎪⎪⎨⎧++=+==mmm 2m 2(2)m 333223(2)3222(2)2q r q r a q r q r a q r a (3.6)(3.6)式的左边是已知的,用类似的方法可求出(3.6)式的右边第一列(也就是(3.5)式右边第二列).依此做下去,即可求出Q 和. R算法3.3.2(MGS 算法)已知)a ,,a ,a (A m 21 =且m.r(A)=记)m ,,2,1i ( a a i (1)i == (5)对 1.m ,2,1k -= ,计算(k)k kk a r =,kk (k)k k r /a q =, )m ,,1k j ( a q r (k)jT k kj +==, 令m).,1,k (j q r a a k kj (k)j1)(k j +=-=+ (6)由m m m 1m m 1)(m )1(m m(m)m q r q r a a =-=---可得(m)m mm a r =,.r /a q mm (m)mm = MGS 算法有时也会得到不太理想的正交阵,但经验表明,MGS 算法在求解最小二乘解时比GS 算法更加稳定. 故就一般问题而言,MGS 算法的计算精度及稳定性要比GS 算法好. 由此可见,计算顺序的不同,有时会大大影响算法的效果.第四章 矩阵的正交分解及其算法矩阵的正交-三角(QR )分解是用正交变换化矩阵A 为上三角矩阵.而矩阵的正交分解是基于矩阵的正交相似变换将矩阵A 化为更简单的对角矩阵.本章分别讨论对称矩阵、非奇异矩阵以及一般矩阵的正交分解问题.4.1对称阵的谱分解及Jacobi 算法(一)对称阵的谱分解定理4.1 设A 为n 阶实对称矩阵,则存在n 阶正交阵U ,使得),,d i a g (D AU U 21Tn λλλ ,==, (4.1)其中n λλλ,,21 ,为A 的特征值,若记)u ,,u ,u (U n 21 =,则i u 是相对于i λn),1,2,(i =的特征向量.另外,(4.1)式还可以写成另一种形式:∑===n1i T i i i Tu u UDU A λ, (4.2)并称(4.2)式为矩阵A 的谱分解式.n λλλ,,21 ,也称为A 的谱.定理4.1的证明可参考[1],这里不再一一陈述.实对称矩阵A 的谱分解的计算问题其实就是计算A 的特征值和特征向量的过程.下面我们介绍求实对称矩阵A 的特征值和特征向量的经典算法——Jacobi 方法.(二)Jacobi 算法在定义3.2中我们定义了Givens 变换阵,当θcos c =,θsin d =时,它实质上是相应平面上的一个旋转变换阵,记为)(G ij θ,它作用在向量w 上,能使w 的坐标)w ,(w 21旋转θ角后变为0). , w w (2221+下面我们将介绍另一种旋转变换阵,它作用在矩阵A 上之后,将使A 的j)(i,位置和i)(j,位置的元素简化为0.定义4.1 在平面旋转变换阵中,选择角度θ,使得)b ()(AG )(G B ij Tijij ∆==θθ (4.3)中的.0b b ji ij ==这样的旋转变换阵)(G ij θ称为Jacobi 变换阵,并记为).(J ij θ当用)(G ij θ左乘矩阵A 时,A 中只有j i,行的元素变化了,其他元素不变;当用)(G Tij θ右乘矩阵A 后,A 中只有j i,列元素变化了,其他元素不变.记θcos c =,θsin d =,则⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧≠=-+-==+-=++=≠-==≠+==j ).i,t (k, a b a)d c ()a a (cd b b a c cda 2a d b ad cda 2a c b )j ,i t ( da ca b b )j ,i t ( da ca b b kt kt ij22ii jj ji ij jj2ij ii 2jj jj 2ij ii 2iiit jt tj jt jt it ti it (4.4) 为使0b b ji ij ==,应选择合适的θ,使得0a )sin cos ()a a (sin cos ij 22ii jj =-+-θθθθ即yx a a a 22c o s 2s i n t a n 2jj ii ij ∆=-==θθθ (4.5)由(4.5)式可解得⎪⎪⎩⎪⎪⎨⎧+-==++==)yx y 1(21s i n d )y x y 1(21c o sc 2222θθ (4.6)其中ij a 2x =,jj ii a a y -=设A 为n 阶实对称矩阵,ij J 是Jacobi 变换阵.令)b (AJ J B ij Tijij ∆==,则B 有以下性质:(1)B 仍是对称阵;(2)0b b ji ij ==,且j)i,t (k, a b kt kt ≠=; (3))j ,i t ( a a b b 2jt 2it 2jt 2it ≠+=+;(4)∑∑∑∑=====n 1i n1j 2ij n 1i n1j 2ija b ,即变换后所有元素的平方和不变;(5)2ij a 2)A (E )B (E -=,其中)A (E 表示A 中所有非对角元素的平方和. 证明 很显然(1)、(2)是成立的,下面我们证明(3)、(4)、(5) (3)由(4.4)式及1d c 22=+知()()().a a da ca da ca b b 2jt 2it 2it jt 2jt it 2jt 2it +=-++=+所以结论成立;(4)由于相似矩阵具有相同的迹,即)A A (tr )AJ A J (tr TTij Tij =,所以∑∑∑∑=========n1i n1j 2ij TT ijTij T ijij T ijTij Tn1i n1j 2ija )A A (tr )AJ A J (tr )AJ J J A J (tr )B B (trb 即正交相似变换下,矩阵所有元素的平方和不变;(5)将∑≠=t k 2kt b )B (E ,∑≠=tk 2kt a )A (E , )j ,i t ( a a b b 2jt 2it 2jt 2it ≠+=+以及0b b ji ij ==代入∑∑∑∑=====n 1i n1j 2ijn1i n1j 2ija b 即可得到2ij a 2)A (E )B (E -= [证毕] 以上性质(5)表明,经Jacobi 变换后,矩阵非对角元素的平方和减少了.对B 继续实行Jacobi 变换,使非对角元素的平方和趋于零.这时A 将近似的化为对角矩阵D .这就是用Jacobi 变换化对称阵A 为对角阵的基本思想.算法4.1.1(Jacobi 算法) 记A A0=)(,n 0I U =)(,并给定精度ε,置1k =,实行以下步骤: (1)对())1(k ij1k a A--=)(,选非对角元素中绝对值最大者: )1(k ij ni j 1)1(k j i a max a 00-≤<≤-=; (2)由)1(k i i 00a -,)1(k j j 00a -和)1(k j i 00a -确定1k -θ,使得1)(k j j 1)(k i i 1)(kj i 1k 000000aaa 2tan2-----=θ;(3)利用(4.6)式计算出1k 1k cos c --=θ,1k 1k sin d --=θ及1k J -; (4)令())k (ijT1k 1k 1k k a JA J A∆---==)()(,利用(4.4)式可计算出)k (ij aT 1k 1k k J U U --=; (5)检验ε<≤<≤(k)ij ni j 1a max 是否成立,若成立,则停止计算;否则令1k k +=,重复(1)—(5).假设N k =时,已经满足了(5)的检验条件,这时认为)n ,,2,1i (a )N (ii i =≈λ是对称阵A 的特征值,)u u u (U U )N (n N 2)N (1N ,,,)( =≈为相应的特征向量组成的正交阵.Jacobi 算法是否收敛的于对角形的问题,是个很重要的问题.下面的定理给出了满意的结论.定理4.2 .0)A(E lim k k =∞→)(证明 设1)(k t t )1(k ija max a -≠-=vv,则有 .)1n (n )A (E n n )A (E )(a)1k (2)1k (2)1(k ij-=-≥---由Jacobi 变换性质(5)及上式可知,)(k A的非对角元素的平方和满足:).A (E ))1n (n 21()A (E ))1n (n 21( )1n (n )A (E 2)A (E )a (2)A (E )A (E (0)k 1)(k 1)(k 1)(k 21)(k ij 1)(k k --≤--≤--≤-=-----)( 由于1)1n (n 21<--,所以当∞→k 时,0))1n (n 21(k →--;)A (E (0)是一固定的常数,所以.0)A (E lim k k =∞→)( [证毕]。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

⎥⎥⎥⎥⎦

⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡n k kn n
k k n X X X X X X Y Y Y μμ
μβββ
21212222121211
11
1
1
1
⨯⨯⨯⨯n k k
n n X Y μ
β
1、多元回归中,根据最小二乘法得到:
='e X
2、多元回归的最小二乘法结果为:
Y X X X ''=-1)(ˆβ
3、从2中推出的两个常用的结论:
(1)多元中,
1ˆβ的计算:k
k X X X Y ββββˆˆˆˆ33221----=
(2)多元样本回归函数的离差形式:
根据3(1)及ki
k i i i X X X Y ββββˆˆˆˆˆ33221++++= ,有:
ki
k i i i x x x y βββˆˆˆˆ3322+++=
4、多元中
∑2
i
e
的矩阵表示:
e e e e e e e e e
n n i
'=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=∑
212
12)( 【 ⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛=n e e e e 21】
同理:
y y
y
i
ˆˆˆ2
'=∑,
y
y y
i
'=∑2

Y
Y Y
i
'=∑2,
Y Y
Y i
ˆˆˆ
2
'=∑ 等。

“变量值的平方和”可以写成“变量矩阵的转置矩阵*变量矩阵”
5、2
Y
n Y Y TSS
-'=的证明:
2
2
2
2
2
2
2
22)
2()(Y
n Y Y Y
n Y Y Y Y Y Y Y
Y Y Y Y TSS i
i
i
i
i
i
-'=-=+-=+-=-=
∑∑∑∑∑∑
或利用4的结论有:
2
2
))(()
()()(Y
n Y Y Y Y Y Y Y Y Y Y Y Y Y Y Y Y Y Y Y Y
TSS i
-'='+'-'-'=-'-'=-'-=-=
∑ 【⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=n Y Y Y Y 2
1,⎥⎥⎥⎥⎦
⎤⎢⎢⎢
⎢⎣⎡=Y Y
Y Y 】
6、
2
ˆY n Y X ESS -''=β
的证明:
2
2
22
2
2
ˆˆ2ˆ)ˆ2ˆ()
ˆ(Y n Y X Y
Y Y Y Y Y Y Y
Y Y ESS i
i
i
i
i -''=+-=+-=-=
∑∑∑∑∑β
【 】
7、从

∑∑∑=
--=
=
2
22
22
ˆ)
()ˆ
(i
i i
i y y Y Y
Y Y TSS
ESS R 到

∑∑∑
+++=
233222
ˆˆˆi
i
ki k i i i i y
y x y x y x R βββ 的证明:
证明上式,只需证明:∑2
ˆi
y
=
∑∑∑
+++i ki k i i i i y x y x y x βββˆˆˆ3322 即可。

根据3(2):
ki
k i i i x x x y βββˆˆˆˆ3322+++= ,将其写成矩阵形式: ⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥
⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡k kn n
n k k n x x x x x x x x x y
y y
βββˆˆˆˆˆˆ3232232221312121
1
)1()
1(1
ˆˆ⨯--⨯⨯k k n n x y β
Y X e X Y X e Y X Y X Y Y ''=''-''=-''=''='=β
βββ
β
ˆˆˆ)(ˆˆˆˆˆX Y
X Y
''='=>
=ββˆˆˆˆY
n Y
e
Y e Y
Y i
i
i
i i
i
==
-=
-=∑∑∑∑∑)(ˆ
矩阵形式:
βˆˆx y
= 【注意
x 矩阵与X
矩阵的不同,这里的】
∑∑∑
∑∑∑∑
+++=⎥
⎥⎥⎥
⎥⎦⎤
⎢⎢⎢⎢⎢⎣
⎡=⎥⎥⎥⎥⎦

⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=''=''-''=-''=''='=i ki k i i i i i ki i i i i k
n kn k k n n k
i y x y x y x y x y x y x y y y x x x x x x x x x y x e x y x e y x y x y y y ββββββββββ
βββ
β
ˆˆˆ)ˆˆˆ()ˆˆˆ(ˆˆˆ)(ˆˆˆˆˆˆ3322323
22
1213323122221
3
22
?
='e
00ˆˆˆ)()()(ˆˆˆˆ3233223322322
12
133********='=⎥⎥⎥⎥
⎥⎦

⎢⎢⎢⎢⎢⎣⎡'=⎥⎥
⎥⎥⎥
⎦⎤⎢⎢⎢⎢⎢⎣⎡---'=⎥⎥
⎥⎥⎥
⎦⎤⎢⎢⎢⎢⎢⎣⎡---'=⎥⎥⎥⎥
⎥⎦

⎢⎢⎢⎢⎢⎣⎡'=⎥⎥⎥⎥⎦

⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡'=''∑
∑∑∑∑∑∑∑∑∑∑∑∑
∑∑β
ββββββ
i ki i i i i i k i ki i i i i i i i k ki i i i i i ki i i i i n kn k k n n e X e X e X e X e X e X e X e X e X e X X e X X e X X e x e x e x e e e x x x x x x
x x x e x
根据
⎥⎥
⎥⎥

⎥⎦

⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=='∑
∑∑∑
000
00
321 i ki i i i i i i e X e X e X e X e X。

相关文档
最新文档