Gram-Schmidt正交化

格拉姆-施密特(Gram-Schmidt)正交化常用于求解向量空间的标准正交基,同时也是一种天然的求解矩阵的QR分解的方法,即将一个矩阵A分解为一个正交矩阵Q和一个上三角矩阵R的乘积,即A=QR。这里我们假设A是一个方阵,当然A不是方阵的时候也是可以进行QR分解的。QR分解可用于线性方程组的求解,并且使得求解的过程更加高效、稳定,这里不细说,我们重点关注Gram-Schmidt正交化的过程和原理。

Classical Gram-Schmidt

我们普遍在课本中主要接触到Gram-Schmidt正交化属于传统的正交化方法,其相对容易理解,但是在计算机中,由于数值计算的精度问题,其计算的误差累积较快,当矩阵维数较高的时候具有较差的稳定性。其误差影响两个方面,一是矩阵Q的正交性,矩阵的正交性使得其逆矩阵可以直接通过矩阵的转置得到,如果正交性被破坏,那么会造成计算的错误;二是矩阵R元素的精度问题,越往后面计算所得的元素误差越大。具体可以看一下下面的网页中的demo。总的来说传统方法的误差水平大约为计算机精度的平方根。

https://nbviewer.jupyter.org/github/mitmath/18335/blob/master/notes/Gram-Schmidt.ipynb

传统的Gram-Schmidt正交化(CGS)方法类似于拼积木,如下图的例子。其先拿出第一个向量作为基准,接着拿出第二个向量,并通过对照第一个基准找到垂直于第一个基准的第二个基准,并依次类推。因此,在构建第n个基准时,n+1及以后的向量是不可见的。因为在计算的过程中不可避免地引入误差,而后面的向量只有在其被计算的时候才可见,这时误差可能已经积累到一定程度了,所以总的来说其稳定性较差。

对于上图,我们可以将其表示为以下的计算过程,最终我们可以直接得到一个正交矩阵Q和一个上三角矩阵R,这也就是为什么Gram-Schmidt天然就是一种QR分解的方法。

Modified Gram-Schmidt

前面已经说了,传统的Gram-Schmidt正交化方法误差累积速度较快,原因在于在当前计算的向量之后的向量是不可见的,使得后面的计算使用的是已经积累了相当误差的结果。因此,一种改进的Gram-Schmidt正交化(MGS)方法被提了出来,其实际上与CGS是等价的,但其对计算的过程进行了重排,使得从一开始所有的向量都是可见的,这样大部分的计算都在误差尚未积累到较大程度的时候就已经被执行,如下图所示。实际上,MGS的计算量越往后是逐步减少的,因此前面计算积累的误差的影响就不会剧烈地扩散开,所以MGS可以使求得的矩阵Q的正交性更好,同时矩阵R的误差也被控制在计算机精度的水平。

与CGS不同,MGS的主要改进思路是,首先选定一个向量作为第一个基准,然后将其余所有向量都投影到该基准的正交空间中,这时第一个基准与正交空间中所有的向量都是垂直的。因此,在此之后,我们只需在该正交空间中,对剩下的向量重复前面的工作,那么最后所有的向量都是相互正交的了。同时,我们也可以得到相应的正交矩阵Q和上三角矩阵R

对于上图,我们可表示成以下的计算过程。最终我们得到的矩阵和上面的CGS所得到的其实是一样的。

改进的Gram-Schmidt正交化方法虽然相比于传统方法在矩阵Q的正交性方面虽然有了不少改善,但实际上还是有点不太理想的,具体可以看一下上面的那个链接。相比之下,Householder变换方法在正交性方法能够达到非常高的精度,当然对于求解R的精度与MGS也是相同的水平的,而且其计算复杂度通常会比MGS更低,因此目前在QR分解中更多使用的是Householder变换方法。除此之外,还有一种Givens变换方法。具体关于QR分解的知识可参考这个文档,这里不细说。

https://wenku.baidu.com/view/6872ac728e9951e79b892700.html

Python实现

# -*- coding: utf-8 -*-
import numpy as npdef cgs(A):m, n = A.shapeQ = np.copy(A)R = np.zeros([n, n], dtype='float64')for i in range(n):R[:i, i] = A[:, i].T.dot(Q[:, :i])Q[:, [i]] = A[:, [i]] - Q[:, :i].dot(R[:i, [i]])R[i, i] = np.linalg.norm(Q[:, i])if np.abs(R[i, i]) < 1e-15:Q[:, i] = np.zeros(m, dtype='float64')R[i, i] = 0.else:Q[:, i] /= R[i, i]return Q, Rdef mgs(A):m, n = A.shapeQ = np.copy(A)R = np.zeros([n, n], dtype='float64')for i in range(n):R[i, i] = np.linalg.norm(Q[:, i])if np.abs(R[i, i]) < 1e-15:Q[:, i] = np.zeros(m, dtype='float64')R[i, i] = 0.else:Q[:, i] /= R[i, i]R[[i], i+1:] = Q[:, [i]].T.dot(Q[:, i+1:])Q[:, i+1:] -= Q[:, [i]].dot(R[[i], i+1:])return Q, R

图解格拉姆-施密特正交化和改进的格拉姆-施密特正交化相关推荐

  1. matlab格拉姆施密特,改进的格拉姆-施密特正交化(modified Gram-Schmidt Process)

    最近在重新学习线性代数,学习的教材是MIT Gilbert Strang 教授的<INTRODUCTION TO LINEAR ALGEBRA>,在第4.4章节格拉姆-施密特正交化时,书中 ...

  2. C++实现Schmidt施密特正交化算法(附完整源码)

    C++实现Schmidt施密特正交化算法 C++实现Schmidt施密特正交化算法完整源码(定义,实现,main函数测试) C++实现Schmidt施密特正交化算法完整源码(定义,实现,main函数测 ...

  3. 施密特正交化_考研数学答疑210施密特正交化

    线代 同学的疑问: 这句话怎么去理解? 我作的答疑: 由施密特正交化的过程知,施密特正交化在实对称矩阵对角化中使用,在相同特征值下的特征向量之间使用.实矩阵A可施密特正交化化为对角阵,那么A一定是实对 ...

  4. matlab施密特正交化,浅谈压缩感知(十九):MP、OMP与施密特正交化

    浅谈压缩感知(十九):MP.OMP与施密特正交化 关于MP.OMP的相关算法与收敛证明,这里仅简单陈述算法流程及二者的不同之处. 主要内容: MP的算法流程及其MATLAB实现OMP的算法流程以及MA ...

  5. 【线性代数】施密特正交化方法——Python实现

    文章目录 思想 Python代码 思想 施密特正交化方法: 将n维子空间中的任意一组基向量变换成标准正交向量. 假设有两个向量a⃗\vec{a}a和b⃗\vec{b}b,若要使两向量正交,则a⃗\ve ...

  6. 矩阵论——施密特正交化求行列式QR分解

    矩阵论--施密特正交化&求行列式&QR分解 前言 施密特正交化(Schimidt Orthogonalization ) 代码实现 example 求行列式det 代码实现 examp ...

  7. Self-Orthogonality Module:一个即插即用的核正交化模块

    作者丨苏剑林 单位丨追一科技 研究方向丨NLP,神经网络 个人主页丨kexue.fm 前些天刷 arXiv 看到新文章 Self-Orthogonality Module: A Network Arc ...

  8. 字符串匹配之KMP(KnuthMorrisPratt)算法(图解)

    文章目录 最长相等前后缀 next数组 概念 代码实现 图解GetNext中的回溯 改进 代码实现 代码 复杂度分析 最长相等前后缀 给出一个字符串 ababa 前缀集合:{a, ab, aba, a ...

  9. 施密特触发器(Schmitt Trigger)?

    施密特触发器(Schmitt Trigger),简单的说就是具有滞后特性的数字传输门. (一)施密特触发器结构举例 (二)施密特触发器具体分析 (三)施密特触发器电路用途 (四)施密特触发器相关部分总 ...

  10. 什么是施密特触发器(Schmitt Trigger)?

    http://hi.baidu.com/hieda/blog/item/c996d9cc5d1a8c1400e92877.html 施密特触发器(Schmitt Trigger),简单的说就是具有滞后 ...

最新文章

  1. 如何写新的Python OP
  2. 设置自动关门时长_小米苹果全适配,绿米D100全自动指纹锁新鲜上手
  3. Dockerfile实践优化建议
  4. shell操作典型案例--FTP操作
  5. STM32 基础系列教程 5 – 系统定时器
  6. HTML DOM节点的属性获取
  7. 获取当前程序运行的主机名称
  8. Jmeter+ForEach控制器+BeanShell取样器+BeanShell PostProcessor爬取网站信息储存csv
  9. 【机器学习算法应用和学习_1】1.1 机器学习框架
  10. unix系统安装及应用
  11. linux DSA 开发上手笔记(一)
  12. Windows下让Tomcat6定时重启服务的方法
  13. [Bullet3]常见物体和初始化
  14. android apk 重新签名工具,安卓apk重新签名教程,快来定制自己的apk吧
  15. SNAP Java API处理Sentinel-1数据
  16. LeetCode-Python-275. H指数 II
  17. centos7 挂载 硬盘 shell 懒人系列-2
  18. PHP5.4发布:新特性与改动
  19. 关于微信投票的刷票分析
  20. js 创建标签 追加标签

热门文章

  1. CC2530F256RHAR 射频芯片 无线收发器芯片 ZigBee 解决方案
  2. java中类成员的限定词_JAVA中类成员的限定词有以下几种、private ,public ,
  3. NeRF神经辐射场学习笔记(二)——Pytorch版NeRF实现以及代码注释
  4. 台式电脑计算机怎么用,怎么用键盘开机电脑_台式电脑键盘怎么开机-win7之家
  5. 嵌入式视频处理考虑(二)
  6. 持续学习:(Elastic Weight Consolidation, EWC)Overcoming Catastrophic Forgetting in Neural Network
  7. DDI(DNS、DHCP和IPAM)解决方案的全球与中国市场2022-2028年:技术、参与者、趋势、市场规模及占有率研究报告
  8. 否则在Python中使用for / while语句
  9. 有没有集工作记录、项目时间线于一身的便签软件?
  10. 线性代数 --- 向量的内积与正交(垂直),Orthogonal Vectors