数值代数（一）QR分解

Laplaces 1988

已于 2023-07-01 20:27:25 修改

阅读量8.8k

点赞数 8

分类专栏：笔记文章标签： matlab

于 2022-01-20 19:34:17 首次发布

本文链接：https://blog.csdn.net/weixin_45612813/article/details/122607837

版权

笔记专栏收录该内容

2 篇文章 0 订阅

订阅专栏

这学期刚学完数值代数，在这里总结一下学过的算法，附MATLAB代码。

QR分解有两种形式：完全QR分解（reduced QR factorization）、约化QR分解（full QR factorization）。

完全QR分解： $A=QR$ ，其中 $A\in \mathbb{C}^{m\times n}(m\geqslant n)$ ， $Q\in \mathbb{C}^{m\times m}$ 为酉矩阵， $R\in \mathbb{C}^{m\times n}$ 为上三角矩阵。

约化QR分解： $A=\hat{Q}\hat{R}$ ，其中 $A\in \mathbb{C}^{m\times n}(m\geqslant n)$ ， $\hat{Q}\in \mathbb{C}^{m\times n}$ 具有正交列， $\hat{R}\in \mathbb{C}^{n\times n}$ 为上三角矩阵。

实际上，两者是等价的。因为 $Q$ 比 $\hat{Q}$ 多的几列可以通过扩充标准正交基得到， $R$ 比 $\hat{R}$ 多的几行全为0。我认为完全QR分解最大的优点就是 $Q$ 是酉矩阵，很好处理。

以下是三中QR分解的经典算法，其中豪斯霍尔德算法工作量最少且向后稳定。

1.经典格拉姆-施密特正交化（classical Gram-Schmidt iteration）

基本思想：

找到单位正交序列 $q_{1},q_{2},...$ 满足 $<q_{1},...,q_{j}>\supseteq <a_{1},...,a_{j}>,j=1,...,n$

这等价于

利用这种思想对矩阵 $A$ 进行约化QR分解。

代码实现

function [Q,R]=classical_gramschmidt(A)
[~,n] = size(A);
V = zeros(size(A));
Q = zeros(size(A));
R = zeros(n);
for j = 1:n
    V(:,j) = A(:,j);
    for i = 1:j-1
        R(i,j) = Q(:,i)'*A(:,j);
        V(:,j) = V(:,j)-R(i,j)*Q(:,i);
    end
    R(j,j) = norm(V(:,j),2);
    Q(:,j) = V(:,j)/R(j,j);
end
end

2.修正的格拉姆-施密特正交化（modified Gram-Schmidt iteration）

基本思想

设 $P_{j}$ 为投影到正交于 $range(<q_{1},...,q_{j-1}>)$ 的正交投影算子，这样做可以保证 $q_{j}\perp <q_{1},...,q_{j-1}>$ ，

由正交投影算子的相关知识可以得到， $P_{j}=I-\hat{Q}_{j-1}\hat{Q}_{j-1}^{*}$ ，其中

代码实现

function [Q,R]=modified_gramschmidt(A)
[~,n] = size(A);
Q = zeros(size(A));
R = zeros(n);
V=A;
for i = 1:n
    R(i,i) = norm(V(:,i),2);
    Q(:,i) = V(:,i)/R(i,i);
    for j = i+1:n
        R(i,j) = Q(:,i)'*V(:,j);
        V(:,j) = V(:,j)-R(i,j)*Q(:,i);
    end
end
end

3.豪斯霍尔德三角形化（Householder）

function [Q,A] = householder(A)
[m,n] = size(A);
Q = eye(m);
for k = 1:n
    x = A(k:m,k);
    v = sign(x(1))*norm(x,2)*eye(m-k+1,1)+x;
    v = v/norm(v,2);
    F = eye(m-k+1)-2*v*v'/(v'*v);
    Q = [eye(k-1) zeros(k-1,m-k+1);zeros(m-k+1,k-1) F]*Q; 
    A(k:m,k:n) = A(k:m,k:n)-2*v*(v'*A(k:m,k:n));
end
Q = Q';
end