Gram-Schmidt Orthogonalization -- MIT 18.06 Notes

1. The Gram-Schmidt Algorithm

The Gram-Schmidt algorithm takes a set of $n$ independent vectors $\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_n$ (the columns of $A$) and produces a set of $n$ orthonormal vectors $\mathbf{q}_1, \mathbf{q}_2, \dots, \mathbf{q}_n$ (the columns of $Q$).

For each column $j$:

1. Start with $\mathbf{a}_j$.
2. Subtract off its projections onto all previously computed $\mathbf{q}_i$ (for $i < j$). The result is a vector $\mathbf{v}$ that is orthogonal to all earlier $\mathbf{q}_i$.
3. Normalize $\mathbf{v}$ by its length to produce the unit vector $\mathbf{q}_j$.

The inner products $\mathbf{q}_i^T \mathbf{a}_j$ are collected into the upper-triangular matrix $R$, satisfying $A = QR$. $R$ is upper triangular because $\mathbf{q}_i^T \mathbf{a}_j = 0$ when $i > j$ (later $\mathbf{q}$'s are orthogonal to earlier $\mathbf{a}$'s -- that is the whole point of the algorithm).

---

2. The Complete 9-Line MATLAB Code

[m,n] = size(A);
Q = zeros(m,n); R = zeros(n,n);
for j = 1:n
    % Gram-Schmidt orthogonalization
    v = A(:,j);                 % begins as column j of A
    for i = 1:j-1
        R(i,j) = Q(:,i)' * A(:,j);
        % modify A(j) to v for more accuracy
        v = v - R(i,j) * Q(:,i);
        % subtract the projection (qi^T aj) qi
    end
    % v is now perpendicular to all of q1,...,qj-1
    R(j,j) = norm(v);
    Q(:,j) = v / R(j,j);        % normalize v to be the next unit vector qj
end

After the loop, $Q$ contains orthonormal columns and $R$ is upper triangular such that $A = QR$.

If you "undo" the last step and the intermediate steps, column $j$ is:

$$R(j,j)\,\mathbf{q}_j = \mathbf{a}_j - \sum_{i=1}^{j-1} R(i,j)\,\mathbf{q}_i$$

Moving the sum to the left gives column $j$ of the matrix multiplication $A = QR$:

$$\mathbf{a}_j = \sum_{i=1}^{j} R(i,j)\,\mathbf{q}_i$$

---

3. Modified Gram-Schmidt

The crucial change from $\mathbf{a}_j$ to $\mathbf{v}$ in line 4 of the code gives modified Gram-Schmidt.

In exact arithmetic, $R(i,j) = \mathbf{q}_i^T \mathbf{a}_j$ is the same as $\mathbf{q}_i^T \mathbf{v}$ (where $\mathbf{v}$ has already had its projections onto earlier $\mathbf{q}_1, \dots, \mathbf{q}_{i-1}$ subtracted). This is because the new $\mathbf{q}_i$ is orthogonal to those earlier directions.

In floating-point arithmetic, however, orthogonality is not perfect. Computing $R(i,j)$ from the current $\mathbf{v}$ (which has already been partially orthogonalized) rather than from the original $\mathbf{a}_j$ gives significantly better numerical stability. Everybody uses $\mathbf{v}$ at that step in the code. This is the standard modified Gram-Schmidt algorithm.

---

4. A 2x2 Example

$$A = \begin{bmatrix} 4 & -2 \\ 3 & 1 \end{bmatrix}$$

The Gram-Schmidt process factors $A$ into $QR$:

$$A = QR = \frac{1}{5}\begin{bmatrix} 4 & -3 \\ 3 & 4 \end{bmatrix} \begin{bmatrix} 5 & -1 \\ 0 & 2 \end{bmatrix}$$

---

Step 1: Normalize $\mathbf{a}_1$ to $\mathbf{q}_1$

$$\mathbf{a}_1 = \begin{bmatrix} 4 \\ 3 \end{bmatrix}, \qquad \|\mathbf{a}_1\| = \sqrt{4^2 + 3^2} = \sqrt{25} = 5$$ $$\mathbf{q}_1 = \frac{\mathbf{a}_1}{\|\mathbf{a}_1\|} = \frac{1}{5}\begin{bmatrix} 4 \\ 3 \end{bmatrix}$$

The diagonal entry $R(1,1) = \|\mathbf{a}_1\| = 5$.

---

Step 2: Subtract the projection of $\mathbf{a}_2$ onto $\mathbf{q}_1$

$$\mathbf{a}_2 = \begin{bmatrix} -2 \\ 1 \end{bmatrix}$$

First compute the inner product (this becomes $R(1,2)$):

$$R(1,2) = \mathbf{q}_1^T \mathbf{a}_2 = \frac{1}{5}(4)(-2) + \frac{1}{5}(3)(1) = -\frac{8}{5} + \frac{3}{5} = -1$$

Now subtract the projection of $\mathbf{a}_2$ onto $\mathbf{q}_1$:

$$\mathbf{v} = \mathbf{a}_2 - R(1,2)\,\mathbf{q}_1 = \begin{bmatrix} -2 \\ 1 \end{bmatrix} - (-1)\frac{1}{5}\begin{bmatrix} 4 \\ 3 \end{bmatrix} = \begin{bmatrix} -2 \\ 1 \end{bmatrix} + \begin{bmatrix} 4/5 \\ 3/5 \end{bmatrix} = \begin{bmatrix} -6/5 \\ 8/5 \end{bmatrix} = \frac{2}{5}\begin{bmatrix} -3 \\ 4 \end{bmatrix}$$

---

Step 3: Normalize $\mathbf{v}$ to $\mathbf{q}_2$

$$\|\mathbf{v}\| = \sqrt{\left(-\frac{6}{5}\right)^2 + \left(\frac{8}{5}\right)^2} = \sqrt{\frac{36}{25} + \frac{64}{25}} = \sqrt{\frac{100}{25}} = \sqrt{4} = 2$$

This is $R(2,2) = 2$.

$$\mathbf{q}_2 = \frac{\mathbf{v}}{\|\mathbf{v}\|} = \frac{1}{2}\begin{bmatrix} -6/5 \\ 8/5 \end{bmatrix} = \frac{1}{5}\begin{bmatrix} -3 \\ 4 \end{bmatrix}$$

---

Verification

$$Q = \begin{bmatrix} 4/5 & -3/5 \\ 3/5 & 4/5 \end{bmatrix}, \qquad R = \begin{bmatrix} 5 & -1 \\ 0 & 2 \end{bmatrix}$$ $$QR = \frac{1}{5}\begin{bmatrix} 4 & -3 \\ 3 & 4 \end{bmatrix} \begin{bmatrix} 5 & -1 \\ 0 & 2 \end{bmatrix} = \frac{1}{5}\begin{bmatrix} 4(5)+(-3)(0) & 4(-1)+(-3)(2) \\ 3(5)+4(0) & 3(-1)+4(2) \end{bmatrix} = \frac{1}{5}\begin{bmatrix} 20 & -10 \\ 15 & 5 \end{bmatrix} = \begin{bmatrix} 4 & -2 \\ 3 & 1 \end{bmatrix} = A$$

---

5. Explanation of the Entries in $R$

Entry	Value	Meaning
$R(1,1) = 5$	$\|\mathbf{a}_1\|$	The length of the first column of $A$. When we normalize $\mathbf{a}_1$ to get $\mathbf{q}_1$, this length becomes the first diagonal entry of $R$.
$R(2,2) = 2$	$\|\mathbf{v}\|$	The length of the orthogonalized vector $\mathbf{v}$ after subtracting the projection of $\mathbf{a}_2$ onto $\mathbf{q}_1$. This is the second diagonal entry of $R$.
$R(1,2) = -1$	$\mathbf{q}_1^T \mathbf{a}_2$	The inner product between $\mathbf{q}_1$ and $\mathbf{a}_2$. This is the projection coefficient -- how much of $\mathbf{a}_2$ lies in the direction of $\mathbf{q}_1$. It appears in the upper-triangular part of $R$ because $\mathbf{q}_1$ comes before $\mathbf{q}_2$ in the ordering.
$R(2,1) = 0$	$\mathbf{q}_2^T \mathbf{a}_1$	Zero because $\mathbf{q}_2$ is orthogonal to $\mathbf{a}_1$ by construction (equivalently, $\mathbf{a}_1$ has no component in the $\mathbf{q}_2$ direction). This is what makes $R$ upper triangular.