9 行 MATLAB 实现 Gram-Schmidt

1. 格拉姆-施密特算法 (Gram-Schmidt Algorithm)

格拉姆-施密特算法 (Gram-Schmidt algorithm) 将一组 $n$ 个独立向量 $\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_n$（即 $A$ 的列）变换为一组 $n$ 个标准正交向量 $\mathbf{q}_1, \mathbf{q}_2, \dots, \mathbf{q}_n$（即 $Q$ 的列）。

对于每一列 $j$：

1. 从 $\mathbf{a}_j$ 开始。
2. 减去它在所有先前已计算的 $\mathbf{q}_i$（其中 $i < j$）上的投影 (projection)。结果得到一个与所有先前 $\mathbf{q}_i$ 正交的向量 $\mathbf{v}$。
3. 用 $\mathbf{v}$ 的长度将其归一化 (normalize)，得到单位向量 $\mathbf{q}_j$。

内积 $\mathbf{q}_i^T \mathbf{a}_j$ 被收集到上三角矩阵 $R$ 中，满足 $A = QR$。$R$ 是上三角矩阵，因为当 $i > j$ 时 $\mathbf{q}_i^T \mathbf{a}_j = 0$（后来的 $\mathbf{q}$ 与先前的 $\mathbf{a}$ 正交——这正是该算法的全部意义所在）。

---

2. 完整的 9 行 MATLAB 代码 (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

循环结束后，$Q$ 包含标准正交列，$R$ 是上三角矩阵，满足 $A = QR$。

如果你"撤销"最后一步和中间步骤，则第 $j$ 列为：

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

将求和项移到左边，就得到矩阵乘法 $A = QR$ 的第 $j$ 列：

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

---

3. 修正的格拉姆-施密特 (Modified Gram-Schmidt)

代码第 4 行中从 $\mathbf{a}_j$ 到 $\mathbf{v}$ 的关键改变就是修正的格拉姆-施密特 (modified Gram-Schmidt)。

在精确算术 (exact arithmetic) 中，$R(i,j) = \mathbf{q}_i^T \mathbf{a}_j$ 与 $\mathbf{q}_i^T \mathbf{v}$ 相同（其中 $\mathbf{v}$ 已经减去了在先前 $\mathbf{q}_1, \dots, \mathbf{q}_{i-1}$ 上的投影）。这是因为新的 $\mathbf{q}_i$ 与那些先前方向正交。

然而，在浮点运算 (floating-point arithmetic) 中，正交性 (orthogonality) 并不完美。从当前的 $\mathbf{v}$（已经过部分正交化 (orthogonalized)）而不是从原始的 $\mathbf{a}_j$ 计算 $R(i,j)$，可以显著提高数值稳定性 (numerical stability)。所有人在该步骤中都使用 $\mathbf{v}$。这就是标准的修正格拉姆-施密特算法 (modified Gram-Schmidt algorithm)。

---

4. 一个 2x2 示例 (A 2x2 Example)

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

格拉姆-施密特过程将 $A$ 分解为 $QR$：

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

---

第 1 步：将 $\mathbf{a}_1$ 归一化为 $\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}$$

对角元素 $R(1,1) = \|\mathbf{a}_1\| = 5$。

---

第 2 步：减去 $\mathbf{a}_2$ 在 $\mathbf{q}_1$ 上的投影

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

首先计算内积（这将变为 $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$$

现在减去 $\mathbf{a}_2$ 在 $\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}$$

---

第 3 步：将 $\mathbf{v}$ 归一化为 $\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$$

这就是 $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. $R$ 中各项的解释 (Explanation of the Entries in $R$)

项 (Entry)	值 (Value)	含义 (Meaning)
$R(1,1) = 5$	$\|\mathbf{a}_1\|$	$A$ 第一列的长度。当我们将 $\mathbf{a}_1$ 归一化为 $\mathbf{q}_1$ 时，这个长度成为 $R$ 的第一个对角元素。
$R(2,2) = 2$	$\|\mathbf{v}\|$	减去 $\mathbf{a}_2$ 在 $\mathbf{q}_1$ 上的投影后，正交化向量 $\mathbf{v}$ 的长度。这是 $R$ 的第二个对角元素。
$R(1,2) = -1$	$\mathbf{q}_1^T \mathbf{a}_2$	$\mathbf{q}_1$ 与 $\mathbf{a}_2$ 之间的内积。这是投影系数 (projection coefficient)——$\mathbf{a}_2$ 在 $\mathbf{q}_1$ 方向上的分量大小。它出现在 $R$ 的上三角部分，因为在排序中 $\mathbf{q}_1$ 在 $\mathbf{q}_2$ 之前。
$R(2,1) = 0$	$\mathbf{q}_2^T \mathbf{a}_1$	零，因为根据构造 $\mathbf{q}_2$ 与 $\mathbf{a}_1$ 正交（等价地说，$\mathbf{a}_1$ 在 $\mathbf{q}_2$ 方向上没有分量）。这就是 $R$ 成为上三角矩阵的原因。