【矩阵的约化】用matlab实现基于Givens变换与Householder变换的QR分解

本文深入探讨了Givens变换与Householder变换在QR分解中的具体应用及其数学原理，通过MATLAB编程实现了这两种变换。本文还提供了直观的代码示例和对比分析，旨在帮助读者更好地理解概念及计算过程。

m0_74150590

2152人浏览 · 2024-12-14 21:20:04

m0_74150590 · 2024-12-14 21:20:04 发布

善用目录，可以根据需求直接跳转至代码部分~

欢迎讨论，如有错误麻烦指正(*^▽^*)

本篇不进行详细解释，深入学习可以参考教材

在理解算式之前，先明确目标: 得到一个上三角矩阵方便线性方程组的计算，由此衍生出QR分解。

一、Givens变换

借用了解析几何里面的旋转矩阵：

$T=\begin{pmatrix} cos\theta & sin\theta\\ -sin\theta & cos\theta \end{pmatrix}$

对二维向量x, y，y=Tx. 即x顺时针旋转θ角到y.

我们的目标是得到一个上三角矩阵，那我们希望，矩阵X的经过变化得到的矩阵Y在对角线下放的元素为0。根据旋转矩阵，T一次只能变化矩阵X一列的元素值的一个值，那么我们就把目标拆解为从第一列对角线以下的非零元素开始，一个一个地，逐步将矩阵下方的元素化为0。

如何利用旋转矩阵将元素化为0呢，依旧以2维向量为例，设 $x=(x_1,x_2)^T$ ,

$T=\begin{pmatrix} cos\theta & sin\theta\\ -sin\theta & cos\theta \end{pmatrix}=\begin{pmatrix} \frac{x_1}{\sqrt{x^2_1+x^2_1}} & \frac{x_2}{\sqrt{x^2_1+x^2_1}} \\ -\frac{x_2}{\sqrt{x^2_1+x^2_1}} & \frac{x_1}{\sqrt{x^2_1+x^2_1}} \end{pmatrix}$

$y=Tx=\begin{pmatrix} \sqrt{x^2_1+x^2_2}\\ 0\end{pmatrix}$

当我们简记cos为c，sin为s，将2维向量拓展到多维向量，我们就可以得到将某个元素变为0的T：（T为分块矩阵形式）

$T=\begin{pmatrix} I&&\\ &c&&s\\ &&I&\\ &-s&&c\\ &&&&I \end{pmatrix}$

其中E为单位阵。对于给定n维数向量x，T中两行分别为第j行和第i行，经过Tx得到的n维向量第j个元素为0. （即矩阵的第i行第j列元素化为0）这样的矩阵就叫做Givens矩阵.

QR分解中，R为最后得到的上三角矩阵， $Q=(T_{12}T_{13}T{14}...T_{1n}T_{21}...T_{n-2,n}T_{n-1,n})^T$

二、Householder变换

设 $v$ 为单位向量，矩阵

$H=I-vv^T$

则称 $H$ 为household矩阵.

这个矩阵依旧来源于几何意义：镜面映射矩阵。根据这个几何意义，我们就知道，对于任意的两个向量 $x, y$ ，若x的范数等于y的范数，则可以构造household矩阵 $H$ 实现 $y=Hx$ ，即 $x$ 经过镜面映射后变为 $y$ 。又 $v$ 为单位向量，我们只需要取

$v=\frac{x-y}{||x -y||}$

那household变换又与上三角矩阵有什么关系呢？记上三角矩阵为 $Y=(y_1,y_2,...,y_n)$ ，每一次变换就是把原矩阵的 $x_i$ 变换为 $y_i$ .

$\begin{pmatrix} \times&\times&\times&\times\\ \times&\times&\times&\times\\ \times&\times&\times&\times\\ \times&\times&\times&\times\\ \end{pmatrix} \rightarrow \begin{pmatrix} \times&\times&\times&\times\\ 0&\times&\times&\times\\ 0&\times&\times&\times\\ 0&\times&\times&\times\\ \end{pmatrix} \rightarrow$

$\begin{pmatrix} \times&\times&\times&\times\\ 0&\times&\times&\times\\ 0&\times&\times&\times\\0&\times&\times&\times\\ \end{pmatrix} \rightarrow \begin{pmatrix} \times&\times&\times&\times\\ 0&\times&\times&\times\\ 0&0&\times&\times\\ 0&0&\times&\times\\ \end{pmatrix} \rightarrow$ $\begin{pmatrix} \times&\times&\times&\times\\ 0&\times&\times&\times\\ 0&0&\times&\times\\ 0&0&0&\times\\ \end{pmatrix}$

QR分解中，R为最后得到的上三角矩阵， $Q=H_1H_2H_3...H_n$

三、注意

1. 在得到上三角矩阵的过程中，Givens变换是一个元素一个元素的变化，householder矩阵是一列一列的进行变化。

2. Givens在变化中，一直都是使用原矩阵；而household变换在计算过程中需要不断定义新的矩阵。

四、代码实现

1. 基于Givens变换的QR分解

（1）代码

%基于 Givens 变换的QR分解
%以3阶方阵为例
A = [2 4 1;
    1 5 3;
    3 2 6]
%选取a11与a21（第1行与第2行，变换矩阵记为T12）
x11 = A(1,1);
x21 = A(2,1);
c1 = x11/sqrt(x11^2+x21^2);
s1 = x21/sqrt(x11^2+x21^2);
T12 = [c1 s1 0;
    -s1 c1 0;
    0 0 1];
A1 = T12*A; %将A的a21化为了0

%选取新的a11与a31（第1行与第3行，变换矩阵记为T13）
x12 = A1(1,1);
x22 = A1(3,1);
c2 = x12/sqrt(x12^2+x22^2);
s2 = x22/sqrt(x12^2+x22^2);
T13 = [c2 0 s2;
    0 1 0;
    -s2 0 c2];
A2 = T13*A1; %将A的a31化为了0

%选取新的a22与a32（第2行与第3行，变换矩阵记为T23）
x13 = A2(2,2);
x23 = A2(3,2);
c3 = x13/sqrt(x13^2+x23^2);
s3 = x23/sqrt(x13^2+x23^2);
T23 = [1 0 0;
    0 c3 s3;
    0 -s3 c3];
A3 = T23*A2;%将A的a32化为了0

R = A3
Q = T12*T13*T23

Q*R

（2）运行结果

2. 基于Householder变换的QR分解

（1）代码

%基于Householder变换的QR分解
%以3阶方阵为例
A = [9 3 1;
    8 4 -2;
    2 1 1]
e1 = [1 0 0]';
e2 = [1 0]';

%取第1列，对第1列进行变换
%计算v向量
a1 = A(:,1)
v1 = (a1-norm(a1,2)*e1)/norm(a1-norm(a1,2)*e1,2)
%计算H矩阵
H1 = eye(3) - 2*v1*v1';;%其中eye为单位矩阵
A1 = H1*A

%取新的第2列，对第2列进行变换
%计算v向量
A10 = A1(2:3,2:3);%重新定义新的矩阵,记为A10
a2 = A10(:,1)
v2 = (a2-norm(a2,2)*e2)/norm(a2-norm(a2,2)*e2,2)
%计算H矩阵
H2 = eye(2) - 2*v2*v2'
H2 = vertcat([0,0],H2)';
H2 = vertcat([1,0,0],H2)

A2 = H2*H1*A;

R = A2
Q = H1*H2

Q*R%验证

（2）运行结果对照解释

魔乐社区

魔乐社区（Modelers.cn) 是一个中立、公益的人工智能社区，提供人工智能工具、模型、数据的托管、展示与应用协同服务，为人工智能开发及爱好者搭建开放的学习交流平台。社区通过理事会方式运作，由全产业链共同建设、共同运营、共同享有，推动国产AI生态繁荣发展。

更多推荐

全家桶集齐！Qwen3.5四款小模型上线魔乐社区，附昇腾全套实践教程

魔乐社区

Pont - 搭建前后端之桥：高效、灵活的接口管理工具

Pont 是一款强大的数据服务层解决方案，它能够帮助开发者快速搭建前后端之间的桥梁，实现接口的高效管理和代码自动生成。无论是新手还是有经验的开发者，都能通过 Pont 轻松处理接口文档、生成类型安全的 API 代码，从而显著提升开发效率。[![Pont 工具标志](https://raw.gitcode.com/gh_mirrors/po/pont/raw/3f1b7d4bbba3fd2dda