计算方法(三)矩阵分解1-正交分解(QR分解)
为什么80%的码农都做不了架构师?>>>
正交分解
矩阵的正交分解又称为QR分解,是将矩阵分解为一个正交矩阵Q和一个上三角矩阵的乘积的形式。
任意实数方阵A,都能被分解为 。这里的Q为正交单位阵,即 R是一个上三角矩阵。这种分解被称为QR分解。 QR分解也有若干种算法,常见的包括Gram–Schmidt、Householder和Givens算法。 QR分解是将矩阵分解为一个正交矩阵与上三角矩阵的乘积。用一张图可以形象地表示QR分解:
为啥我们需要正交分解呢?
实际运用过程中,QR分解经常被用来解线性最小二乘问题,这个问题我们后面讲述。
提到正交分解就不得不讨论(Householder transformation Householder变换)豪斯霍尔德变换和(Schmidt orthogonalization Schmidt正交化)施密特正交化
Schmidt正交化
定理1 设A是n阶实非奇异矩阵,则存在正交矩阵Q和实非奇异上三角矩阵R使A有QR分解;且除去相差一个对角元素的绝对值(模)全等于1的对角矩阵因子外,分解是唯一的.
定理2 设A是m×n实矩阵,且其n个列向量线性无关,则A有分解A=QR,其中Q是m×n实矩阵,且满足QHTQ=E,R是n阶实非奇异上三角矩阵该分解除去相差一个对角元素的绝对值(模)全等于1的对角矩阵因子外是唯一的.用Schmidt正交化分解方法对矩阵进行QR分解时,所论矩阵必须是列满秩矩阵。
算法步骤
- 写出矩阵的列向量;
- 列向量按照Schmidt正交化正交;
- 得出矩阵的Q′,R′;
- 对R′的列向量单位化得到Q,R′的每行乘R′每列的模得푹
matlab代码
function[X,Q,R] = QRSchmidt(A,b)
%方阵的QR的Gram-Schmidt正交化分解法,并用于求解AX=b方程组[m,n]=size(A);
if m~=n%如果不是方阵,则不满足QR分解要求disp('不满足QR分解要求');
end
Q=zeros(m,n);
X=zeros(n,1);
R=zeros(n);
for k=1:nR(k,k)=norm(A(:,k));if R(k,k)==0break;endQ(:,k)=A(:,k)/R(k,k);for j=k+1:nR(k,j)=Q(:,k)'*A(:,j); A(:,j)=A(:,j)-R(k,j)*Q(:,k);end
if nargin==2b=Q'* b;X(n)=b(n)/R(n,n);for i=n-1:-1:1X(i)=(b(i)-sum(R(i,i+1:n).*X(i+1:n)'))/R(i,i);end
elseX=[];
end
end
Householder变换
设A为任一n阶方阵,则必存在n阶酉矩阵Q和n阶上三角阵R,使得A=QR
设w∈Cn是一个单位向量,令
则称H是一个Householder矩阵或Householder变换。则对于任意的 存在Householder矩阵H,使得Hx=-au。其中
酉矩阵(unitary matrix)
若n阶复矩阵A满足
则称A为酉矩阵,记之为其中,Ah是A的共轭转置
酉矩阵性质
如果A是酉矩阵
- 也是酉矩阵;
- det(A)=1;
- 充分条件是它的n个列向量是两两正交的单位向量。
算法步骤
- 将矩阵A按列分块写成A=(α1,α2,...,αn).如果α1≠0,则可得,存在n阶householder矩阵H1使得于是有
如果α1=0,则直接进行下一步,此时相当于取,而a1=0.- 将矩阵An-1按列分块写成An-1=(αi,α2,... ,αn-1)。如果α1≠0,则可得,存在n-1阶householder矩阵H’2使得 于是有
此时,令
则H2是n阶Householder矩阵,且使
如果α1=0,则直接进行下一步- 对n-2阶矩阵继续进行类似的变换,如此下去,之多在第n-1步,我们可以找到Householder矩阵H1,H2,...,Hn-1使得
令,则Q是酉矩阵之积,从而必有酉矩阵并且A=QR
matlab代码
function[ X,Q,R ] = QRHouseholder(A,b)
%用Householder变换将方阵A分解为正交Q与上三角矩阵R的乘积,并用于求解AX=b方程组
[n,n]=size(A);
E=eye(n);
X=zeros(n,1);
R=zeros(n);
P1=E;
for k=1:n-1%构造w,使Pk=I-2ww's=-sign(A(k,k))* norm(A(k:n,k));R(k,k)=-s;if k==1w=[A(1,1)+s,A(2:n,k)']';elsew=[zeros(1,k-1),A(k,k)+s,A(k+1:n,k)']';R(1:k-1,k)=A(1:k-1,k);endif norm(w)~=0w=w/norm(w);endP=E-2*w*w';A=P*A;P1=P*P1;R(1:n,n)=A(1:n,n);
end
Q=P1';
if nargin==2b=P1*b;X(n)=b(n)/R(n,n);for i=n-1:-1:1X(i)=(b(i)-sum(R(i,i+1:n).*X(i+1:n)'))/R(i,i);end
elseX=[];
end
matlab自带方法
%产生一个3*3大小的魔方矩阵
A=magic(3)
[Q,R]=qr(A)
使用Eigen C++ Eigen提供了几种矩阵分解的方法
分解方式 | Method | 矩阵满足条件 | 计算速度 | 计算精度 |
---|---|---|---|---|
PartialPivLU | partialPivLu() | Invertible | ++ | + |
FullPivLU | fullPivLu() | None | - | +++ |
HouseholderQR | householderQr() | None | ++ | + |
ColPivHouseholderQR | colPivHouseholderQr() | None | + | ++ |
FullPivHouseholderQR | fullPivHouseholderQr() | None | - | +++ |
LLT | llt() | Positive definite | +++ | + |
LDLT | ldlt() | Positive or negative semidefinite | +++ | ++ |
其中HouseholderQR、ColPivHouseholderQR、FullPivHouseholderQR是我们目前要用到的QR分解方法
C++的QR分解代码为
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
int main() {Matrix3d A;A<<1,1,1,2,-1,-1,2,-4,5;HouseholderQR<Matrix3d> qr;qr.compute(A);MatrixXd R = qr.matrixQR().triangularView<Upper>();MatrixXd Q = qr.householderQ();std::cout << "QR2(): HouseholderQR---------------------------------------------"<< std::endl;std::cout << "A "<< std::endl <<A << std::endl << std::endl;std::cout <<"qr.matrixQR()"<< std::endl << qr.matrixQR() << std::endl << std::endl;std::cout << "R"<< std::endl <<R << std::endl << std::endl;std::cout << "Q "<< std::endl <<Q << std::endl << std::endl;std::cout <<"Q*R" << std::endl <<Q*R << std::endl << std::endl;return 0;
}
输出
好了大功告成,为什么我要写计算方法的文章呢,虽然现在有很多的库和包给我们调用,但是我们也不能忘了代码的本质是为了解决复杂的数学问题,从根源上去理解一种计算方法有助于我们对自身代码的优化,比如这些方法我们可以把它写到FPGA和CUDA等并行或者分布式的计算当中,加速我们的计算方法,这比直接单机去调用这些库会超乎想象的快。
转载于:https://my.oschina.net/VenusV/blog/3030201
计算方法(三)矩阵分解1-正交分解(QR分解)相关推荐
- qr分解实验 matlab,QR分解与最小二乘
主要内容: 1.QR分解定义 2.QR分解求法 3.QR分解与最小二乘 4.Matlab实现 一.QR分解 R分解法是三种将矩阵分解的方式之一.这种方式,把矩阵分解成一个正交矩阵与一个上三角矩阵的积. ...
- 豪斯霍尔德qr分解java_[转]QR分解和酉矩阵
预备知识 平面旋转与 Householder 矩阵是特殊的酉矩阵,它们在建立某些基本的矩阵分解过程中起着重要的作用. 平面旋转 设 1⩽i 为平面旋转或者 Givens 旋转. 容易验证对任何一对指数 ...
- matlab qr分解的作用,QR分解在MatLab
我的任務是擬合一個多項式的數據.我想使用Gram-Schimdt正交化過程來實現QR算法.它是建立在這樣的功能: function [ Q,R ] = QRDec(A) n = length(A(1, ...
- java 矩阵分解_计算方法(三)矩阵分解1-正交分解(QR分解)
正交分解 矩阵的正交分解又称为QR分解,是将矩阵分解为一个正交矩阵Q和一个上三角矩阵的乘积的形式. 任意实数方阵A,都能被分解为 .这里的Q为正交单位阵,即 R是一个上三角矩阵.这种分解被称为QR分解 ...
- QR算法的Matlab 程序,三种实现矩阵QR分解的算法与程序
To learn, to share, to debate, then comes progress. ------------------------------------------------ ...
- [矩阵的QR分解系列四] QR(正交三角)分解
QR分解 简介 QR分解 定义 存在和唯一性 存在性证明 唯一性证明 分解方法 施密特(Schmidt)方法 吉文斯(Givens)方法 豪斯霍尔德(Householder)方法 例子 施密特(Sch ...
- crout分解计算例题_矩阵与数值计算(2)——矩阵三角分解LU、PALU、Cholesky三角分解、QR分解...
前言 矩阵分解是设计算法的主要技巧,通过分解可以将复杂问题转化为几个简单问题求解,通常完成这一转化任务的主要技巧就是矩阵分解.例如,我们知道上三角矩阵和下三角矩阵是容易求解的,或者对角矩阵是最理想的求 ...
- 矩阵的QR分解以及在最小二乘法中的应用
一.最小二乘法 最小二乘法是一种数学优化方法,通过最小化误差的平方和来拟合数据点. 以线性回归模型为例,如果我们用最小二乘法来求解线性回归的系数,可得: err(yi−y^)=1n∑i=1n( ...
- [矩阵的QR分解系列一] 施密特(Schmidt)正交规范化
施密特正交规范化 简介 规范化步骤 例子 引用 之前介绍的矩阵的三角分解系列介绍了利用矩阵初等变换解决了矩阵三角化问题以及具体的三角分解.但是以初等变换工具的三角分解方法并不能消除病态线性方程组不稳定 ...
最新文章
- PyTorch分布式训练
- 手写 单隐藏层神经网络_反向传播(Matlab实现)
- 解决 Microsoft Excel has stopped working
- java 获取服务器硬件_dell服务器远程获取硬件状态
- X-Mas Musings –在Grails集成测试中不要使用随机服务器端口
- 帖子回复——无限级分类
- 数据库---三大设计范式
- 那些年开发中遇到的坑。。。
- 经典排序算法(十八)--Proxmap Sort
- 小程序中间放大轮播图_微信小程序实现类3D轮播图
- ZEMAX实例学习2:双透镜(a doublet)
- c语言编程中分数怎么表示,用C语言编程平均分数
- 百度网盘获取下载链接
- ⭐App爬虫系列⭐:获取王者荣耀全英雄的名称、类型、热度、胜率、登场率、Ban率
- 常见的测试用例设计方法8---正交试验法
- 理解标准差、标准化、协方差、正态分布
- 5个最优秀的Java和C#代码转换工具
- 现身说法:37岁老码农找工作!
- 为什么你还没有买新能源汽车?
- JPEG压缩算法详解(转载)
热门文章
- 论文阅读【EMScore: Evaluating Video Captioning via Coarse-Grained and Fine-Grained Embedding Matching】
- 《公民的不服从》---梭罗(1) 英文翻译3
- 一目了然-火焰图初探
- MFC CListCtrl右键菜单
- 什么是前端开发及学习路线
- 双非计算机考研可以调剂的学校,考研调剂时别盲目选择名校!最好的调剂院校其实是你的本科院校!...
- oracle去重复合计,如何在excel合并同类项数据并求和(去除重复项)
- LINUX安装nodejs(centos7)
- Fiddler2 下断点修改HTTP报文
- spring 容器中bean的扩展点记录 —— 个人学习记录