Microsoft PowerPoint - new第四讲.ppt [兼容模式]

Size: px

Start display at page:

Download "Microsoft PowerPoint - new第四讲.ppt [兼容模式]"

凉抽王
5 years ago
Views:

1 高等数值算法与应用 ( 四 ) Advanced Numerical Algorithms & Applications 计算机科学与技术系喻文健

2 内容概要矩阵分解及其应用六大分解简介矩阵与线性方程组求解基本理论 LU 分解, Cholesky 分解及其应用 QR Q 分解与线性最小二乘问题特征值分解奇异值分解 (SVD) 及其应用稀疏矩阵的直接解法

3 Spectrum Decomposition and Algorithms 3

4 本节主要内容矩阵的特征值与特征向量有关的计算问题特征值计算方法幂法反幂法 ( 应用 :PageRank) 收缩技术 QR 算法演示程序 eigshow 4

5 n n 矩阵 A 的特征值分解基本概念与重要性质 n n n Ax = x det( I A) = λ λ λ a jj λ + = j = 特征值 λ, 特征向量 x, 特征方程, 0 谱 λ(a), 特征子空间实矩阵 : 特征值可能是虚数, 它一定与其共轭成对出现实特征值对应实特征向量, 虚特征值对应非实特征向量实对称矩阵 A: 特征值均为实数, 可正交对角化 A = QΛQ AQ = QΛ λ ( A) = λ( A ), 同一特征值还有左特征向量 y A 若为对角阵上 ( 下 ) 三角阵, 特征值为其对角线元素若为分块对角阵分块上 ( 下 ) 三角阵, 特征值为对角块矩阵特征值的合并 = λy A A A m m A A ( A) = ( A ) m λ λ, A = j = jj Amm 5

6 特征值分解基本概念与重要性质 ( 续 ) 倍数 ca 的特征值, 平移 A+pI 的特征值幂 A k 的特征值, k 为正整数多项式 P(A) 的特征值, 为特征值的多项式函数值若 A 非奇异, A - 的特征值相似变换不改变特征值 B = X AX, λ( A) = λ( B) 设有 m 个不同的特征值 λ,, λ m, 若 λj 是特征方程的 n j 重根, 称 n j 为其代数重数, 其几何重数 d j 为特征子空间维数 j, nj dj ; 若两者总相等, 则有 n 个线性无关特征向量, 这种矩阵称为非亏损矩阵, 否则为亏损矩阵非亏损矩阵 A 可对角化, X AX = Λ 6

7 特征值分解计算有关的问题矩阵 A 是实的, 还是复的? 规模大小? 既要特征值, 又要特征向量? 要所有特征值, 还是一部分 ( 最大最小 )? 特征值 ( 向量 ) 计算问题的敏感性同一个矩阵的不同特征值 ( 向量 ) 的敏感性可能不同对非亏损矩阵, 设 X AX = Λ, 矩阵 A 变化 ΔA, μ 为扰动后矩阵的特征值, λ k 为与它最接近的 A 的特征值, 则 μ λk ( A) cond ( X) ΔA 若特征向量接近线性相关, 则特征值计算非常敏感对于正规阵 A (A A=AA, 特例是实对称阵 ), 有正交特征向量组,cond(X) =, y 不敏感为左特征向量 y x 亏损矩阵一般分析 : Δλ < ΔA = ΔA y x cosθ 7

8 特征值分解及其计算特征值计算的方法求解特征方程? 计算特征多项式系数工作量很大, 且对矩阵扰动敏感高阶多项式方程求根的工作量很大, 两种误差累计使解不准例 :, 比小一点, 准确特征值为 ±, 但, 算出两个根均为因此, 反而是求矩阵特征值的方法被用于求解多项式方程是 C n 的特征多项式因此, 通过求 C n 的特征值即求解了多项式方程幂法 (power iteration) 求最大的特征值 ( 主特征值 ) 及对应的特征向量条件 : 主特征值是唯一的, 对实矩阵它一定是实的

9 特征值特征向量的计算特征向量的计算幂法与反幂法幂法的应用 :Google 的 PageRank 算法加速幂法的收敛 ( 瑞利商位移技术 ) 反幂法求最小的特征值, 对 A - 应用幂法与位移技术结合, 即对使用, 精确求某一特征值计算所有特征值的主要方法. 若已知某个特征值及特征向量, 可采用收缩技术. QR 算法 3. 分治算法 4. Krylov 子空间迭代法 : Lanczos 算法, Arnoldi 算法

10 计算所有特征值的方法收缩技术?

11 计算所有特征值的方法

12 计算所有特征值的方法 QR 算法二十世纪十大算法之一, 计算中小规模矩阵所有特征值的稳定有效的方法实 Schur 分解 : S A = QSQ Q Q AQ = S 为分块上三角阵, 对角块为阶或阶方阵 QR 算法就是不断做正交相似变换, 得到 Schur 分解 A0 = A, QkRk = Ak, Ak+ = RQ k k, k=0,,, A Q AQ = 基本 QR 算法, 所以 A 的特征值与 A 的相同 k k+ = k k k 可证明, 序列 {A k } 基本收敛到拟上三角阵 S, 对角元按绝对值从大到小排列条件 : 实矩阵 A 为非亏损阵, 等模特征值只有实重特征值, 或多重复的共轭特征值两种情况实际应用, 先化简为上 Hessenberg 阵, 再使用位移技术

13 计算所有特征值的方法实用的 QR 算法 0 对上 Hessenberg 矩阵 A k 执行 QR 迭代对 A k 进行 QR 分解 : 算 A k+ : k+ 采用 Givens 旋转, 以 4 阶矩阵为例 A k PA PPA PPPA k k 3 k 3 k PPPAP A 3 k P k + PPPAP A k+ 仍是上 Hessenberg 阵. 随 k 0, 下三角元素 0

14 计算所有特征值的方法实用的 QR 算法一般矩阵化为上 Hessenberg 阵, 使用 Householder 变换最终为上 Hessenberg 阵好处 : 减小后续 QR 迭代每步计算量加速收敛过程演示程序 eigsvdgui 矩阵元素颜色什么意思?

15 计算所有特征值的方法实用的 QR 算法带位移的 QR 迭代迭代矩阵序列仍正交相似, 加快收敛速度, 使迭代过程对更多的矩阵能收敛简单策略取 s k 为 A k (n, n), 还有更复杂的策略和双位移技术实对称矩阵, Wilkinson 证明了单位移的 QR 迭代过程收敛 ; 实非对称矩阵, 还没有这样的证明, 理论上可能不收敛非对称矩阵对称矩阵

16 计算所有特征值的方法实用的 QR 算法的算法描述

17 Matlab 中有关命令及其他 eig 命令 - 主要对稠密矩阵 d=eig(a); 返回所有特征值, 使用 qr 算法 [V, D]=eig(A); 返回特征值和完备特征向量组, qr 算法若 A 稀疏且对称, 也可返回所有特征值, 对其他稀疏矩阵, 或要算特征向量, 则需使用 eigs 命令可能还会对矩阵元素进行比例调整, 保证数值稳定性 eigs 命令针对稀疏矩阵, Lanczos 算法, Arnoldi 算法 d = eigs(a, k); 返回最大的 k 个特征值, k 缺省值为 6 [V, D]=eigs(A, k); 返回最大的 k 个特征值与相应特征向量 hess, schur 命令

18 期末 Project 安排要求根据所学内容编写一个 Matlab 演示软件形象地表现概念算法或应用, 利于使用者学习巩固知识, 中文界面且交互性效率较好选题可根据 Heath 网站改编, 或根据教学内容自编由单独的.m 程序文件实现, 设计使用者的思考编程题特别欢迎有独特的创意 ( 若改编 NCM 需明确说明修改计划 ) 时间安排例 : floatshow.m, qformula.m, revste.m, matrix_cond.m 选题 : 第周之前 ( 月 3 日 deadline), 在课程讨论提交选题内容, 获得同意, 各人之间不要重复完成 / 演示阶段 : 第 6 周, 做课堂演示, 交实验报告占成绩比重 :40% 8

19 Singular Value Decomposition (SVD) 9

20 Singular Value Decomposition (SVD) 矩阵的奇异值分解 A 为 m n 矩阵, 则 A = UΣV U, V 均为正交阵, Σ=(σ ij ) 为 m n 的对角阵, 且对角元素 0, 称为奇异值记 Σ 的第 i 个对角元为 σ i, 通常调整 U, V 列的顺序使得奇异值按降序排列 : σ σ σ 3 = U Σ V = 只讨论 m n 的情况 Why?

21 Singular Value Decomposition (SVD) 奇异值分解有关性质 A 的分解 : A = UΣV A = VΣ U U, V 列向量为 u i, v i, 称为左右奇异向量 Av = σ u, 即 Av = σ u, Av = σ u, L 演示程序 eigshow i i i σ A = U Σ V A A = V Σ Σ V = V V n n σ n 假设 m n, 有 n 个奇异值所以奇异值是 A A特征值的算术平方根, 是唯一的存在性的证明 : 设 AA非零特征值为 σ,,, σ σr 形成 [ ] r O V AV = O AA= V V O O V Σ = r V AAV Σ r I r 令 U = AVΣ r, 则 UU = I r, 则可将 U 看作正交阵的 r 列 r

22 Singular Value Decomposition (SVD) 奇异值分解存在性的证明 A = UΣV [ ] r O V AA= V V O O V U AV U = = O AVΣ r 为正交阵 r 列根据补齐正交单位向量得到正交阵 U = [ U U ] U AV U = r = [ ] U A V V = U AV U AV = U AV O U U AV U AV U AV O AVΣ U AV = Σ V A AV r = Σ r = U AV U U O = Σ 所以, U AV = Σ r O = Σ SVD 的存在性得证! O O r

23 Singular Value Decomposition (SVD) 奇异向量的意义 = = A UΣV AV U Σ 对应于非零奇异值的矩阵 U 的列, u,, u r 是 span(a) 的单位正交基, U 的剩下的列是 span(a) 补空间的单位正交基 0 奇异值对应的矩阵 V 的列, v r+,, v n 是 A 的核空间 { x n : Ax = O } 的单位正交基与特征值分解的比较 (m=n) = = A XΛX AX X 特征值分解仅对非亏损阵存在, 非零特征值对应的矩阵 X 的列向量为 span(a) 的基, 但一般不正交当 A 为对称半正定时, 其特征值分解存在, 特征值 0 λ A = QΛQ 也是 SVD, 若正定, cond( A ) = max = σ max λ min σ 其他实对称矩阵, A = QΛQ 适当修改即得 SVD 分解 Λ min

24 Singular Value Decomposition (SVD) 奇异值与矩阵的 - 范数条件数 A A = UΣV, 其中 Σ=diag( σ, σ,, σ n ) 最大 Ax UΣV x ΣVx Σy = max = max = max = max = σ x O x O x x x O O Vx y y 若非奇异 (m=n), min Ax A Σy A = 其奇异值均 >0 = min = x O x /σ y O y n cond( A ) = σ σ σ 若 m>n, 仍有 A = σ. 若列满秩 ( 奇异值 >0), 奇异值 n σn + = = = A ( A A) A ( VΣΣV ) VΣU V ( ΣΣ ) Σ U A + = /σ n cond( A ) = σ σ n 其他情况, 规定 cond( A ) =

25 如何计算奇异值分解利用求特征值的 QR 算法等效于求 A A 矩阵的特征值实际上不直接计算 A A, 因为对 A 施加左右正交变换不改变它的奇异值所以, 先将 A 化简 (Householder 变换 ) 为简单的 B, 再对 B B 执行 QR 迭代算法 HA A = a LLan α β OO O βn αn H a HA = α e ( m) n α b HA 0 0 H α β a LLan a LLan ( n ) Hb = βe = H LH A L = B n B B k QR 迭代过程 : k 为三对角矩阵, 迭代中保持其形状演示程序 : eigsvdgui 5

26 SVD 的应用任意矩阵的 - 范数, - 条件数, cond( A ) max = σ σmin 确定矩阵的秩 : 不为零的奇异值的数目解线性最小二乘问题 Ax b A 列满秩, 奇异值均 >0 AV = UΣ span(a) 的维数 UΣV x b Σ V ΣV x U b x = U b = Σ A 列不满秩, 存在 0 奇异值 min U b- ΣV x 令 y=v x min σ y = u b, i =, L,r U b- Σy 取最小值要求 i i i 若 y 的其他分量 =0, 是 - 范数最小的解. 由于 y = x, r r u = = = i b 得到原问题最小范数解 : x Vy viyi vi x V U b i= i= σ i 6

27 SVD 的应用更一般的矩阵伪逆 ( 推广到列不满秩矩阵 ) 标量 σ 的伪逆 : /σ, 或 0 Σ + 对角阵的伪逆 : 矩阵转置, 然后每个对角元求伪逆一般的矩阵 A: A = UΣV, 则 A + = VΣ + U 易证明, 它与 A 列满秩时伪逆定义兼容 + = A = ( A A ) A, 任意情况下, Ax = b的最小 - 范数解为 x = A + b 其他伪逆的性质 : + + A = A ( ) + + A ( A ) = ( ) A + = / σ / r

28 SVD 的应用低秩矩阵近似 (Low rank matrix approximation) σ v v n A = [ u ] [ ] u σ u u E n = σ σ = σ i = vn vn 其中 E = u v, i =,, n, 它秩为, 奇异值 σ =, 其余 =0 i i i m n n i i 存储 Ei 只要 m+n 个单元, 计算 Ex i 只需 m+n 次乘法 k 忽略小奇异值对应求和项, 称 A = σi Ei Ei 为 A 的秩 k 近似, 可证明, 以 Frobenius 范数度量, 在所有秩为 k 的矩阵中, 它最接近 A 在图像处理数据压缩信息抽取等方面都有应用 i = 若 A 实对称, 从合同标准型导出类似公式, 其中 = 演示程序 svd_appl.m E i 更简单 PCA, PFA

29 SVD 的应用完全最小二乘问题 (otal least squares) 普通线性最 x 小二乘问题 :[ a a ] n = Ax = b x n a span(a) 假设 b的数据有误差, 而模型数据, 即 A准确, 数据点到曲线竖直距离最小 b 完全最小二乘问题 : 所有数据都有误差或不确定性, 包括 A 中的数据, 要求 b 和 A到拟合 â a 平面总的 - 范数距离最小 Âx âa 设 y = Ax ˆ a, x 为待求的参数, 则 span(a) ^ ˆ A y [ A b ] 解法 : 用 SVD 求 [A b] 的秩 n 近似列秩 : n (n+) ˆ ˆ ˆ A x = O 若 [ A b] = UΣV, 则 A y = UΣV y, ˆΣ 由将 Σ中 σ n + 变为 0 得到 x 与成比例 v n+ a b Ax

30 SVD 的应用完全最小二乘问题的例子用模型函数 f () t = xt t 拟合数据 : y 常规最小二乘问题拟合结果 : x = 05. 完全最小二乘拟合 : A b = 3 = [ ] vn+(: n) x = = ( n + ) v n+ v ˆ y ˆ [ b ] 所有过原点直线中, 数据点 tˆ yˆ t y 所有过原点直线中, 数据点到该线垂线距离平方和最小!

31 Matlab 命令与方法比较 Matlab 中 svd 命令 s = svd(x); [U,S,V] = svd(x); 返回奇异值, 或完整分解 [U,S,V] = svd(x,0); [U,S,V] = svd(x, econ ); 分解的简化形式,S 为方阵,U 或 V 为正交阵的一部分方法间比较解线性最小二乘问题, 正交变换法最稳定准确, 但计算量可能比法方程方法大 (Househoulder 变换, m>>n) 对 ( 近似 ) 秩亏损问题, 只能用列主元正交变换法, 或基于 SVD 的方法, 后者更可靠, 计算量则更大奇异值特征值分解联系紧密, 前者似乎更易用一些

32 Some useful Matlab commands () Constants: t realmin, realmax, eps, Inf, -Inf, NaN N Logical: isinf, isfinite, isnan Floating point arithmetic: single, double, vpa Operators: \ (back slash) Matrices: zeros, ones, eye, diag, rand, cond, condest, rank, det, inv, eig, tril, triu, sparse, issparse, speye, spones, spdiags, sprand, nnz, full, eigs Graphs: spy Algorithms: lu, chol, qr, eig, svd, hess, schur, polyfit 3

33 本节内容小结矩阵的奇异值分解的概念分解的存在性, 唯一性奇异向量的意义奇异值与 - 范数,- 条件数计算奇异值的算法奇异值分解的应用 - 范数条件数, 矩阵的秩解秩亏损的线性最小二乘问题, 最小 - 范数解一般的伪逆定义低秩矩阵近似完全最小二乘问题的求解 33

34 Assignment 阅读课本第三章, 第四章有关内容阅读文献 : ( 教学资源 -- top 0 algorithm ) B. N. Parlett, he QR Algorithm, Computing in Science and Engineering, Vol., No., pp , 000 作业题见网络学堂 34

슬라이드 1

슬라이드 1 08-09 年度第一学期 0050 00503 计算方法 (B) 童伟华管理科研楼 05 室 E-mail: togwh@ustc.edu.c 中国科学技术大学数学科学学院 http://math.ustc.edu.c/ 第七章计算矩阵的特征值与特征向量特征值与特征向量在实际工程计算中, 经常会遇到特征值和特征向量的计算, 如 : 机械结构或电磁振动中的固有值问题 ; 物理学中的各种临界值等