QR decomposition

Numerical Analysis

QR分解(1)

QR分解是线性代数中一种常用的矩阵分解，即把矩阵分解为一个正规正交矩阵和一个上梯形矩阵——QR分解是求一般矩阵全部特征值的最有效并广泛应用的方法，同时也是许多迭代算法的基础。

QR分解的导出

QR分解可以通过两种基本的正交变换来导出——吉文斯变换和豪斯霍尔德变换，下面将详细地给出两种正交变换推导QR分解的过程。

吉文斯变换法

  1. 吉文斯变换的定义

  2. 吉文斯变换的特性

    • 通过计算我们可以得到$P(i,j,\theta)\times P(i,j,\theta)^{T}=I$，因此吉文斯变换矩阵是一个正交矩阵。

    • 当吉文斯矩阵左乘一个方阵x时（$P(i,j,\theta)x$），只会改变x的第i个和第j个分量。

    • 直观的看，吉文斯变换是一个平面旋转矩阵。由上一条性质，我们知道吉文斯矩阵只改变x的i,j分量，这种改变的几何意义是：它会把x在i,j平面的投影顺时针旋转$\theta$度，以二维的吉文斯矩阵为例：

    • 吉文斯矩阵可以将给定向量的某个分量化为零，由第二个特性，我们可以得到以下结论：

 也可以直观的理解为将给定向量在i,j平面的分量旋转至与i坐标轴重合，显然这种旋转存在且旋转角度就是上式所给出的。

  1. 基于吉文斯变换的QR分解

    思路：由于QR分解的形式是分解为一个上梯形矩阵和一个正交矩阵，因此可以考虑用若干个吉文斯矩阵的乘积将给定矩阵化为上三角矩阵。

引理：设x是任意给定的非零向量，则存在有限个吉文斯矩阵的乘积P，使得$Px=||x||_{2}\vec e_{1}$，其中$\vec e_{i}=(1,0,\dots,0)$

证明：

根据吉文斯矩阵的最后一条特性，依次构造

根据引理我们可进行基于吉文斯变换的QR分解：

设A是一个$m\times n$的矩阵且它的列向量线性无关，则利用吉文斯变换可以把A逐步化为上梯形矩阵。

继续上述变换过程，直至$s=min\{m-1,n\}$步，最后得到

由于上述的过程中每个$P_{i}$都是正交矩阵，所以$P=P_{s}P_{s-1}\dots P_{1}$也是正交矩阵，所以式（7）可以写为

  • 当m=n时，R是N阶上三角矩阵，这时式（8）可以写为

​ 其中$Q=P^{T}$为n阶正交矩阵，R为可逆的上三角矩阵

  • 当m>n时，R是$m\times n$阶的上梯形矩阵，此时

其中$R_{1}$为可逆的上三角矩阵

至此，我们实现了矩阵A基于吉文斯变换的QR分解。

吉文斯变换法QR分解的程序

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
namespace QR_decomposition
{
    class Program
    {
        static void Main(string[] args)
        {
            //用于验证的demo，自行修改
            Matrix m = new Matrix(3, 3);
            Matrix[] result_QR = m.QR_Decomposition();
            result_QR[0].DisPlay();
            Console.WriteLine();
            result_QR[1].DisPlay();
            Matrix n = result_QR[0] * result_QR[1];
            Console.WriteLine();
            n.DisPlay();

        }
        //对象封装
        class Matrix
        {
            //数据和矩阵形状的信息
            public int row = 0;
            public int column = 0;
            private float[,] matrix;
            //不同的构造方法，输入初始化和自动初始化
            public Matrix(int row, int column, int option)
            {
                this.row = row;
                this.column = column;
                matrix = new float[row, column];
                if (option == 1 && row == column)
                {
                    for (int i = 0; i < row; i++)
                    {
                        for (int j = 0; j < row; j++)
                        {
                            if (i == j)
                            {
                                matrix[i, i] = 1;
                            }
                            else
                            {
                                matrix[i, j] = 0;
                            }
                        }
                    }
                }

            }
            public Matrix(int row, int column)
            {
                this.row = row;
                this.column = column;
                Console.WriteLine($"Please enter a matrix with {this.row}rows,{this.column}colums");
                matrix = new float[row, column];
                try
                {
                    for (int i = 0; i < row; i++)
                    {
                        for (int j = 0; j < column; j++)
                        {
                            string mid = Console.ReadLine();
                            var rows = mid.Split(' ');
                            foreach (var item in rows)
                            {
                                matrix[i, j++] = float.Parse(item);
                            }

                        }

                    }
                    Console.WriteLine("\r\n");
                }
                catch (Exception e)
                {
                    Console.WriteLine($"something wrong!\n{e.Message}");
                }
            }
            //打印矩阵
            public void DisPlay()
            {
                for (int i = 0; i < this.row; i++)
                {
                    for (int j = 0; j < this.column; j++)
                    {
                        Console.Write(string.Format("{0:F3}", matrix[i, j]) + "   ");
                    }
                    Console.WriteLine("\r");
                }
            }
            
            //构造吉文斯矩阵的函数
            public Matrix RotateMatrix(int order, float para1, float para2, int index1, int index2)
            {
                Matrix Rotate_matrix = new Matrix(order, order, 1);
                float cos = (float)(para1 / (Math.Sqrt(Math.Pow(para1, 2) + Math.Pow(para2, 2))));
                float sin = (float)(para2 / (Math.Sqrt(Math.Pow(para1, 2) + Math.Pow(para2, 2))));
                Rotate_matrix.matrix[index1, index1] = cos;
                Rotate_matrix.matrix[index2, index2] = cos;
                Rotate_matrix.matrix[index1, index2] = sin;
                Rotate_matrix.matrix[index2, index1] = -sin;
                return Rotate_matrix;
            }
            //矩阵转置
            public void Transpose()
            {
                float[,] temp=new float[row,column];
                for (int i = 0; i <this. row; i++)
                {
                    for (int j = 0; j <this. column; j++)
                    {
                        temp[i, j] = matrix[j, i];
                    }
                }
                matrix = temp;

            }
            //矩阵乘法的复写
            public static Matrix operator *(Matrix p1, Matrix p2)
            {
                Matrix result = new Matrix(p1.row, p2.column, 0);
                for (int i = 0; i < p1.row; i++)
                {

                    for (int j = 0; j < p2.column; j++)
                    {
                        float temp = 0;
                        for (int k = 0; k < p1.column; k++)
                        {
                            temp += p1.matrix[i, k] * p2.matrix[k, j];

                        }
                        result.matrix[i, j] = temp;
                    }

                }
                return result;
            }
            //QR分解的主要过程
            public Matrix[] QR_Decomposition()
            {
				//申明两个矩阵用于存放结果Q矩阵和R矩阵
                Matrix result_Q = new Matrix(this.column, this.column, 1);
                Matrix result_R = new Matrix(this.row, this.column, 0);
                result_R.matrix = this.matrix;
                //按照算法先列后行进行扫描
                for (int i = 0; i < column; i++)
                {
                    for (int j = i + 1; j < row; j++)
                    {
                        //遇到元素是零的情况，跳过此次循环
                        if (result_R.matrix[j, i] == 0) continue;
                        Matrix rotate_matrix = RotateMatrix(column, result_R.matrix[i, i], result_R.matrix[j, i], i, j);
                        result_Q = rotate_matrix * result_Q;
                        result_R = rotate_matrix * result_R;
                    }
                }
                result_Q.Transpose();
                Matrix[] result = { result_Q, result_R };
                return result;
            }
        }

    }

}

程序的验证

对于上面的程序，我们用两个例子来进行验证：

  • 输入矩阵

    我们将得到输出如下

    mark

    经过验证程序输出的QR分解是正确的，并且通过计算QR矩阵的乘积，我们得到了原来输入的矩阵，这足以证明程序的正确性

  • 再次输入矩阵

    输出结果如下

    mark

    经过检验，结果依然是正确的。

本文标题:QR decomposition

文章作者:ZhuMing

发布时间:2019-08-08, 16:22:38

最后更新:2019-08-11, 11:26:10

原始链接:http://zhumingcj.github.io/2019/08/08/qr-decomposition_1/

许可协议: "署名-非商用-相同方式共享 4.0" 转载请保留原文链接及作者。

文章目录
  1. 1. QR分解(1)
    1. 1.1. QR分解的导出
      1. 1.1.1. 吉文斯变换法
      2. 1.1.2. 吉文斯变换法QR分解的程序
        1. 1.1.2.1. 程序的验证