Math-Model（五）正交分解(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

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

%产生一个3*3大小的魔方矩阵
A=magic(3)
[Q,R]=qr(A)

#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;
}