【matlab】QR分解

QR分解

给定一个m×n的矩阵A，其中m≥n，即矩阵A是高矩阵或者是方阵，QR分解将矩阵A分解为两个矩阵Q和R的乘积，其中矩阵Q是一个m×n的各列正交的矩阵，即QTQ=I，矩阵R是一个n×n的上三角矩阵，其对角线元素为正。

如果矩阵A是方阵，且各列线性无关，那么Q是一个正交矩阵，即QTQ=QQT=I。

QR分解有多种算法实现，包括Gram-Schmidt正交化方法、Householder变换方法和Givens旋转方法等，下面我们介绍Gram-Schmidt正交化方法和Householder变换方法，并在MATLAB平台上使用这两种算法来实现QR分解。

Gram-Schmidt算法

对于给定的n维向量a1，a2，……，an，Gram-Schmidt算法可以解决将其标准正交化的问题，即将一个线性无关的向量组转化为一个正交向量组，使得每个向量都与前面的向量正交（垂直），并且可以检验a1，a2，……，an是否是线性相关。

Gram-Schmidt算法的步骤如下：

  • 初始化n维向量q1，q2，……，qn，其中q1=a1/||a1||2。
  • 对于每个向量ai，i=2：n，进行正交化处理：qi= ai-( q1Tai)q1-…-( qi-1Tai)qi-1。
  • 如果qi=0，说明ai是a1，a2，……，ai-1的一个线性组合，可以结束算法了。
  • 否则将qi进行单位化，qi=qi/||qi||2。

如果步骤③没有结束，那么说明a1，a2，……，an是线性无关的，而且得到了一个正交向量组q1，q2，……，qn。

Gram-Schmidt算法实现的QR分解

对于给定矩阵A，其列向量线性无关，Gram-Schmidt算法实现的QR分解步骤如下：

  • 对列向量a1，a2，……，an按照Gram-Schmidt方法进行正交化。
  • 对上一步得到的正交化向量组进行单位化得到各列正交的矩阵Q。
  • 根据A=QR，QTQ=I→R=QTA，得到上三角矩阵R

MATLAB验证Gram-Schmidt算法实现QR分解稳定性

通过直观的方法来观察到Gram-Schmidt QR分解的正交性偏差，理论上通过Gram-Schmidt算法后可以得到列向量线性无关的各列正交的矩阵Q，即QTQ=I，我们可以直接计算QTQ，看看计算结果与单位矩阵I的差距

左图是QTQ的计算结果，有图是单位矩阵I，可见由于浮点数存储的舍入误差，随着k增大，积累的误差越大，矩阵Q逐渐失去正交性

clc,clear;
load MatrixA.mat;
[m,n]=size(A);
Q=zeros(m,n);
R=zeros(n,n);
%% Gram-Schmidt QR分解
for k=1:n
    R(1:k-1,k)=Q(:,1:k-1)'*A(:,k);  %求出R(1,K) - R(K-1,K)
    v=A(:,k)-Q(:,1:k-1)*R(1:k-1,k); %求出正交化向量q
    R(k,k)=norm(v);                 %求出R（K,K）
    Q(:,k)=v/R(k,k);                %单位化向量q
end
%% 正交性偏差
figure(1);
E = zeros(1,n);
for k=2:n
    max = 0;
    for i=1:k-1
        temp = abs(Q(:,i)' *  Q(:,k));
        if temp > max
            max = temp;
        end
    end
    E(1,k)=max;
end    
plot(E)
%% 比较QTQ和I
QTQ=Q'*Q;
figure(2);
for i=1:n
    for j=1:n
        scatter3(i,j,QTQ(i,j),'red');
        hold on;
    end
end
zlim([0,1]);
I=eye(n);
figure(3);
for i=1:n
    for j=1:n
        scatter3(i,j,I(i,j),'red');
        hold on;
    end
end
zlim([0,1]);

Householder变换

Householder变换是一种镜面反射变换，householder变换矩阵为H = I - 2wwT，如何理解这个变换矩阵呢，考虑向量w，那么有：

Hw = (I - 2wwT)w = w - 2w(wTw) = - w

这说明对于平行于w的向量，householder变换的作用是将其反向，再考虑与向量w垂直的向量v，即wTv=0，那么有：

Hv = (I - 2wwT)v = v - 2w(wTv) = v

这说明对于垂直于w的向量，householder变换的作用就是对其不起任何作用，那么对于一个普通的向量v来说，平行于w的分量被householder反向，垂直于w的分量不变，那么最终的效果就是将向量v作关于法向量为w的平面的镜像对称

 

基于Householder变换的QR分解

 因为H=H-1，所以A=H1H2,…,Hn-1R，即Q= H1H2,…,Hn-1，再根据A=QR，QTQ=I→R=QTA。

再来比较一下QTQ与单位矩阵I的差距，结果如图所示，左边的是计算出来的QTQ，右边是单位矩阵I

结果QTQ和I基本一样，可见相比其他分解方法，Householder算法能够减小舍入误差的累积，提高计算结果的稳定性。此外，该算法的时间复杂度较低，具备较高的计算效率。

clc,clear;
load MatrixA.mat;
[m,n]=size(A);
Q=zeros(m,n);
R=zeros(n,n);
%% Householder QR分解
[Q,R]=qr(A);    % matlab库函数就是用的Householder
%% 正交性偏差
figure(1);
E = zeros(1,n);
for k=2:n
    max = 0;
    for i=1:k-1
        temp = abs(Q(:,i)' *  Q(:,k));
        if temp > max
            max = temp;
        end
    end
    E(1,k)=max;
end    
plot(E)
%% 比较QTQ和I
QTQ=Q'*Q;
figure(2);
for i=1:n
    for j=1:n
        scatter3(i,j,QTQ(i,j),'red');
        hold on;
    end
end
zlim([0,1]);
I=eye(n);
figure(3);
for i=1:n
    for j=1:n
        scatter3(i,j,I(i,j),'red');
        hold on;
    end
end
zlim([0,1]);

判断矩阵是否可逆

判断矩阵是否可逆有以下几种方法：

  • 存在一个矩阵B，使得AB=BA=I，确实可逆。
  • 矩阵行列式不为0，可逆。
  • 矩阵满秩，可逆。
  • 线性方程组Ax=0只有0解，可逆。
  • 线性方程组Ax=b只有特解，可逆。

实际上如果一个方阵可以进行QR分解，那么这个方阵也是可逆的。

所以我们直接尝试对矩阵B进行QR分解，如果可以进行QR分解，那么矩阵B可逆。那么我们可以先假设矩阵B是可以进行QR分解，然后我们对矩阵B进行QR分解，显然矩阵B是可以进行QR分解的，这说明矩阵B是可逆的。

求逆

我们之前使用过高斯消元法来求解矩阵的逆，实际上也可以使用QR分解求矩阵的逆。由A = QR，QTQ = I，则A-1 = (QR)-1 = R-1Q-1 = R-1QT。

那么A-1就可以通过R-1QT得到，但是实际上我们并不需要计算R-1，让x= R-1QT，那么我们目标就是要得到x的结果，因为RR-1QT=QT，即Rx=QT，那么我们就需要求解这个线性方程组，由于R是上三角矩阵，所以直接回代就可以求出x，即求出R-1QT，即求出了A-1。

我们先用Gram-Schmidt算法实现的QR分解求解矩阵B的逆，将其与用MATLAB内置的求逆函数结果进行比较，结果如图所示，红色的圆圈是matlab内置的求逆函数计算出来的结果，绿色实心点是我们QR分解求出来的结果，如果二者重合说明计算结果相同。

可以看到基本上绿色的点都和红色的圆圈重合了，可见Gram-Schmidt算法QR分解求逆效果不错。

clc,clear;
load MatrixB.mat;
[m,n]=size(B);
Q=zeros(m,n);
R=zeros(n,n);
%% Gram-Schmidt QR分解
for k=1:n
    R(1:k-1,k)=Q(:,1:k-1)'*B(:,k);  %求出R(1,K) - R(K-1,K)
    v=B(:,k)-Q(:,1:k-1)*R(1:k-1,k); %求出正交化向量q
    R(k,k)=norm(v);                 %求出R（K,K）
    Q(:,k)=v/R(k,k);                %单位化向量q
end
%% 求逆
inverseQR=R\Q';
inverse=inv(B);
%% 画图比较
for i=0:n-1
    for j=1:n
        scatter(i*n+j,inverse(i+1,j),'red');
        hold on;
        scatter(i*n+j,inverseQR(i+1,j),'green','.');
        hold on;
    end
end

我们再用之前的高斯消元法求解矩阵B的逆，将其与用MATLAB内置的求逆函数结果进行比较，结果如图所示

可见高斯消元法求逆的结果也很好，基本上绿色的点都和红色的圆圈重合了。

clc,clear;
load MatrixB.mat;
b=eye(50);
B_b=[B,b];
[n,m]=size(B_b);
for i=1:n
    for j=m:-1:i
        B_b(i,j)=B_b(i,j)/B_b(i,i);
    end
    for j=i+1:n
        for k=m:-1:i
            B_b(j,k)=B_b(j,k)-B_b(j,i)*B_b(i,k);
        end
    end
%     fprintf('第%d次消元\n',i);
%     disp(rats(A_b));
end
for i=n-1:-1:1
    for j=i:-1:1
        for k=m:-1:n+1
            B_b(j,k)=B_b(j,k)-B_b(j,i+1)*B_b(i+1,k);
        end
        B_b(j,i+1)=0;
    end
%     fprintf('第%d次回代\n',n-i);
%     disp(rats(A_b));
end
gaussInverse=B_b(:,end-49:end);
inverse=inv(B);
%% 画图比较
for i=0:n-1
    for j=1:n
        scatter(i*n+j,inverse(i+1,j),'red');
        hold on;
        scatter(i*n+j,gaussInverse(i+1,j),'green','.');
        hold on;
    end
end

再用householder算法实现的QR分解求解矩阵B的逆，将其与用MATLAB内置的求逆函数结果进行比较，结果如图所示。

可见householder实现的QR分解求逆结果效果很好，基本上和matlab内置求逆函数结果相同，速度上也不慢。

clc,clear;
load MatrixB.mat;
[m,n]=size(B);
Q=zeros(m,n);
R=zeros(n,n);
%% Householder QR分解
[Q,R]=qr(B);    % matlab库函数就是用的Householder
%% 求逆
inverseQR=R\Q';
inverse=inv(B);
%% 画图比较
for i=0:n-1
    for j=1:n
        scatter(i*n+j,inverse(i+1,j),'red');
        hold on;
        scatter(i*n+j,inverseQR(i+1,j),'green','.');
        hold on;
    end
end

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：/a/217028.html

如若内容造成侵权/违法违规/事实不符，请联系我们进行投诉反馈qq邮箱809451989@qq.com，一经查实，立即删除！

相关文章

初识动态规划算法(题目加解析)

初识动态规划算法(题目加解析)

文章目录 什么是动态规划正文力扣题第 N 个泰波那契数三步问题使用最小花费爬楼梯 总结 什么是动态规划 线性动态规划：是可以用一个dp表来存储内容，并且找到规律存储,按照规律存储。让第i个位置的值等于题目要求的答案 >dp表：dp表就是用一…
阅读更多...
【数据结构】——栈|队列（基本功能）

【数据结构】——栈|队列（基本功能）

目录 栈 基本概念 栈的常见基本操作 栈的存储 ✌栈的基本操作实现 栈的构建 栈的初始化 入栈 打印栈 出栈 获取栈顶元素 获取栈的有效元素个数 判断栈是否为空 销毁栈 队列 基本概念 队列的常见基本操作 ✌队列的基本操作实现 队列的构建 初始化 入队列 出…
阅读更多...
不再只是android，华为自爆Harmony将对标iOS

不再只是android，华为自爆Harmony将对标iOS

今年10月，华为官方宣布，鸿蒙OS 4升级设备数量已突破1亿，成为史上升级最快的鸿蒙OS版本。 日前，据数码博主“定焦数码”消息，大厂技术员工做适配，通过线下沟通时，华为反复提到一个问题&#xff…
阅读更多...
实战技巧：为Android应用设置独立的多语言

实战技巧：为Android应用设置独立的多语言

原文链接 实战技巧：为Android应用设置独立的多语言 通常情况下多语言的设置都在系统设置中，应用需要做的就是提供本应用所使用的字串的多语言翻译，使用时使用R.string.app_name类似的引用，然后系统会根据用户在系统设置中的选项来…
阅读更多...
不瞒各位，不安装软件也能操作Xmind文档

不瞒各位，不安装软件也能操作Xmind文档

大家好，我是小悟 作为搞技术的一个人群，时不时就要接收产品经理发过来的思维脑图，而此类文档往往是以Xmind编写的，如果你的电脑里面没有安装Xmind的话，不好意思，是打不开这类后缀结尾的文档。 打不开的话…
阅读更多...
【雷电模拟器桥接问题解决方法】

【雷电模拟器桥接问题解决方法】

1.ROOT权限开启 2.开启网络桥接模式，选择静态IP设置，点击安装桥接网卡，填写IP地址（注意：IP地址要与host主机在同一IP段内） 3.重启后 adb shell就能进入到模拟器控制台中了，如果出现以下内容&…
阅读更多...
进程程序替换和shell实现

进程程序替换和shell实现

先前fork说创建子进程执行代码，如何让子进程执行和父进程完全不一样的代码?程序替换。 一 单进程替换演示 1 execl函数使用 最近转到在vs code下写代码，之前也在xhell下用过execl函数，所以才想写篇博客总结总结，没想到在vs code…
阅读更多...
（C语言）计算n的阶乘

（C语言）计算n的阶乘

要求使用双精度 #include<stdio.h> double factorial(int n) {if(n 1)return 1;return n * factorial(n-1); } int main() {int n ;double res;scanf("%d",&n);res factorial(n);printf("%lf",res); return 0; } 运行截图&#xff1a; 注&am…
阅读更多...
oops-framework框架 之 界面管理(三)

oops-framework框架 之 界面管理(三)

引擎&#xff1a; CocosCreator 3.8.0 环境&#xff1a; Mac Gitee: oops-game-kit 注&#xff1a; 作者dgflash的oops-framework框架QQ群&#xff1a; 628575875 回顾 在上文中主要通过oops-game-kit大家了一个新的模版项目&#xff0c; 主要注意项是resources目录下的两个文…
阅读更多...
Python Opencv实践 - Yolov3目标检测

Python Opencv实践 - Yolov3目标检测

本文使用CPU来做运算&#xff0c;未使用GPU。练习项目&#xff0c;参考了网上部分资料。 如果要用TensorFlow做检测&#xff0c;可以参考这里 使用GPU运行基于pytorch的yolov3代码的准备工作_little han的博客-CSDN博客文章浏览阅读943次。记录一下自己刚拿到带独显的电脑&a…
阅读更多...
卷积神经网络（CNN）：艺术作品识别

卷积神经网络（CNN）：艺术作品识别

文章目录 一、前言一、设置GPU二、导入数据1. 导入数据2. 检查数据3. 配置数据集4. 数据可视化 三、构建模型四、编译五、训练模型六、评估模型1. Accuracy与Loss图2. 混淆矩阵3. 各项指标评估 一、前言 我的环境&#xff1a; 语言环境&#xff1a;Python3.6.5编译器&#xf…
阅读更多...
继承 多态 拆箱装箱 128陷阱 枚举类

继承 多态 拆箱装箱 128陷阱 枚举类

继承 在java里一个类只能继承一个类&#xff0c;但可以被多个类继承&#xff1b;c里一个类可以继承多个类&#xff1b; 子类可以使用父类的方法&#xff1b; 在java中&#xff0c;Object是所有类的父类&#xff1b; equals方法比较的是对象是否指向同一个地方&#xff0c;这个方…
阅读更多...
原生横向滚动条 吸附 页面底部

原生横向滚动条 吸附 页面底部

效果图 /** 横向滚动条 吸附 页面底部 */ export class StickyHorizontalScrollBar {constructor(options {}) {const { el, style } optionsthis.createScrollbar(style)this.insertScrollbar(el)this.setScrollbarSize()this.onEvent()}/** 创建滚轴组件元素 */createS…
阅读更多...
Windows下打包C++程序无法执行：无法定位程序输入点于动态链接库

Windows下打包C++程序无法执行：无法定位程序输入点于动态链接库

1、问题描述 环境&#xff1a;CLionCMakeMinGW64遇到问题&#xff1a;打包的exe无法运行&#xff0c;提示无法定位程序输入点于动态链接库。 2、解决思路 ​ 通过注释头文件的方式&#xff0c;初步定位问题是因为使用了#include <thread> 多线程库引起的。而且exe文件…
阅读更多...
外包干了2个月，技术倒退2年。。。。。

外包干了2个月，技术倒退2年。。。。。

先说一下自己的情况&#xff0c;本科生&#xff0c;20年通过校招进入深圳某软件公司&#xff0c;干了接近4年的功能测试&#xff0c;今年国庆&#xff0c;感觉自己不能够在这样下去了&#xff0c;长时间呆在一个舒适的环境会让一个人堕落!而我已经在一个企业干了四年的功能测试…
阅读更多...
如何创建maven项目的多模块项目

如何创建maven项目的多模块项目

Maven多模块项目是指一个Maven项目中包含多个子模块&#xff0c;每个子模块又是一个独立的Maven项目&#xff0c;但它们之间可以存在依赖关系。Maven多模块项目可以方便地管理多个子模块的依赖和构建过程&#xff0c;同时也可以提高项目的可维护性和可扩展性。创建maven项目的父…
阅读更多...
ChatGPT发布一年后，搜索引擎的日子还好吗？

ChatGPT发布一年后，搜索引擎的日子还好吗？

导读&#xff1a;生成式AI&#xff0c;搜索引擎的终结者还是进化加速器 ChatGPT发布刚刚一年&#xff0c;互联网世界已经换了人间。 2023年&#xff0c;以ChatGPT和大模型为代表的生成式AI浪潮对全球互联网、云计算、人工智能领域都带来巨大冲击。而且生成式AI在各行各业的应用…
阅读更多...
深入理解JVM虚拟机第二十七篇：详解JVM当中InvokeDynamic字节码指令，Java是动态类型语言么？

深入理解JVM虚拟机第二十七篇：详解JVM当中InvokeDynamic字节码指令，Java是动态类型语言么？

😉😉 学习交流群： ✅✅1：这是孙哥suns给大家的福利！ ✨✨2：我们免费分享Netty、Dubbo、k8s、Mybatis、Spring...应用和源码级别的视频资料 🥭🥭3：QQ群：583783824 📚📚 工作微信：BigTreeJava 拉你进微信群，免费领取！ 🍎🍎4：本文章内容出自上述：Sp…
阅读更多...
[ROS2] --- ROS diff ROS2

[ROS2] --- ROS diff ROS2

1 ROS存在的问题 一旦Ros Master主节点挂掉后&#xff0c;就会造成整个系统通信的异常,通信基于TCP实现&#xff0c;实时性差、系统开销大对Python3支持不友好&#xff0c;需要重新编译消息机制不兼容没有加密机制、安全性不高 2 ROS and ROS2架构对比 ROS和ROS2架构如下图所…
阅读更多...
Redis实战篇笔记（最终篇）

Redis实战篇笔记（最终篇）

Redis实战篇笔记&#xff08;七&#xff09; 文章目录 Redis实战篇笔记&#xff08;七&#xff09;前言达人探店发布和查看探店笔记点赞点赞排行榜 好友关注关注和取关共同关注关注推送关注推荐的实现 总结 前言 本系列文章是Redis实战篇笔记的最后一篇&#xff0c;那么到这里…
阅读更多...
最新文章

Copyright @ 2022~2023