zl程序教程

您现在的位置是:首页 >  其它

当前栏目

4.8 QR分解四:向量方程法

向量 分解 方程 4.8 QR
2023-09-14 09:06:54 时间

算法

  除了三大主流算法外,还有向量方程法,向量方程法,顾名思义,就是通过解方程的方式求矩阵的QR分解。其原理是把矩阵 A A A按列向量分块,得到一系列方程,首先有:
A = ( a 1 , a 2 , ⋯   , a n ) Q = ( q 1 , q 2 , ⋯   , q n ) R = ( r 11 r 12 ⋯ r 1 n 0 r 22 ⋯ r 2 n ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ r n n ) ( a 1 , a 2 , ⋯   , a n ) = ( q 1 , q 2 , ⋯   , q n ) ( r 11 r 12 ⋯ r 1 n 0 r 22 ⋯ r 2 n ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ r n n ) A=(\bold{a}_1,\bold{a}_2,\cdots,\bold{a}_n)\\ Q=(\bold{q}_1,\bold{q}_2,\cdots,\bold{q}_n)\\ R=\begin{pmatrix} r_{11} & r_{12} & \cdots & r_{1n}\\ 0 & r_{22} & \cdots & r_{2n}\\ \vdots & \vdots &\ddots &\vdots &\\ 0 & 0 & \cdots & r_{nn}\\ \end{pmatrix}\\ (\bold{a}_1,\bold{a}_2,\cdots,\bold{a}_n)=(\bold{q}_1,\bold{q}_2,\cdots,\bold{q}_n)\begin{pmatrix} r_{11} & r_{12} & \cdots & r_{1n}\\ 0 & r_{22} & \cdots & r_{2n}\\ \vdots & \vdots &\ddots &\vdots &\\ 0 & 0 & \cdots & r_{nn}\\ \end{pmatrix} A=(a1,a2,,an)Q=(q1,q2,,qn)R= r1100r12r220r1nr2nrnn (a1,a2,,an)=(q1,q2,,qn) r1100r12r220r1nr2nrnn
  所以展开后就有以下向量方程:
r 11 q 1 = a 1 r 12 q 1 + r 22 q 2 = a 2 ⋮ r 1 n q 1 + r 2 n q 2 + ⋯ + r n n q n = a n r_{11}\bold{q}_1=\bold{a}_1\\ r_{12}\bold{q}_1+r_{22}\bold{q}_2=\bold{a}_2\\ \vdots\\ r_{1n}\bold{q}_1+r_{2n}\bold{q}_2+\cdots+r_{nn}\bold{q}_n=\bold{a}_n\\ r11q1=a1r12q1+r22q2=a2r1nq1+r2nq2++rnnqn=an
  这个向量方程要从上到下解。首先第一个方程貌似有很多未知数,但是要知道正交矩阵的列向量,长度必须为1的,也就是要单位化,那么 q 1 \bold q_1 q1就是单位化后的 a 1 \bold a_1 a1了。对于第二个方程,需要左乘 q 1 T \bold q_1^T q1T,方程就变成了:
r 12 q 1 T q 1 + r 22 q 1 T q 2 = q 1 T a 2 r_{12}\bold q_1^T\bold{q}_1+r_{22}\bold q_1^T\bold{q}_2=\bold q_1^T\bold{a}_2\\ r12q1Tq1+r22q1Tq2=q1Ta2
  因为正交矩阵的列之间互相正交,所以 q 1 T q 2 \bold q_1^T\bold{q}_2 q1Tq2就是0了,于是有:
r 12 q 1 T q 1 = q 1 T a 2 r_{12}\bold q_1^T\bold{q}_1=\bold q_1^T\bold{a}_2\\ r12q1Tq1=q1Ta2
  这样就求出了 r 12 r_{12} r12,再代入,继续求出 r 22 r_{22} r22。对于最后一个方程:
r 1 n q 1 + r 2 n q 2 + ⋯ + r n n q n = a n r_{1n}\bold{q}_1+r_{2n}\bold{q}_2+\cdots+r_{nn}\bold{q}_n=\bold{a}_n\\ r1nq1+r2nq2++rnnqn=an
  左乘 q i T \bold q_i^T qiT,得到:
r 1 n q i T q 1 + r 2 n q i T q 2 + ⋯ + r i n q i T q i + ⋯ + r n n q i T q n = q i T a n ⇒ r i n q i T q i = q i T a n r_{1n}\bold q_i^T\bold{q}_1+r_{2n}\bold q_i^T\bold{q}_2+\cdots+r_{in}\bold q_i^T\bold{q}_i+\cdots+r_{nn}\bold q_i^T\bold{q}_n=\bold q_i^T\bold{a}_n\\ \Rightarrow r_{in}\bold q_i^T\bold{q}_i=\bold q_i^T\bold{a}_n r1nqiTq1+r2nqiTq2++rinqiTqi++rnnqiTqn=qiTanrinqiTqi=qiTan
  但是得不到 r n n r_{nn} rnn,这个值需要前面的求出来后再代入,才能得到。代入公式为:
r n n q n = a n − r 1 n q 1 − + r 2 n q 2 + ⋯ + r n − 1 , n q n − 1 r_{nn}\bold{q}_n=\bold{a}_n-r_{1n}\bold{q}_1-+r_{2n}\bold{q}_2+\cdots+r_{n-1,n}\bold{q}_{n-1} rnnqn=anr1nq1+r2nq2++rn1,nqn1
  这个算法比较机械,也比较容易实现。但是仔细研究,就会发现,具体算法和施密特正交化差不多。事实上,QR分解是唯一的,四种算法分解的结果都是一样的。

Python代码

  知道了算法之后,直接写代码迭代就完事了。代码如下:

    # 向量方程法QR分解
    def qr_equations(self):
        n = len(self.__vectors)
        q = Matrix([None] * n)
        r = Matrix([[0] * n for _ in range(n)])
        # 按列循环
        for column in range(n):
            # 先将之前的列代入
            r_n = r.__vectors[column]
            for i in range(0, column):
                q_i = q.__vectors[i]
                a_n = self.__vectors[column]
                r_n[i] = inner_product(q_i, a_n) / vector_len(q_i)
            # 最后求得自己
            a_i = self.__vectors[column]
            for i in range(0, column):
                a_i = sub(a_i, mul_num(q.__vectors[i], r_n[i]))
            # 求rnn 与qn
            r_n[column] = vector_len(a_i)
            q_i = mul_num(a_i, 1/r_n[column])
            q.__vectors[column] = q_i
        return q, r

测试数据

  向量方程法,可以分解任意阶矩阵,不限于方阵,如以下矩阵的分解:
( 0 3 1 1 4 − 2 2 1 − 2 1 1 1 ) = ( 0 0.691 0.425 0.408 0.653 − 0.473 0.816 − 0.307 − 0.143 0.408 − 0.038 0.759 ) ( 2.449 2.858 − 2.041 0 4.34 − 0.038 0 0 2.415 ) \begin{pmatrix}0 & 3 & 1\\ 1 & 4 & -2\\ 2 & 1 & -2\\ 1 & 1 & 1\\ \end{pmatrix}=\begin{pmatrix}0 & 0.691 & 0.425\\ 0.408 & 0.653 & -0.473\\ 0.816 & -0.307 & -0.143\\ 0.408 & -0.038 & 0.759\\ \end{pmatrix} \begin{pmatrix}2.449 & 2.858 & -2.041\\ 0 & 4.34 & -0.038\\ 0 & 0 & 2.415\\ \end{pmatrix} 012134111221 = 00.4080.8160.4080.6910.6530.3070.0380.4250.4730.1430.759 2.449002.8584.3402.0410.0382.415