Alongside SVD and LU decomposition, another highly useful way to factorize a matrix is the QR decomposition. For an $m \times n$ matrix $A$ (often with $m \ge n$ and linearly independent columns), the QR decomposition finds two specific matrices, $Q$ and $R$ , such that:

$A = QR$

Here, $Q$ is an $m \times n$ matrix with orthonormal columns, and $R$ is an $n \times n$ upper triangular matrix. Let's examine these components more closely.

The Orthogonal Matrix Q

The matrix $Q$ has a special property: its columns form an orthonormal set. Recall from Chapter 4 that this means two things:

Each column vector has a Euclidean norm (length) of 1 ( $||\mathbf{q}_i||_2 = 1$ ).
Any two distinct column vectors are orthogonal to each other (their dot product is zero, $\mathbf{q}_i \cdot \mathbf{q}_j = \mathbf{q}_i^T \mathbf{q}_j = 0$ for $i \neq j$ ).

Combining these, we can state the defining property of a matrix $Q$ with orthonormal columns concisely using matrix multiplication:

$Q^T Q = I$

where $I$ is the $n \times n$ identity matrix.

Matrices satisfying this condition are called orthogonal matrices. (Technically, if $Q$ were square, $Q^T Q = I$ would also imply $Q Q^T = I$ , making $Q^T$ the inverse of $Q$ , $Q^T = Q^{-1}$ ).

A significant geometric implication of orthogonal matrices is that multiplication by $Q$ (or $Q^T$ ) preserves lengths and angles. Transformations represented by orthogonal matrices are rigid motions, like rotations or reflections. In the context of QR decomposition, the columns of $Q$ form an orthonormal basis for the column space of the original matrix $A$ . This process is conceptually similar to the Gram-Schmidt procedure, which takes a set of linearly independent vectors (the columns of $A$ ) and produces an orthonormal set ( $Q$ ) that spans the same space.

The Upper Triangular Matrix R

The matrix $R$ is an $n \times n$ upper triangular matrix. This means all entries below the main diagonal are zero:

R = \begin{pmatrix} r_{11} & r_{12} & \cdots & r_{1n} \\ 0 & r_{22} & \cdots & r_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & r_{nn} \end{pmatrix}

The matrix $R$ essentially encodes the information about how the original columns of $A$ are related to the orthonormal basis vectors in $Q$ . Specifically, each column $\mathbf{a}_j$ of $A$ can be expressed as a linear combination of the first $j$ columns of $Q$ , with the coefficients coming from the $j$ -th column of $R$ .

The upper triangular structure of $R$ is computationally advantageous, particularly when solving systems of linear equations, as we'll see next.

The QR decomposition factorizes matrix A into an orthogonal matrix Q (columns are orthonormal vectors) and an upper triangular matrix R.

Applications in Least Squares Problems

One of the most important applications of QR decomposition is in solving linear least-squares problems. Recall that these problems often arise in machine learning when fitting models, like linear regression, where we want to find a vector $\mathbf{x}$ that minimizes the difference between $A\mathbf{x}$ and a target vector $\mathbf{b}$ , i.e., minimize $||A\mathbf{x} - \mathbf{b}||_2^2$ .

The standard approach involves solving the normal equations:

$A^T A \mathbf{x} = A^T \mathbf{b}$

However, computing $A^T A$ can sometimes lead to numerical stability issues, especially if the columns of $A$ are close to being linearly dependent. The matrix $A^T A$ might become ill-conditioned (have a large condition number), making the solution sensitive to small errors.

QR decomposition provides a more numerically stable alternative. If we substitute $A = QR$ into the original least-squares objective (or directly into $A\mathbf{x} = \mathbf{b}$ for systems where a solution exists), we get:

$QR\mathbf{x} = \mathbf{b}$

Now, multiply both sides by $Q^T$ :

$Q^T Q R \mathbf{x} = Q^T \mathbf{b}$

Since $Q^T Q = I$ , this simplifies beautifully to:

$R \mathbf{x} = Q^T \mathbf{b}$

This gives us a system of linear equations with an upper triangular matrix $R$ . Such systems are very easy and numerically stable to solve using a process called back-substitution. We first solve for the last component of $\mathbf{x}$ , then substitute that value back into the second-to-last equation to solve for the next component, and so on, working our way up.

Because this method avoids forming $A^T A$ , it generally has better numerical properties than using the normal equations directly, making it the preferred method in many numerical libraries.

Computing the QR Decomposition

While the conceptual idea behind QR decomposition can be understood through the Gram-Schmidt process, practical implementations typically use more numerically stable algorithms like Householder reflections or Givens rotations. These methods construct $Q$ and $R$ iteratively. Fortunately, you rarely need to implement these yourself; numerical libraries like NumPy and SciPy provide efficient and reliable functions to compute the QR decomposition, which we will utilize in the hands-on sections.