QR_Decomposition_2

QR分解（2）

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

豪斯霍尔的变换法

1. 豪斯霍尔德变换定义

   $$设u\in R^{n}且||u||_{2}=1，矩阵H=I-2uu^{T}称为豪斯霍尔德矩阵$$

2. 豪斯霍尔德矩阵的性质

   • H矩阵是一个正交矩阵

     $$证：\\ HH^{T}=(I-2uu^{T})(I-2uu^{T})\\ =I-4uu^{T}+4uu^{T}uu^{T}\\ =I-4uu^{T}+4u(u^{T}u)u^{T}\\ =I-4uu^{T}+4uu^{T} =I$$

   • 直观的来看，豪斯霍尔德变换是一个镜像变换，如图

     

     豪斯霍尔德矩阵的作用是把向量x映射到以π为镜面的像，其中平面π的法向量为u

基于豪斯霍尔德变换的QR分解

引理：

$$设\ x\in R^{n}是任意的给定非零向量，v\in R^{n}是任意给定的单位向\\量，则存在初等反射矩阵H=I-2uu^{T},使得Hx=\sigma v，其中\sigma为常数$$

证明：

$$\because v为单位向量\\ \therefore x^{T}H^{T}Hx=\sigma^{2}v^{T}v=\sigma^{2}=||x||_{2}^{2}\Rightarrow\\ \sigma=\pm\sqrt{||x||_{2}}\\ 又\because Hx=(I-2uu^{T})x=x-2uu^{T}x=x-2(u^{T}x)u=\sigma v\\ \therefore u的方向与x-\sigma v方向相同\\ \because u是单位向量\\ \therefore u=\frac{x-\sigma x}{||x-\sigma x||_{2}}\\ \therefore u存在且唯一，从而H存在$$

注：其中计算σ时，为避免舍入误差，应取$\sigma=-sgn(x_{1})||x||_{2}$

实际上根据豪斯霍尔德变换的几何意义，我们可以把x看做原向量，v 看成变换后镜像向量的一个单位向量，则u向量方向一定为原向量（x）减去原向量（想）模长倍的镜像向量（v）。单位化后即为所要求的u。

QR分解的过程：

$$设要分解的矩阵为：\\ \quad\\ A= \begin{pmatrix} a_{11}&a_{12}&\dots&a_{1n}\\ a_{21}&a_{22}&\dots&a_{2n}\\ \vdots&\vdots&&\vdots\\ a_{m1}&a_{m2}&\dots&a_{mn} \end{pmatrix}= (a_{1},a_{2},\dots,a_{n})\\ \quad\\ 记为\\ A^{(0)}=(a_{1}^{(0)},a_{2}^{(0)},\dots,a_{n}^{(0)})=A\\ \quad\\ 第一步：\\ 为把矩阵A的第一列化为(\sigma,0,\dots,0)^{T},取x=a_{1}^{(0)},v=e_{1}=(1,0,\dots,0)^{T}\\ 根据式（4），取\\ u_{1}=\frac{a_{1}^{(0)}-\sigma_{1} e_{1}}{||a_{1}^{(0)}-\sigma _{1}e_{1}||_{2}}=\frac{w_{1}}{||w_{1}||_{2}}\\ 其中\sigma_{1}=-sgn(a_{1}^{(0)})||a_{1}^{(0)}||_{2},\ w_{1}=a_{1}^{(0)}-\sigma_{1}e_{1}\\ 而\\ ||w_{1}||_{2}=||a_{1}^{(0)}-\sigma_{1}e_{1}||_{2}=\sqrt{2(\sigma_{1}^{2}-\sigma_{1}a_{11}^{(0)}})\\ 2u_{1}u_{1}^{T}=\frac{2w_{1}w_{1}^{T}}{||w_{1}||_{2}^{2}}=\frac{w_{1}w_{1}^{T}}{\sigma_{1}(\sigma_{1}-a_{11}^{(0)})}=\frac{w_{1}w_{1}^{T}}{\alpha_{1}},\alpha_{1}=\sigma_{1}(\sigma_{1}-a_{11}^{(0)})\\ 令H_{1}=I-2u_{1}u_{1}^{T}=I-\alpha^{-1}w_{1}w_{1}^{T},用H_{1}左乘A^{(0)}得\\ A^{(1)}=H_{1}A^{(0)}=(H_{1}a_{1}^{(0)},H_{1}a_{2}^{(0)},\dots,H_{1}a_{n}^{(0)})\\ =(\sigma_{1} e_{1},a_{2}^{(1)},a_{3}^{(1)},\dots,a_{n}^{(1)})\\ = \begin{pmatrix} \sigma_{1}&a_{12}^{(1)}&a_{13}^{(1)}&\dots&a_{1n}^{(1)}\\ 0&a_{22}^{(1)}&a_{23}^{(1)}&\dots&a_{2n}^{(1)}\\ \vdots&\vdots&\vdots&&\vdots\\ 0&a_{m2}^{(1)}&a_{m3}^{(1)}&\dots&a_{mm}^{(1)} \end{pmatrix}\\ \quad\\ 其中\\ a_{j}^{(1)}=H_{1}a_{j}^{(0)}=a_{j}^{(0)}-\beta_{1j}(a_{1}^{(0)}-\sigma_{1}e_{1})\\ = \begin{pmatrix} a_{1j}^{(0)}-\beta_{1j}(a_{11}^{(0)}-\sigma_{1})\\ a_{2j}^{(0)}-\beta_{1j}a_{21}^{(0)}\\ \vdots\\ a_{mj}^{(0)}-\beta_{1j}a_{mj}^{(0)} \end{pmatrix}\\ 由此可得\\ \left\{ \begin{array}{} a_{1j}^{(1)}=a_{1j}^{(0)}-\beta_{1j}(a_{11}^{(0)}-\sigma_{1}),\quad j=2,3,\dots,n\\ a_{ij}^{(1)}=a_{ij}^{(0)}-\beta_{1j}a_{i1}^{(0)},\quad i=2,3,\dots,m;j=2,3,\dots,n \end{array} \right.\\ 其中\\ \beta_{ij}=\alpha_{1}^{-1}((a_{11}^{(0)}-\sigma_{1})a_{1j}^{(0)}+\sum_{i=2}^{m}a_{i1}^{(0)}a_{ij}^{(0)}),\quad j=2,3,\dots,n\\ \quad\\ 第二步:\\ 记a_{2}^{(1)}=(a_{12}^{(1)},a_{22}^{(1)},\dots,a_{m2}^{(1)})^{T}=(a_{12}^{(1)},\widetilde a_{2}^{(1)T})^{T}\\ 重复第一步取\\ u_{2}=\frac{\widetilde a_{2}^{(1)}-\sigma_{2}e_{1}^{(1)}}{||\widetilde a_{2}^{(1)}-\sigma_{2}e_{1}^{(1)}||_{2}}=\frac{w_{2}}{||w_{2}||_{2}}\\ \sigma_{2}=-sgn(a_{22}^{(1)})||\widetilde a_{2}^{(1)}||_{2};w_{2}=\widetilde a_{2}^{(1)}-\sigma_{2}e_{1}^{(1)}\\ 构造\\ \widetilde {H_{2}}=I-2u_{2}u_{2}^{T}=I-\alpha^{-1}w_{2}w_{2}^{T} 令\\ H_{2}= \begin{pmatrix} 1&0\\ 0&\widetilde H_{2} \end{pmatrix}\\ 用H_{2}左乘A_{(2)}得\\ A^{(2)}=H_{2}A^{(1)}=(\sigma_{1}e_{1},a_{12}^{(1)}e_{1}+\sigma_{2}e_{2},a_{3}^{(2)},\dots,a_{n}^{(2)})\\ =\begin{pmatrix} \sigma_{1}&a_{12}^{(1)}&a_{13}^{(1)}&a_{14}^{(1)}&\dots&a_{1n}^{(1)}\\ 0&\sigma_{2}&a_{23}^{(2)}&a_{24}^{(2)}&\dots&a_{2n}^{(2)}\\ 0&0&a_{33}^{(2)}&a_{34}^{(2)}&\dots&a_{3n}^{(2)}\\ \vdots&\vdots&\vdots&\vdots&&\vdots\\ 0&0&a_{m3}^{(2)}&a_{m4}^{(2)}&\dots&a_{mn}^{(2)} \end{pmatrix}\\ 继续重复这个过程直至s=min{m-1,n}步，最后得到\\ A^{(s)}=H_{s}H_{s-1}\dots H_{1}A=R\\ 记\\ H=H_{s}H_{s-1}\dots H_{1}A\\ 则\\ HA=R \Rightarrow A=QR$$

• 当m=n时,R是n阶可逆上三角矩阵

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

$\bigstar$ 矩阵的QR分解具有唯一性

$$证明\\ 设存在两种分解 \quad A=QR=\widetilde Q \widetilde R\\ \because A^{T}A=R^{T}QQR=R^{T}R\\ A^{T}A=\widetilde R^{T} \widetilde Q \widetilde Q \widetilde R=\widetilde R^{T} \widetilde R\\ 又\because A^{T}A是对称正定的，且矩阵的楚列斯基分解是唯一的\\ \therefore R=\widetilde R,\quad Q=\widetilde Q$$

$\bigstar$ QR分解常用于求解病态线性方程组Ax=b，即化为

$$y=Q^{T}b,\quad Rx=y$$

通过这种方法求解线性方程是非常稳定的，得到的结果往往比三角分解好很多。

豪斯霍尔德法QR分解的程序

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();
        }
		//单位阶跃函数
        public static int Sgn(float num)
        {
            if (num < 0)
            {
                return -1;
            }
            else if (num == 0)
            {
                return 0;
            }
            else
            {
                return 1;
            }
        }
    \\类的封装
        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)
                {
                    for (int i = 0; i < row; i++)
                    {
                        for (int j = 0; j < column; 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 static Matrix operator -(Matrix p1, Matrix p2)
            {
                if (p1.row == p2.row && p1.column == p2.column)
                {
                    Matrix result = new Matrix(p1.row, p1.column, 0);
                    for (int i = 0; i < p1.row; i++)
                    {
                        for (int j = 0; j < p1.column; j++)
                        {
                            result.matrix[i, j] = p1.matrix[i, j] - p2.matrix[i, j];
                        }
                    }
                    return result;
                }
                else
                {
                    throw new Exception();

                }
            }
            public static Matrix operator *(float multip, Matrix p1)
            {

                Matrix result = new Matrix(p1.row, p1.column, 0);
                for (int i = 0; i < p1.row; i++)
                {

                    for (int j = 0; j < p1.column; j++)
                    {

                        result.matrix[i, j] = p1.matrix[i, j] * multip;

                    }

                }
                return result;

            }
            public static Matrix operator *(Matrix p1, Matrix p2)
            {
                if (p1.column != p2.row)
                {
                    //Matrix a = new Matrix(1, 2) { matrix = p1.matrix };
                    throw new Exception();
                }
                else
                {
                    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;
                }
            }
			//矩阵转置
            public void Transpose()
            {
                this.row = column + row;
                this.column = row - column;
                this.row = row - column;
                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 Matrix ReflectMatrix(int order, Matrix column_vec, int suborder)
            {
                Matrix Reflect_matrix = new Matrix(order, order, 1);
                \\生成单位向量
                Matrix unit_vec = new Matrix(column_vec.row, column_vec.column, 1);
                \\生成单位矩阵
                Matrix unit_matrix = new Matrix(suborder, suborder, 1);
                
                \\需要使用临时变量存储转置矩阵（临时变量初始化不能直接赋值column_vec）,原因是：matrix是引用类型，赋值的方式只是使变量指向同一个对象，进行转置操作时会丢失原来的矩阵，下同
                    
                Matrix temp_1 = new Matrix(column_vec.row, column_vec.column, 0);
                temp_1.matrix = column_vec.matrix;
                temp_1.Transpose();
                \\根据算法，计算相应参数
                float sigma = -Sgn(column_vec.matrix[0, 0]) * (float)Math.Sqrt((temp_1 * column_vec).matrix[0, 0]);
                float alpha = sigma * (sigma - column_vec.matrix[0, 0]);
                Matrix omega = column_vec - sigma * unit_vec;
                Matrix temp_2 = new Matrix(omega.row, omega.column, 0);
                temp_2.matrix = omega.matrix;
                temp_2.Transpose();
                if (order == suborder)
                {
                    Reflect_matrix = unit_matrix - (1 / alpha) * omega * temp_2;
                }
                else
                {
                    Matrix temp_matrix = unit_matrix - (1 / alpha) * omega * temp_2;
                    for (int i = order - suborder; i < order; i++)
                    {
                        int start_value = order - suborder;
                        for (int j = start_value; j < order; j++)
                        {
                            Reflect_matrix.matrix[i, j] = temp_matrix.matrix[i - start_value, j - start_value];
                        }
                    }
                }
                return Reflect_matrix;
            }
			//QR分解主程序
             public Matrix[] QR_Decomposition()
            {
				//准备储存结果的容器
                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++)
                {
                    Matrix column_vec = new Matrix(row - i, 1, 0);
                    for (int j = 0; j < row - i; j++)
                    {
                        column_vec.matrix[j, 0] = result_R.matrix[j + i, i];
                    }
                    Matrix reflect_matrix = ReflectMatrix(row, column_vec, row - i);
                    result_Q = reflect_matrix * result_Q;
                    result_R = reflect_matrix * result_R;
                }
                result_Q.Transpose();
                Matrix[] result = { result_Q, result_R };
                return result;
            }
			//仅用于矩阵的打印
            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");
                }
            }
        }
    }

同样使用两个例子验证算法的正确性：

1. 输入矩阵

   $$A= \begin{pmatrix} 3&5&5\\ 0&3&4\\ 4&0&5 \end{pmatrix}$$

   得到输出如下图

   

2. 输入矩阵

   $$B= \begin{pmatrix} 3 &14 &9\\ 6 &43 &3\\ 6 &22 &15\\ \end{pmatrix}$$

   得到结果如下图：

   

   经过检验，基本可以确定程序是正确的。