矩阵的正交分解又称为QR分解,是将矩阵分解为一个正交矩阵Q和一个上三角矩阵的乘积的形式。
任意实数方阵A,都能被分解为 。这里的Q为正交单位阵,即 R是一个上三角矩阵。这种分解被称为QR分解。 QR分解也有若干种算法,常见的包括Gram–Schmidt、Householder和Givens算法。 QR分解是将矩阵分解为一个正交矩阵与上三角矩阵的乘积。用一张图可以形象地表示QR分解:
为啥我们需要正交分解呢? 实际运用过程中,QR分解经常被用来解线性最小二乘问题,这个问题我们后面讲述。
提到正交分解就不得不讨论(Householder transformation Householder变换)豪斯霍尔德变换和(Schmidt orthogonalization 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分解时,所论矩阵必须是列满秩矩阵。
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)==0
break;
end
Q(:,k)=A(:,k)/R(k,k);
for j=k+1:n
R(k,j)=Q(:,k)'*A(:,j);
A(:,j)=A(:,j)-R(k,j)*Q(:,k);
end
if nargin==2
b=Q'* b;
X(n)=b(n)/R(n,n);
for i=n-1:-1:1
X(i)=(b(i)-sum(R(i,i+1:n).*X(i+1:n)'))/R(i,i);
end
else
X=[];
end
end
设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是酉矩阵
也是酉矩阵;
于是有
如果α1=0,则直接进行下一步,此时相当于取
,而a1=0.
于是有
此时,令
则H2是n阶Householder矩阵,且使
如果α1=0,则直接进行下一步
令
,则Q是酉矩阵之积,从而必有酉矩阵并且A=QR
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==1
w=[A(1,1)+s,A(2:n,k)']';
else
w=[zeros(1,k-1),A(k,k)+s,A(k+1:n,k)']';
R(1:k-1,k)=A(1:k-1,k);
end
if norm(w)~=0
w=w/norm(w);
end
P=E-2*w*w';
A=P*A;
P1=P*P1;
R(1:n,n)=A(1:n,n);
end
Q=P1';
if nargin==2
b=P1*b;
X(n)=b(n)/R(n,n);
for i=n-1:-1:1
X(i)=(b(i)-sum(R(i,i+1:n).*X(i+1:n)'))/R(i,i);
end
else
X=[];
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等并行或者分布式的计算当中,加速我们的计算方法,这比直接单机去调用这些库会超乎想象的快。
(adsbygoogle = window.adsbygoogle || []).push({});