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 where forward and back substitution algorithms are able to quickly find a solution. We will multiply the matrix with unitary (also called orthogonal) transformation matrices to move matrix elements to zeros. The sub-matrices of the factoring are the product of the unitary transformation matrices and what is left of the original matrix. Recall that unitary matrices (Special Matrices) have a condition number of 1 (Matrix Condition Number), so the condition of the matrix factors is not degraded by the multiplications. Orthogonal matrix factoring provides accurate procedures for solving rectangular systems of equations and are also useful to other linear algebra applications beyond solving systems of equations.

Two orthogonal matrix factorings are presented in this section—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, a full understanding of 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 algorithm for computing the QR decomposition is presented in The Unitary Transformation Matrices. More complete coverage of the SVD, including algorithms for finding the matrix factors, is in Singular Value Decomposition (SVD).

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).

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

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

\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}

In MATLAB, the solution is x = R \ (Q' * b). This is an efficient and numerically stable procedure since \bf{Q} is orthogonal and \bf{R} is upper triangular, which can be quickly resolved with back substitution.

6.10.2. The Singular Value Decomposition

The singular value decomposition (SVD) has many data analysis applications 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. Square, rectangular, complex, and even singular matrices all have SVD factors. The factoring is \mathbf{A} = \mathbf{U\, \Sigma \,V}^T. Both \bf{U} and \bf{V} are orthogonal matrices. The \bf{\Sigma} matrix is a diagonal matrix of positive, real, scalar values that we call the singular values.

For a square n {\times} n, full-rank, real matrix, the factors look as follows.

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

The size of the factors are as follows.

\bf{A}:
m {\times} n
\bf{U}:
m {\times} m
\bf{\Sigma}:
m {\times} n, The number of nonzero singular values in \bf{\Sigma} is the rank (r) of \bf{A}. If m = n = r then each diagonal element of the square matrix will have a nonzero singular value. If m > r, then m - r values on diagonal of \bf{\Sigma} will be zeros. If n > r, then n - r values on diagonal of \bf{\Sigma} will be zeros.
\bf{V}:
n {\times} n

Another interesting point about the singular values is that they are ordered, (\sigma_1 > \sigma_2 > \cdots > \sigma_n), which is 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}. It is also useful for finding a null solution of a singular matrix because the zero or near zero singular values are last in the diagonal sequence.

Matrix Inverse

A good application of the SVD is to use the SVD to compute the inverse of a matrix. Recall that the inverse of a product of matrices is the product of each matrix inverse in reverse the order. Since \bf{U} and \bf{V} are orthogonal, their inverse is just their transpose. Since \bf{\Sigma} is diagonal and real, the inverse of \bf{\Sigma} is a diagonal matrix holding the reciprocals of each singular value. Any singular values near zero are kept at zero in the inverse. MATLAB’s pinv 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';
>> invA - inv(A)
ans =
    1.0e-14 *  %  very small difference
    0.0194   -0.0069    0.0416   -0.0333
   -0.0097    0.0042   -0.0167    0.0139
    0.0722   -0.0153    0.1443   -0.0888
   -0.0500    0.0167   -0.1110    0.0833
Rank from the SVD

There are multiple ways to compute the rank of a matrix. An example of using Gaussian elimination to find the rank of a matrix is shown in Rank. The MATLAB rank function uses the SVD. The number of singular values that are greater than a small tolerance is the rank of the matrix. According to the MATLAB documentation, the SVD algorithm requires more computation than some alternatives, but it is the most reliable.

In the following example, we start with a full rank matrix and make the third row the difference of the first two rows. The final singular value is close to zero and the rank of the singular matrix is 2.

>> A(3,:) = A(1,:) - A(2,:)
A =
     6     3     6
     4    -3     3
     2     6     3
>> rank(A)
ans =
     2
>> svd(A)
ans =
   11.1167
    6.3576
    0.0000

>> fprintf("%g/n", ans(3))
1.32734e-15
The Null Solution

As mentioned in The Row and Column View, only singular systems have a solution to the equation \mathbf{C}\,\bm{y} = \bm{0}, which is called the null solution. Finding the null space from both elimination and the SVD is covered in detail in Null Space. With our introduction to the SVD here, we will make a quick observation about the columns of the \bf{V} matrix in relation to the null solution.

If we consider the SVD of a singular system as \mathbf{A\,V} = \mathbf{U\,\Sigma}, then we can split the \bf{V}, \bf{\Sigma}, and \bf{U} matrices into blocks related to nonzero and zero valued singular values.

\begin{aligned} &\mathbf{A}\,\mat{\mathbf{V}_r \mathbf{V}_0} = \mat{\mathbf{U}_r \mathbf{U}_0} \mat{\mathbf{\Sigma}; \mathbf{0}} \\ &\mathbf{A}\,\mathbf{V}_r = \mathbf{U}_r\,\mathbf{\Sigma}, \quad \text{and} \quad \mathbf{A}\,\mathbf{V}_0 = \mathbf{U}_0\,\mathbf{0} = \mathbf{0} \end{aligned}

The last n - r columns of \bf{V} is the null solution because \mathbf{A}\,\mathbf{V}_0 = \mathbf{0}. MATLAB’s null function compares the singular values to a threshold to decide if and how many columns of \bf{V} constitutes the null solution. However, determining if a matrix is singular, or just close to singular is sometimes vague. The null function will sometimes not return a solution for nearly singular matrices. The SVD can be used in these cases to find the borderline null solution.

Footnote:

[1]QR can also provide a solution to an under-determined system of equations, but it is usually solved by passing the transpose to QR and then working with transposed matrix factors. See Left-Divide of Under-determined Systems.