[proofplan] We prove by induction on $k$ that the Gram-Schmidt algorithm produces vectors $\mathbf{q}_1, \ldots, \mathbf{q}_k$ that are orthonormal and that $\operatorname{span}(\mathbf{q}_1, \ldots, \mathbf{q}_k) = \operatorname{span}(\mathbf{a}_1, \ldots, \mathbf{a}_k)$. At each inductive step, the residual $\mathbf{d}_k = \mathbf{a}_k - \sum_{i=1}^{k-1}(\mathbf{q}_i^\top \mathbf{a}_k)\mathbf{q}_i$ is orthogonal to $\mathbf{q}_1, \ldots, \mathbf{q}_{k-1}$ because it is the projection of $\mathbf{a}_k$ onto the orthogonal complement of $\operatorname{span}(\mathbf{q}_1, \ldots, \mathbf{q}_{k-1})$. Normalising $\mathbf{d}_k$ (or choosing an arbitrary orthogonal unit vector when $\mathbf{d}_k = \mathbf{0}$) extends the orthonormal system. The factorisation $A = QR$ with $R$ upper triangular then follows from the construction of the entries $R_{ik}$. [/proofplan]

[step:Establish the base case $k = 1$] At $k = 1$, the algorithm sets: - If $\mathbf{a}_1 \neq \mathbf{0}$: $\mathbf{q}_1 = \mathbf{a}_1 / \|\mathbf{a}_1\|$ and $R_{11} = \|\mathbf{a}_1\|$. Then $\|\mathbf{q}_1\| = 1$ and $\mathbf{a}_1 = R_{11} \mathbf{q}_1$, so $\operatorname{span}(\mathbf{q}_1) = \operatorname{span}(\mathbf{a}_1)$. - If $\mathbf{a}_1 = \mathbf{0}$: $R_{11} = 0$ and $\mathbf{q}_1$ is an arbitrary unit vector. Then $\|\mathbf{q}_1\| = 1$ and $\mathbf{a}_1 = \mathbf{0} = 0 \cdot \mathbf{q}_1 = R_{11}\mathbf{q}_1$. Moreover, $\operatorname{span}(\mathbf{a}_1) = \{\mathbf{0}\} \subset \operatorname{span}(\mathbf{q}_1)$. In both cases, the set $\{\mathbf{q}_1\}$ is orthonormal (a single unit vector) and $\mathbf{a}_1 = R_{11}\mathbf{q}_1$. [/step]

[step:State the inductive hypothesis] Assume that after step $k - 1$, the vectors $\mathbf{q}_1, \ldots, \mathbf{q}_{k-1}$ satisfy: 1. **Orthonormality:** $\mathbf{q}_i^\top \mathbf{q}_j = \delta_{ij}$ for all $i, j \in \{1, \ldots, k-1\}$. 2. **Span preservation:** $\operatorname{span}(\mathbf{q}_1, \ldots, \mathbf{q}_j) \supseteq \operatorname{span}(\mathbf{a}_1, \ldots, \mathbf{a}_j)$ for each $j = 1, \ldots, k-1$. 3. **Column reconstruction:** For each $j = 1, \ldots, k-1$, we have $\mathbf{a}_j = \sum_{i=1}^{j} R_{ij} \mathbf{q}_i$ (with $R_{ij} = \mathbf{q}_i^\top \mathbf{a}_j$ for $i < j$ and $R_{jj} = \|\mathbf{d}_j\|$). [/step]

[step:Verify orthogonality of $\mathbf{d}_k$ to the existing basis vectors]The residual at step $k$ is \begin{align*} \mathbf{d}_k = \mathbf{a}_k - \sum_{i=1}^{k-1} (\mathbf{q}_i^\top \mathbf{a}_k)\, \mathbf{q}_i. \end{align*} For any $\ell \in \{1, \ldots, k-1\}$, we compute \begin{align*} \mathbf{q}_\ell^\top \mathbf{d}_k = \mathbf{q}_\ell^\top \mathbf{a}_k - \sum_{i=1}^{k-1} (\mathbf{q}_i^\top \mathbf{a}_k)(\mathbf{q}_\ell^\top \mathbf{q}_i). \end{align*} By the inductive hypothesis (orthonormality), $\mathbf{q}_\ell^\top \mathbf{q}_i = \delta_{\ell i}$, so the sum collapses to a single term: \begin{align*} \mathbf{q}_\ell^\top \mathbf{d}_k = \mathbf{q}_\ell^\top \mathbf{a}_k - (\mathbf{q}_\ell^\top \mathbf{a}_k) \cdot 1 = 0. \end{align*} Therefore $\mathbf{d}_k \perp \mathbf{q}_\ell$ for all $\ell = 1, \ldots, k-1$.[/step]

[guided]The residual $\mathbf{d}_k$ is constructed by subtracting from $\mathbf{a}_k$ its orthogonal projection onto $\operatorname{span}(\mathbf{q}_1, \ldots, \mathbf{q}_{k-1})$. The projection of a vector $\mathbf{v}$ onto an orthonormal set $\{\mathbf{q}_1, \ldots, \mathbf{q}_{k-1}\}$ is $\sum_{i=1}^{k-1} (\mathbf{q}_i^\top \mathbf{v}) \mathbf{q}_i$, so the residual $\mathbf{d}_k = \mathbf{a}_k - \sum_{i=1}^{k-1} (\mathbf{q}_i^\top \mathbf{a}_k) \mathbf{q}_i$ is the component of $\mathbf{a}_k$ orthogonal to this subspace. To verify orthogonality rigorously: take the inner product of $\mathbf{d}_k$ with any $\mathbf{q}_\ell$ ($\ell \le k-1$): \begin{align*} \mathbf{q}_\ell^\top \mathbf{d}_k &= \mathbf{q}_\ell^\top \mathbf{a}_k - \sum_{i=1}^{k-1} (\mathbf{q}_i^\top \mathbf{a}_k) \underbrace{\mathbf{q}_\ell^\top \mathbf{q}_i}_{= \delta_{\ell i}} = \mathbf{q}_\ell^\top \mathbf{a}_k - \mathbf{q}_\ell^\top \mathbf{a}_k = 0. \end{align*} The orthonormality of $\{\mathbf{q}_1, \ldots, \mathbf{q}_{k-1}\}$ is consumed here: without it, the Kronecker delta replacement $\mathbf{q}_\ell^\top \mathbf{q}_i = \delta_{\ell i}$ would fail, and the sum would not telescope.[/guided]

[step:Normalise $\mathbf{d}_k$ and verify orthonormality of the extended set] **Case 1: $\mathbf{d}_k \neq \mathbf{0}$.** Set $\mathbf{q}_k = \mathbf{d}_k / \|\mathbf{d}_k\|$ and $R_{kk} = \|\mathbf{d}_k\|$. Then $\|\mathbf{q}_k\| = 1$, and for each $\ell \in \{1, \ldots, k-1\}$: \begin{align*} \mathbf{q}_\ell^\top \mathbf{q}_k = \frac{\mathbf{q}_\ell^\top \mathbf{d}_k}{\|\mathbf{d}_k\|} = \frac{0}{\|\mathbf{d}_k\|} = 0. \end{align*} Combined with the inductive hypothesis $\mathbf{q}_i^\top \mathbf{q}_j = \delta_{ij}$ for $i, j \le k-1$, the extended set $\{\mathbf{q}_1, \ldots, \mathbf{q}_k\}$ is orthonormal. **Case 2: $\mathbf{d}_k = \mathbf{0}$.** Set $R_{kk} = 0$ and choose $\mathbf{q}_k$ to be any unit vector orthogonal to $\mathbf{q}_1, \ldots, \mathbf{q}_{k-1}$. Such a vector exists because $k - 1 < m$ (the ambient space $\mathbb{R}^m$ has dimension $m \ge n \ge k$), so the orthogonal complement of $\operatorname{span}(\mathbf{q}_1, \ldots, \mathbf{q}_{k-1})$ is nontrivial. By construction, $\|\mathbf{q}_k\| = 1$ and $\mathbf{q}_\ell^\top \mathbf{q}_k = 0$ for $\ell = 1, \ldots, k-1$, so $\{\mathbf{q}_1, \ldots, \mathbf{q}_k\}$ is orthonormal. [/step]

[step:Verify span preservation at step $k$]We show $\operatorname{span}(\mathbf{q}_1, \ldots, \mathbf{q}_k) \supseteq \operatorname{span}(\mathbf{a}_1, \ldots, \mathbf{a}_k)$. By the inductive hypothesis, $\operatorname{span}(\mathbf{q}_1, \ldots, \mathbf{q}_{k-1}) \supseteq \operatorname{span}(\mathbf{a}_1, \ldots, \mathbf{a}_{k-1})$, so it suffices to show $\mathbf{a}_k \in \operatorname{span}(\mathbf{q}_1, \ldots, \mathbf{q}_k)$. From the definition of $\mathbf{d}_k$: \begin{align*} \mathbf{a}_k = \sum_{i=1}^{k-1} (\mathbf{q}_i^\top \mathbf{a}_k)\, \mathbf{q}_i + \mathbf{d}_k. \end{align*} **Case 1: $\mathbf{d}_k \neq \mathbf{0}$.** Then $\mathbf{d}_k = R_{kk}\, \mathbf{q}_k$ with $R_{kk} = \|\mathbf{d}_k\| > 0$. Substituting: \begin{align*} \mathbf{a}_k = \sum_{i=1}^{k-1} R_{ik}\, \mathbf{q}_i + R_{kk}\, \mathbf{q}_k = \sum_{i=1}^{k} R_{ik}\, \mathbf{q}_i, \end{align*} where $R_{ik} = \mathbf{q}_i^\top \mathbf{a}_k$ for $i < k$. This lies in $\operatorname{span}(\mathbf{q}_1, \ldots, \mathbf{q}_k)$. **Case 2: $\mathbf{d}_k = \mathbf{0}$.** Then $\mathbf{a}_k = \sum_{i=1}^{k-1} (\mathbf{q}_i^\top \mathbf{a}_k)\, \mathbf{q}_i \in \operatorname{span}(\mathbf{q}_1, \ldots, \mathbf{q}_{k-1}) \subset \operatorname{span}(\mathbf{q}_1, \ldots, \mathbf{q}_k)$. In this case $R_{kk} = 0$, so $\mathbf{a}_k = \sum_{i=1}^{k} R_{ik}\, \mathbf{q}_i$ still holds (the $k$-th term is zero). In both cases, $\mathbf{a}_k \in \operatorname{span}(\mathbf{q}_1, \ldots, \mathbf{q}_k)$ and the column reconstruction formula $\mathbf{a}_k = \sum_{i=1}^{k} R_{ik}\, \mathbf{q}_i$ is established.[/step]

[guided]The span preservation property has two directions. The inclusion $\operatorname{span}(\mathbf{a}_1, \ldots, \mathbf{a}_k) \subseteq \operatorname{span}(\mathbf{q}_1, \ldots, \mathbf{q}_k)$ is what we just proved. Does the reverse inclusion hold? Not necessarily. When $\mathbf{d}_k = \mathbf{0}$, the vector $\mathbf{a}_k$ already lies in $\operatorname{span}(\mathbf{q}_1, \ldots, \mathbf{q}_{k-1})$, so adding $\mathbf{a}_k$ does not enlarge the span. But the algorithm chooses $\mathbf{q}_k$ orthogonal to $\mathbf{q}_1, \ldots, \mathbf{q}_{k-1}$, introducing a genuinely new direction. Thus $\operatorname{span}(\mathbf{q}_1, \ldots, \mathbf{q}_k)$ may be strictly larger than $\operatorname{span}(\mathbf{a}_1, \ldots, \mathbf{a}_k)$ when columns of $A$ are linearly dependent. This is not a defect of the algorithm: the purpose of the $\mathbf{d}_k = \mathbf{0}$ branch is to ensure that $Q$ always has orthonormal columns, even when $A$ is rank-deficient. The column $\mathbf{q}_k$ is a "filler" that completes the orthonormal system without affecting the factorisation $A = QR$ (because $R_{kk} = 0$ zeroes out the contribution of $\mathbf{q}_k$ in the product $QR$).[/guided]

[step:Assemble the factorisation $A = QR$ with $R$ upper triangular]After completing all $n$ steps, define the matrices: \begin{align*} Q := \begin{pmatrix} \mathbf{q}_1 & \mathbf{q}_2 & \cdots & \mathbf{q}_n \end{pmatrix} \in \mathbb{R}^{m \times n}, \qquad R := (R_{ik}) \in \mathbb{R}^{n \times n}. \end{align*} The column reconstruction formula from the inductive argument gives, for each $k = 1, \ldots, n$: \begin{align*} \mathbf{a}_k = \sum_{i=1}^{k} R_{ik}\, \mathbf{q}_i. \end{align*} This means $R_{ik} = 0$ for $i > k$ (since the sum only runs to $i = k$), so $R$ is upper triangular. In matrix form, the equation $\mathbf{a}_k = Q \mathbf{r}_k$ (where $\mathbf{r}_k$ is the $k$-th column of $R$) for each $k = 1, \ldots, n$ assembles into \begin{align*} A = QR. \end{align*} The columns of $Q$ form an orthonormal set by the induction: $Q^\top Q = I_n$. The matrix $R$ is upper triangular with diagonal entries $R_{kk} = \|\mathbf{d}_k\| \ge 0$ and off-diagonal entries $R_{ik} = \mathbf{q}_i^\top \mathbf{a}_k$ for $i < k$.[/step]

[guided]Let us verify the matrix equation $A = QR$ entry by entry. The $(p, k)$ entry of $QR$ is \begin{align*} (QR)_{pk} = \sum_{i=1}^{n} Q_{pi} R_{ik} = \sum_{i=1}^{n} (\mathbf{q}_i)_p \, R_{ik}. \end{align*} Since $R$ is upper triangular, $R_{ik} = 0$ for $i > k$, so the sum truncates to $i = 1, \ldots, k$: \begin{align*} (QR)_{pk} = \sum_{i=1}^{k} (\mathbf{q}_i)_p \, R_{ik} = \left(\sum_{i=1}^{k} R_{ik}\, \mathbf{q}_i\right)_p = (\mathbf{a}_k)_p = A_{pk}, \end{align*} which is precisely the column reconstruction formula. This confirms $A = QR$. The upper triangularity of $R$ is a consequence of the sequential nature of the algorithm: when processing column $\mathbf{a}_k$, only the previously constructed vectors $\mathbf{q}_1, \ldots, \mathbf{q}_k$ are used, so $\mathbf{a}_k$ is expressed as a linear combination of $\mathbf{q}_1, \ldots, \mathbf{q}_k$ with no contribution from $\mathbf{q}_{k+1}, \ldots, \mathbf{q}_n$. This forces $R_{ik} = 0$ for $i > k$.[/guided]