5.7 QR分解三：Givens旋转

醒过来摸鱼

已于 2023-09-14 14:55:37 修改

阅读量1.7k

点赞数 2

分类专栏：线性代数【更新中】文章标签：线性代数矩阵算法

于 2023-02-03 14:55:20 首次发布

本文链接：https://blog.csdn.net/m0_66201040/article/details/128860777

版权

线性代数【更新中】专栏收录该内容

71 篇文章 42 订阅

订阅专栏

文章目录

Givens矩阵
旋转到标准基
QR分解步骤
Python代码
测试例子

Givens矩阵

Givens矩阵有三个参数组成， $i,j,\theta$ .它长这样：
$G(i,j,\theta)=\begin{pmatrix} I & 0 & 0 & 0 & 0\\ 0 & \cos(\theta) & 0 & \sin(\theta) & 0 \\ 0 & 0 & I & 0 & 0\\ 0 & -\sin(\theta) & 0 & \cos(\theta) & 0\\ 0 & 0 & 0 & 0 & I \end{pmatrix}\\$
如果角度不好计算也可以把 $\cos \theta$ 和 $\sin \theta$ 改成 $c, s$ ，于是Given矩阵就由四个参数定义：
$G(i,j,c,s)=\begin{pmatrix} I & 0 & 0 & 0 & 0\\ 0 & c & 0 & s & 0 \\ 0 & 0 & I & 0 & 0\\ 0 & -s & 0 & c & 0\\ 0 & 0 & 0 & 0 & I \end{pmatrix}\\$
整个矩阵除了对角线和 $a_{ii},a_{ij},a_{ji},a_{jj}$ ，其余位置都是0。对角线上除了第 $i$ 行和第 $j$ 行，其余位置都是1.
Givens矩阵和Householder差不多，都可以把向量变成第一个坐标为它的长度，其余坐标变成0.不同的是这样变，只需要一个Householder矩阵，而Givens变换需要多个Givens矩阵。

旋转到标准基

回忆利用Householder变换进行QR分解就是将向量投影到标准基上。Givens旋转也是一样，不过旋转比较麻烦。把任意向量x旋转到标准基 $e_1$ 上，这个旋转变换是 $n$ 个Givens矩阵的乘积。乘积的每一项定义如下：
$T_{1i}=G(1,i,\frac{\sqrt{\sum_{j=1}^{i-1}x_j^2}}{\sqrt{\sum_{j=1}^{i}x_j^2}},\frac{a_i}{\sqrt{\sum_{j=1}^{i}x_j^2}})$
把所有的这些Givens矩阵倒序乘起来得到一个矩阵：
$T=\prod_{i=n}^2T_{1i}$
这个矩阵乘以向量会把向量旋转到 $e_1$ 方向，也就是第一个坐标为向量的长度，其余坐标为0.说得这么多不如举个例子:
$x=\begin{pmatrix}3\\ 4\\ 3\ 4\\ 5\\ \end{pmatrix}$
分解步骤：
$T_{12}=\begin{pmatrix}0.6 & 0.8 & 0 & 0 & 0\\ -0.8 & 0.6 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 1\\ \end{pmatrix}\\ T_{13}=\begin{pmatrix}0.857 & 0 & 0.514 & 0 & 0\\ 0 & 1 & 0 & 0 & 0\\ -0.514 & 0 & 0.857 & 0 & 0\\ 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 1\\ \end{pmatrix}\\ T_{14}=\begin{pmatrix}0.825 & 0 & 0 & 0.566 & 0\\ 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0\\ -0.566 & 0 & 0 & 0.825 & 0\\ 0 & 0 & 0 & 0 & 1\\ \end{pmatrix}\\ T_{15}= \begin{pmatrix}0.816 & 0 & 0 & 0 & 0.577\\ 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 0\\ -0.577 & 0 & 0 & 0 & 0.816\\ \end{pmatrix}\\ T=T_{15}T_{14}T_{13}T_{12}=\begin{pmatrix}0.346 & 0.462 & 0.346 & 0.462 & 0.577\\ -0.8 & 0.6 & 0 & 0 & 0\\ -0.309 & -0.412 & 0.857 & 0 & 0\\ -0.291 & -0.388 & -0.291 & 0.825 & 0\\ -0.245 & -0.327 & -0.245 & -0.327 & 0.816\\ \end{pmatrix}\\ Tx=\begin{pmatrix}8.66\\ 0\\ 0\\ 0\\ 0\\ \end{pmatrix}$

QR分解步骤

思想和Householder变换一样，按列循环，对于第一列找到一个 $T_1$ 把矩阵的第一列变成上三角矩阵的第一列。对于第二列，就用同样的办法，生成Givens矩阵，连乘起来变成变换阵，把第二行以下的数字变成0,但是这个时候Givens矩阵的 $c, s$ 参数的计算方式要变了，变成这个样子：
$T_{ki}=G(k,i,\frac{\sqrt{\sum_{j=k}^{i-1}x_j^2}}{\sqrt{\sum_{j=k}^{i}x_j^2}},\frac{a_i}{\sqrt{\sum_{j=k}^{i}x_j^2}})$
还是以这个向量为例子，通过一系列Givens矩阵可以将其第二行以下变成0：
$T_{23}=\begin{pmatrix}1 & 0 & 0 & 0 & 0\\ 0 & 0.8 & 0.6 & 0 & 0\\ 0 & -0.6 & 0.8 & 0 & 0\\ 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 1\\ \end{pmatrix}\\ T_{24}=\begin{pmatrix}1 & 0 & 0 & 0 & 0\\ 0 & 0.781 & 0 & 0.625 & 0\\ 0 & 0 & 1 & 0 & 0\\ 0 & -0.625 & 0 & 0.781 & 0\\ 0 & 0 & 0 & 0 & 1\\ \end{pmatrix}\\ T_{25}= \begin{pmatrix}1 & 0 & 0 & 0 & 0\\ 0 & 0.788 & 0 & 0 & 0.615\\ 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 0\\ 0 & -0.615 & 0 & 0 & 0.788\\ \end{pmatrix}\\ T_2=T_{25}T_{24}T_{23}=\begin{pmatrix}1 & 0 & 0 & 0 & 0\\ 0 & 0.492 & 0.369 & 0.492 & 0.615\\ 0 & -0.6 & 0.8 & 0 & 0\\ 0 & -0.5 & -0.375 & 0.781 & 0\\ 0 & -0.384 & -0.288 & -0.384 & 0.788\\ \end{pmatrix}\\ T_2x=\begin{pmatrix}3\\ 8.124\\ 0\\ 0\\ 0\\ \end{pmatrix}$
所以完整的QR分解就是这样按列迭代。但是要注意的是这一系列迭代是这样的过程：
$T_{n-1}\cdots T_1A=R\\ A=T_1^{-1}T_2^{-1}\cdots T_{n-1}^{-1}R=QR$
Givens矩阵是正交阵，他们的乘积也是正交阵，所以有：
$T_i^{-1}=T_i^T$
所以有：
$Q=T_1^TT_2^T\cdots T_{n-1}^T$

Python代码

经过上述讲解，代码写起来就容易了,最终代码如下：

   @staticmethod
    def givens_matrix(n, i, j, c, s):
        array = Matrix.unit_matrix(n)
        array[i][i] = c
        array[j][j] = c

        array[i][j] = -s
        array[j][i] = s

        return Matrix(array)

    @staticmethod
    def get_givens_matrix(vector, line, i):
        sum_vector = 0
        for j in range(line, i):
            sum_vector += vector[j] * vector[j]
        numerator = math.sqrt(sum_vector)
        vector_i_pow = vector[i] * vector[i]
        denominator = math.sqrt(sum_vector + vector_i_pow)
        alpha = numerator / denominator
        beta = math.sqrt(vector_i_pow) / denominator
        return Matrix.givens_matrix(len(vector), line, i, alpha, beta)

    @staticmethod
    def get_givens_rotator(vector, line):
        rotator = Matrix(Matrix.unit_matrix(len(vector)))
        n = len(vector)
        for i in range(n - 1, line, -1):
            givens_matrix = Matrix.get_givens_matrix(vector, line, i)
            rotator = rotator * givens_matrix
        return rotator

    def givens_qr(self):
        n = len(self.__vectors)
        r = self
        q = Matrix(Matrix.unit_matrix(n))
        for i in range(n-1):
            vector = r.__vectors[i]
            givens = Matrix.get_givens_rotator(vector, i)
            # 正交矩阵的逆矩阵就是它的转置矩阵
            q = q * givens.transpose_matrix()
            r = givens * r
        return q, r