20.1. Basic idea: reduction to triangular systems of equations

Consider the problem of solving a system of linear equations [latex]Ax = y[/latex], where [latex]A \in \mathbb{R}^{m \times n}[/latex] and [latex]y \in \mathbb{R}^m[/latex] are given.

The basic idea in the solution algorithm starts with the observation that in the special case when [latex]A[/latex] is upper triangular, that is, [latex]A_{ij} = 0[/latex] if [latex]i j[/latex], then the system can be easily solved by a process known as backward substitution. In backward substitution we simply start solving the system by eliminating the last variable first, then proceed to solve backwards. The process is illustrated in this example, and described in generality here.

20.2. The QR decomposition of a matrix

The QR decomposition allows to express any [latex]m \times n[/latex] matrix [latex]A[/latex] as the product [latex]A = QR[/latex] where [latex]Q[/latex] is [latex]m \times m[/latex] and orthogonal (that is, [latex]Q^TQ = I_m[/latex]) and [latex]R[/latex] is an [latex]m \times m[/latex] upper triangular. For more details on this, see here.

Once the QR factorization of [latex]A[/latex] is obtained, we can solve the system by first pre-multiplying with [latex]Q[/latex] both sides of the equation:

[latex]\begin{align*} QRx = y \iff Rx = Q^Ty. \end{align*}[/latex]

This is due to the fact that [latex]Q^TQ = I_m[/latex]. The new system [latex]Rx = Q^Ty[/latex] is triangular and can be solved by backwards substitution. For example, if [latex]A[/latex] is full column rank, then [latex]R[/latex] is invertible, so that the solution is unique, and given by [latex]x = R^{-1}Q^Ty.[/latex].

Let us detail the process now.

20.3. Using the full QR decomposition

We start with the full QR decomposition of A with column permutations:

[latex]\begin{align*}  AP = QR = \begin{pmatrix} Q_1 & Q_2 \end{pmatrix} \begin{pmatrix} R_1 & R_2 \\ 0 & 0 \end{pmatrix} \end{align*}[/latex]

where

• [latex]Q = [Q_1, Q_2][/latex] is [latex]m \times m[/latex] and orthogonal ([latex]Q^TQ = I_m[/latex]);
• [latex]Q_1[/latex] is [latex]m \times r[/latex], with orthonormal columns ([latex]Q^TQ = I_r[/latex]);
• [latex]Q_2[/latex] is [latex]m \times (m-r)[/latex], with orthonormal columns ([latex]Q^TQ = I_{m - r}[/latex]);
• [latex]r[/latex] is the rank of [latex]A[/latex];
• [latex]R_1[/latex] is [latex]r \times r[/latex] upper triangular, and invertible;
• [latex]R_2[/latex] is a [latex]r \times (n - r)[/latex] matrix;
• [latex]P[/latex] is a [latex]n \times n[/latex] permutation matrix (thus, [latex]P^T = P^{-1}[/latex]).
• The zero submatrices in the bottom (block) row of [latex]R[/latex] have [latex]m - r[/latex] rows.

Using [latex]A = QRP^T[/latex], we can write [latex]Rz = Q^Ty[/latex], where [latex]z := P^Tx[/latex]. Let’s look at the equation in [latex]z[/latex] in expanded form:

[latex]\begin{align*} \begin{pmatrix} R_1 & R_2 \\ 0 & 0 \end{pmatrix} \begin{pmatrix} z_1 \\ z_2 \end{pmatrix} = \begin{pmatrix} Q_1^Ty \\ Q_2^Ty \end{pmatrix}. \end{align*}[/latex]

We see that unless [latex]Q_2^Ty = 0[/latex], there is no solution. Let us assume that  [latex]Q_2^Ty = 0[/latex]. We have then

[latex]\begin{align*} R_1z_1 + R_2z_2 = Q_1^Ty, \end{align*}[/latex]

which is a set of [latex]r[/latex] linear equations in [latex]n[/latex] variables.

A particular solution is obtained upon setting [latex]z_2 = 0[/latex], which leads to a triangular system in [latex]z_1[/latex], with an invertible triangular matrix [latex]R_1[/latex]. So that [latex]z_1 = R_1^{-1}Q_1^Ty[/latex], which corresponds to a particular solution [latex]x_0[/latex] to [latex]Ax = y[/latex] :

[latex]\begin{align*} x_0 := P\begin{pmatrix} R_1^{-1}Q_1^Ty\\ 0. \end{pmatrix} \end{align*}[/latex]

20.4. Set of solutions

We can also generate all the solutions, by noting that [latex]z_2[/latex] is a free variable. We have

[latex]\begin{align*} x = P \begin{pmatrix} R_1^{-1}Q_1^T(y - R_2z_2)\\ 0 \end{pmatrix} = x_0 + Lz_2, \end{align*}[/latex]

where

[latex]\begin{align*} L := -P\begin{pmatrix} R_1^{-1}Q_1^TR_2\\ 0. \end{pmatrix}. \end{align*}[/latex]

The set of solutions is the affine set [latex]\[ x_0 + \textbf{range}(L) \][/latex].