为公司做网站,适合两个人运动前看的电影,wordpress修改字体插件,免费建立手机网站QR分解
给定一个mn的矩阵A#xff0c;其中m≥n#xff0c;即矩阵A是高矩阵或者是方阵#xff0c;QR分解将矩阵A分解为两个矩阵Q和R的乘积#xff0c;其中矩阵Q是一个mn的各列正交的矩阵#xff0c;即QTQI#xff0c;矩阵R是一个nn的上三角矩阵#xff0c;其对角线元素为…QR分解
给定一个m×n的矩阵A其中m≥n即矩阵A是高矩阵或者是方阵QR分解将矩阵A分解为两个矩阵Q和R的乘积其中矩阵Q是一个m×n的各列正交的矩阵即QTQI矩阵R是一个n×n的上三角矩阵其对角线元素为正。
如果矩阵A是方阵且各列线性无关那么Q是一个正交矩阵即QTQQQTI。
QR分解有多种算法实现包括Gram-Schmidt正交化方法、Householder变换方法和Givens旋转方法等下面我们介绍Gram-Schmidt正交化方法和Householder变换方法并在MATLAB平台上使用这两种算法来实现QR分解。
Gram-Schmidt算法
对于给定的n维向量a1a2……anGram-Schmidt算法可以解决将其标准正交化的问题即将一个线性无关的向量组转化为一个正交向量组使得每个向量都与前面的向量正交垂直并且可以检验a1a2……an是否是线性相关。
Gram-Schmidt算法的步骤如下
初始化n维向量q1q2……qn其中q1a1/||a1||2。对于每个向量aii2n进行正交化处理qi ai-( q1Tai)q1-…-( qi-1Tai)qi-1。如果qi0说明ai是a1a2……ai-1的一个线性组合可以结束算法了。否则将qi进行单位化qiqi/||qi||2。
如果步骤③没有结束那么说明a1a2……an是线性无关的而且得到了一个正交向量组q1q2……qn。
Gram-Schmidt算法实现的QR分解
对于给定矩阵A其列向量线性无关Gram-Schmidt算法实现的QR分解步骤如下
对列向量a1a2……an按照Gram-Schmidt方法进行正交化。对上一步得到的正交化向量组进行单位化得到各列正交的矩阵Q。根据AQRQTQI→RQTA得到上三角矩阵R
MATLAB验证Gram-Schmidt算法实现QR分解稳定性
通过直观的方法来观察到Gram-Schmidt QR分解的正交性偏差理论上通过Gram-Schmidt算法后可以得到列向量线性无关的各列正交的矩阵Q即QTQI我们可以直接计算QTQ看看计算结果与单位矩阵I的差距 左图是QTQ的计算结果有图是单位矩阵I可见由于浮点数存储的舍入误差随着k增大积累的误差越大矩阵Q逐渐失去正交性
clc,clear;
load MatrixA.mat;
[m,n]size(A);
Qzeros(m,n);
Rzeros(n,n);
%% Gram-Schmidt QR分解
for k1:nR(1:k-1,k)Q(:,1:k-1)*A(:,k); %求出R(1,K) - R(K-1,K)vA(:,k)-Q(:,1:k-1)*R(1:k-1,k); %求出正交化向量qR(k,k)norm(v); %求出RK,KQ(:,k)v/R(k,k); %单位化向量q
end
%% 正交性偏差
figure(1);
E zeros(1,n);
for k2:nmax 0;for i1:k-1temp abs(Q(:,i) * Q(:,k));if temp maxmax temp;endendE(1,k)max;
end
plot(E)
%% 比较QTQ和I
QTQQ*Q;
figure(2);
for i1:nfor j1:nscatter3(i,j,QTQ(i,j),red);hold on;end
end
zlim([0,1]);
Ieye(n);
figure(3);
for i1:nfor j1:nscatter3(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即wTv0那么有
Hv (I - 2wwT)v v - 2w(wTv) v
这说明对于垂直于w的向量householder变换的作用就是对其不起任何作用那么对于一个普通的向量v来说平行于w的分量被householder反向垂直于w的分量不变那么最终的效果就是将向量v作关于法向量为w的平面的镜像对称 基于Householder变换的QR分解 因为HH-1所以AH1H2,…,Hn-1R即Q H1H2,…,Hn-1再根据AQRQTQI→RQTA。
再来比较一下QTQ与单位矩阵I的差距结果如图所示左边的是计算出来的QTQ右边是单位矩阵I 结果QTQ和I基本一样可见相比其他分解方法Householder算法能够减小舍入误差的累积提高计算结果的稳定性。此外该算法的时间复杂度较低具备较高的计算效率。
clc,clear;
load MatrixA.mat;
[m,n]size(A);
Qzeros(m,n);
Rzeros(n,n);
%% Householder QR分解
[Q,R]qr(A); % matlab库函数就是用的Householder
%% 正交性偏差
figure(1);
E zeros(1,n);
for k2:nmax 0;for i1:k-1temp abs(Q(:,i) * Q(:,k));if temp maxmax temp;endendE(1,k)max;
end
plot(E)
%% 比较QTQ和I
QTQQ*Q;
figure(2);
for i1:nfor j1:nscatter3(i,j,QTQ(i,j),red);hold on;end
end
zlim([0,1]);
Ieye(n);
figure(3);
for i1:nfor j1:nscatter3(i,j,I(i,j),red);hold on;end
end
zlim([0,1]);
判断矩阵是否可逆
判断矩阵是否可逆有以下几种方法
存在一个矩阵B使得ABBAI确实可逆。矩阵行列式不为0可逆。矩阵满秩可逆。线性方程组Ax0只有0解可逆。线性方程组Axb只有特解可逆。
实际上如果一个方阵可以进行QR分解那么这个方阵也是可逆的。
所以我们直接尝试对矩阵B进行QR分解如果可以进行QR分解那么矩阵B可逆。那么我们可以先假设矩阵B是可以进行QR分解然后我们对矩阵B进行QR分解显然矩阵B是可以进行QR分解的这说明矩阵B是可逆的。
求逆
我们之前使用过高斯消元法来求解矩阵的逆实际上也可以使用QR分解求矩阵的逆。由A QRQTQ I则A-1 (QR)-1 R-1Q-1 R-1QT。
那么A-1就可以通过R-1QT得到但是实际上我们并不需要计算R-1让x R-1QT那么我们目标就是要得到x的结果因为RR-1QTQT即RxQT那么我们就需要求解这个线性方程组由于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);
Qzeros(m,n);
Rzeros(n,n);
%% Gram-Schmidt QR分解
for k1:nR(1:k-1,k)Q(:,1:k-1)*B(:,k); %求出R(1,K) - R(K-1,K)vB(:,k)-Q(:,1:k-1)*R(1:k-1,k); %求出正交化向量qR(k,k)norm(v); %求出RK,KQ(:,k)v/R(k,k); %单位化向量q
end
%% 求逆
inverseQRR\Q;
inverseinv(B);
%% 画图比较
for i0:n-1for j1:nscatter(i*nj,inverse(i1,j),red);hold on;scatter(i*nj,inverseQR(i1,j),green,.);hold on;end
end
我们再用之前的高斯消元法求解矩阵B的逆将其与用MATLAB内置的求逆函数结果进行比较结果如图所示 可见高斯消元法求逆的结果也很好基本上绿色的点都和红色的圆圈重合了。
clc,clear;
load MatrixB.mat;
beye(50);
B_b[B,b];
[n,m]size(B_b);
for i1:nfor jm:-1:iB_b(i,j)B_b(i,j)/B_b(i,i);endfor ji1:nfor km:-1:iB_b(j,k)B_b(j,k)-B_b(j,i)*B_b(i,k);endend
% fprintf(第%d次消元\n,i);
% disp(rats(A_b));
end
for in-1:-1:1for ji:-1:1for km:-1:n1B_b(j,k)B_b(j,k)-B_b(j,i1)*B_b(i1,k);endB_b(j,i1)0;end
% fprintf(第%d次回代\n,n-i);
% disp(rats(A_b));
end
gaussInverseB_b(:,end-49:end);
inverseinv(B);
%% 画图比较
for i0:n-1for j1:nscatter(i*nj,inverse(i1,j),red);hold on;scatter(i*nj,gaussInverse(i1,j),green,.);hold on;end
end
再用householder算法实现的QR分解求解矩阵B的逆将其与用MATLAB内置的求逆函数结果进行比较结果如图所示。 可见householder实现的QR分解求逆结果效果很好基本上和matlab内置求逆函数结果相同速度上也不慢。
clc,clear;
load MatrixB.mat;
[m,n]size(B);
Qzeros(m,n);
Rzeros(n,n);
%% Householder QR分解
[Q,R]qr(B); % matlab库函数就是用的Householder
%% 求逆
inverseQRR\Q;
inverseinv(B);
%% 画图比较
for i0:n-1for j1:nscatter(i*nj,inverse(i1,j),red);hold on;scatter(i*nj,inverseQR(i1,j),green,.);hold on;end
end
