QR分解的四种常见方法

饕餮之徒和頑丫头

于 2023-10-10 16:09:30 发布

阅读量1.3k

点赞数

文章标签：算法

本文链接：https://blog.csdn.net/m0_74981381/article/details/133613808

版权

一、gram-schmidt正交法之CGS正交法

1）相较于MGS的特点：全部一起消从而得到 $\alpha _{i}$

2）原理解释：给定矩阵 $A_{mn}$ ，将其看成由n个列向量组成的向量矩阵 $A_{mn}=\left \{ \alpha_{1} ,\alpha _{2,...,}\alpha _{n}\right \}$ ，分别给 $\alpha _{1},\alpha _{2},...,\alpha _{n}$ 进行正交化、单位化

正交化：

$\alpha _{1}=\alpha _{1}$

$\alpha _{2}=\alpha _{2}-\frac{\alpha _{2}'\alpha _{1}}{\alpha _{1}'\alpha _{1}}\alpha _{1},...$

单位化： $\alpha_{1}=\alpha _{1}/norm(\alpha _{1}),\alpha _{2}=\alpha _{2}/norm(\alpha _{2}),...$

3）代码示例

二、gram-schmidt正交法之MGS正交法

1）MGS(modified gram-schmidt)相比于CGS的特点：一个一个消得 $\alpha _{i}$

2）原理解释：

上帝给我们一个需要单位正交化的矩阵Anxn=(a1,a2,...,an)和标准正交基q1、q2、q3...qn,假设dimension(A)=3,则有以下为步骤：

第一步，a1=a1

第二步，

$\alpha$ 12=q1*a2
a2=a2-q1* $\alpha$ 12
$\alpha$ 22=q2*a2
a2=a2-q2* $\alpha$ 22

第三步

$\alpha$ 13=q1'*a3
a3=a3-q1* $\alpha$ 13
$\alpha$ 23=q2'*a3
a3=a3-q2* $\alpha$ 23
$\alpha$ 33=q3'*a3
a3=a3-q3* $\alpha$ 33

第四步 ai=ai/norm(ai),i=1、2、3

所以得到单位正交阵（a1,a2,a3),上面的标准正交基我们在之前的学习中遇到的是（e1，e2，e3，...，en）=diag(n)

3）代码示例：

n=15;
A=rand(n);
st_b=diag(n);%standard_base
alpha=zeros(n);
for j=2:n    
    for i=1:j
        alpha(i,j)=st_b(i)'*A(:,j);
        A(:,j)=A(:,j)-q1*alpha(i,j);
    end
end
for i=1:n
    A(:,i)=A(:,i)/norm(A(:,i));
end

三、householder 镜像法householder reflections

原理说明：

对于一个基 $w$ (注意 $nrom(w)=1$ ),若想要得到 $x$ 关于该基 $w$ 的超平面镜像映射，（将矩阵反射到以标准基为基的向量空间里面）可使用 $H=I-2ww^{T}$ ，使得 $x=\left ( x_{1},x_{2},x_{3},...,x_{n} \right )$ 变为 $x=\left ( x_{1}-nrom(x_{}),0,0,...,0 \right )$

在上机时，为了简化，用 $v/norm(v)$ （其中 $v=x-norm(x)e_{1}$ ）代替 $w$ ，得 $2ww^{T}=2vv^{T}/(v^{T}v)$

代码：在下列代码中，我用 $wvv^{T}$ （其中 $w=2/(v^{T}v)$ ）代替 $2vv^{T}/(v^{T}v)$

四、Givens旋转法Givens rotations

原理说明：现在有两个向量x1,x2,我想要将他们旋转到标准基的向量空间里，则使用矩阵G即可

$G=\begin{pmatrix} cos(\theta ) & sin(\theta )\\ -sin(\theta ) & cos(\theta ) \end{pmatrix}$ , $cos(\theta )=x_{1}/\sqrt{x_{1}^{2}+x_{2}^{2}},sin(\theta )=x_{2}/\sqrt{x_{1}^{2}+x_{2}^{2}}$

$G\begin{pmatrix} x_{1}\\ x_{2}\end{pmatrix}=\begin{pmatrix} x_{1}^{hat}\\ 0\end{pmatrix}$

代码：

n=10;
m=n;
A=orth(rand(m,n));
A_origional=A;
C_s=zeros(m,n);
for j=1:n-1
    for i=m:-1:(j+1)
        if abs(A(i,j))>=abs(A(j,j))
           t=A(j,j)/A(i,j);
           sign=A(i,j)/norm(A(i,j));
           s=sign*1/sqrt(1+t^2);
           c=t*s;
        else
            t=A(i,j)/A(j,j);
            sign=A(j,j)/norm(A(j,j));
            c=sign*1/sqrt(1+t^2);
            s=t*c;
        end
        A(j,:)=c*A(j,:)+s*A(i,:);
        A(i,:)=-s*A(j,:)+c*A(i,:);
        C_s(i,j)=c;
        if s>0 
            C_s(m-i+1,n-j+1)=1;
        else
            C_s(m-i+1,n-j+1)=-1;
        end
    end
end
Q=A_origional/A;//这是使用Q=A*R^(-1)
G=eye(m,n);
E=eye(m,n);
ei=zeros(n,1);
ej=zeros(n,1);
for j=1:(n-1)//此循环使用了Q=G1*G2*G3...*G((n-1)*(n+2)/2)
    ej(j)=1;
    for i=m:-1:(j+1)
        s=sqrt(1-C_s(i,j)^2)*C_s(m-i+1,n-j+1);
        ei(i)=1;
        G0=E+s*(ei*ej'-ej*ei')+(c-1)*(ei*ei'+ej*ej');
        G=G*G0;
        ei(i)=0;
    end
    ej(j)=0;
end

我有个问题想请教正在看这篇文章的你：我在使用GIVENS求QR分解时，Q=A/R求得的Q不等于Q=G1*G2*G3...Q((n-1)*(n+2)/2)得到的Q,是我无意中放大误差了还是代码出现错误？