6.10. Orthogonal Matrix Decompositions

As we learned in the previous sections, Gaussian elimination-based algorithms are the foundation for solving square systems of linear equations. We now present an alternative strategy that provides numerically stable methods to convert rectangular matrices into triangular form. We will multiply the matrix by unitary (orthogonal) transformation matrices to move the matrix elements to zero. The sub-matrices of the factoring are the product of the unitary transformation matrices and the reshaped original matrix. Recall that unitary matrices (Special Matrices) have a condition number of 1 ( Matrix Condition Number), so the multiplications do not degrade the condition of the matrix factors. Orthogonal matrix factoring provides accurate procedures for solving rectangular systems of equations and is also useful in other linear algebra applications beyond solving systems of equations. In this section, we cover the application of orthogonally factored matrices. The factoring algorithms are covered in Orthogonal Matrix Factoring.

This section presents two orthogonal matrix factorings—QR decomposition and the singular value decomposition (SVD). Both decompositions are used to solve rectangular systems of equations. We need to know about vector projections (Projections Onto a Line) before learning how the QR decomposition is found. Similarly, understanding the SVD requires knowledge of eigenvalues and eigenvectors. So this section covers what we need to know about QR and SVD to use them. The QR decomposition algorithm is presented in QR Decomposition. More complete coverage of the SVD is in Singular Value Decomposition (SVD), and the algorithm for finding the SVD matrix factors is in Calculating the SVD with Unitary Transformations.

6.10.1. QR Decomposition

The QR decomposition finds sub-matrices \(\bf{Q}\) and \(\bf{R}\), where \(\bf{Q}\) contains orthogonal column vectors, and \(\bf{R}\) is an upper triangular matrix, such that \(\mathbf{A} = \mathbf{Q}\,\mathbf{R}\). QR decomposition works for both square and rectangular matrices. If the size of \(\bf{A}\) is \(m{\times}n\), the \(\bf{Q}\) factor will be \(m{\times}m\), while \(\bf{R}\) will be \(m{\times}n\). The MATLAB command, [Q,R] = qr(A) finds the factors \(\bf{Q}\) and \(\bf{R}\) of \(\bf{A}\). Function qrFactor is a simple implementation of the algorithm (QR Decomposition).

The QR algorithm finds orthogonal basis vectors from the columns of \(\bf{A}\). The QR algorithm is faster than Gram–Schmidt ( Gram–Schmidt), especially for larger matrices, and has better numerical stability [YETURU20, STRANG16]_.

\[\begin{split}\m{A}_{n{\times}n} = \m{Q\, R} = \begin{bmatrix} \vertbar{} & \vertbar{} & {} & \vertbar{} \\ \bm{q}_1 & \bm{q}_2 & \cdots{} & \bm{q}_n \\ \vertbar{} & \vertbar{} & {} & \vertbar{} \end{bmatrix} \begin{bmatrix} r_{11} & r_{12} & \cdots{} & r_{1n} \\ 0 & r_{22} & \cdots{} & r_{2n} \\ \vdots{} & \vdots{} & \ddots{} & \vdots{} \\ 0 & 0 & \cdots{} & r_{nn} \end{bmatrix}\end{split}\]

One application of QR is to solve a system of square or over-determined equations. [1]

\[\begin{split}\begin{aligned} &\mathbf{A}\,\bm{x} = \mathbf{Q\,R}\,\bm{x} = \bm{b} \\ & \bm{x} = \mathbf{R}^{-1}\,\mathbf{Q}^T\,\bm{b} \end{aligned}\end{split}\]

In MATLAB, the solution x = R \ (Q' * b) is an efficient and numerically stable procedure since \(\bf{Q}\) is orthogonal and \(\bf{R}\) is upper triangular.

6.10.2. Singular Value Decomposition

The singular value decomposition (SVD) has many data analysis applications [MARTIN12] and provides solutions to many linear algebra problems. The linear algebra applications of the SVD include low-rank matrix approximations, finding the rank of a matrix, finding the inverse of a square matrix, finding the pseudo-inverse of rectangular matrices, finding orthogonal basis vectors for the columns of a matrix, finding the condition number of a matrix, and finding the null space (the solution to \(\mathbf{C}\,\bm{y} = \bm{0}\), where \(\bf{C}\) is a singular matrix), and the left null space of a matrix.

The SVD factoring can be found for any matrix. Proof of the existence of the SVD for all matrices is given in The Classic SVD of The Classic SVD. Square, rectangular, complex, and even singular matrices have SVD factors. The factoring is \(\mathbf{A} = \mathbf{U\, \Sigma \,V}^T\). The \(\bf{\Sigma}\) matrix is a diagonal matrix of positive, real, scalar values that we call the singular values. If \(\bf{A}\) is real, both \(\bf{U}\) and \(\bf{V}\) will be real and orthogonal. If \(\bf{A}\) is complex, both \(\bf{U}\) and \(\bf{V}\) will be unitary (the complex version of orthogonal). We will regard \(\bf{A}\) as real because applications of the SVD almost always use real matrices.

The factors look as follows for a square \(n{\times}n\), full-rank matrix.

\[\begin{split}\begin{aligned} \mathbf{A} &= \begin{bmatrix} \vertbar{} & \vertbar{} & {} & \vertbar{} \\ \bm{u}_1 & \bm{u}_2 & \cdots{} & \bm{u}_n \\ \vertbar{} & \vertbar{} & {} & \vertbar{} \end{bmatrix} \begin{bmatrix} \sigma_1 & 0 & \cdots{} & 0 \\ 0 & \sigma_2 & \cdots{} & 0 \\ \vdots{} & \vdots{} & \ddots{} & \vdots{} \\ 0 & 0 & \cdots{} & \sigma_n \end{bmatrix} \begin{bmatrix} \horzbar{} & \bm{v}_1^T & \horzbar{} \\ \horzbar{} & \bm{v}_2^T & \horzbar{} \\ {} & \vdots{} & {} \\ \horzbar{} & \bm{v}_n^T & \horzbar{} \end{bmatrix} \\ &= \mathbf{U\,\Sigma\,V}^T \end{aligned}\end{split}\]

If \(\bf{A}\) is \(m {\times} n\) the factor sizes are \(\bf{U}\): \(m {\times} m\), \(\bf{\Sigma}\): \(m {\times} n\), and \(\bf{V}\): \(n {\times} n\).

Another interesting point about the singular values is that they are ordered (\(\sigma_1 > \sigma_2 > \cdots > \sigma_p\)), where \(p = \min\{m, n\}\). Ordered singular values are useful in data analysis because the larger singular values and corresponding \(\bm{u}\) and \(\bm{v}\) vectors contribute the most to the values in \(\bf{A}\). We will explore several applications of the SVD to linear algebra and data analysis in Singular Value Decomposition (SVD).

Matrix Inverse

An application of the SVD related to solving systems of equations is to compute the inverse of a matrix and the pseudo-inverse of rectangular matrices. Recall that the inverse of a product of matrices is the product of each matrix inverse in reverse order. Since \(\bf{U}\) and \(\bf{V}\) are orthogonal, their inverse is their transpose (Spectral Theorem in Matrix Transpose Properties). Since \(\bf{\Sigma}\) is diagonal and real, the inverse of \(\bf{\Sigma}\) is a diagonal matrix holding the reciprocals of each singular value. Refer to equation (8.4) in Properties of Eigenvalues and Eigenvectors. Any singular values near zero are kept at zero in the inverse. MATLAB’s pinv function can handle the simple task of inverting \(\bf{\Sigma}\) for us.

\[\mathbf{A}^{-1} = \mathbf{V\,\Sigma}^{-1}\,\mathbf{U}^H\]
>> A = randi(20, 4) - 8; %  A random 4-by-4 matrix
>> [U,S,V] = svd(A);
>> invA = V*pinv(S)*U';
>> norm(eye(4) - A*invA, 'fro')
ans =  %  very small difference
   1.7215e-15

Footnote: