Gram-Schmidt process ($u_j=v_j-{\sum }_{i<j}⟨v_j,q_i⟩q_i,q_j=u_j/\Vert u_j\Vert$)

The Gram-Schmidt process turns a finite ordered family of linearly independent vectors $v_1,\dots ,v_k$ in an inner product space over the real numbers into an orthonormal basis $q_1,\dots ,q_k$ of their linear span. It keeps every successive linear span, so the first $j$ new vectors describe exactly the same directions as the first $j$ inputs. Its idea is geometric: remove the orthogonal projection onto directions already found, then divide the remaining vector by its norm. This makes a complicated basis easier to use without changing the vector subspace it describes.^[ref]

Start with $u_1=v_1$ and $q_1=v_1/\Vert v_1\Vert$. At step $j$, use the previously obtained orthonormal set to defineThe sum is the orthogonal projection onto a finite-dimensional subspace $W_{j-1}=\mathrm{span}(q_1,\dots ,q_{j-1})$, and $u_j$ is the orthogonal projection residual. Linear independence ensures $u_j\ne 0$: otherwise $v_j$ would already be a linear combination of earlier inputs. If the input family is linearly dependent, use only the retained orthonormal vectors in the sum and omit an input whenever its remainder is the zero vector; do not divide by zero. The retained outputs still form an orthonormal basis of the original linear span, but their number can be smaller than $k$.

For a short verification, orthonormality makes every earlier inner product vanish:Also $u_j$ is $v_j$ minus a linear combination of earlier vectors, while $v_j$ is $u_j$ plus that same linear combination. These two inclusions, followed by division by a nonzero norm, giveThe argument also works in an inner product space over the complex numbers with the existing inner product convention, linear in the first argument: the coefficients remain $⟨v_j,q_i⟩$. With the convention linear in the second argument, write $⟨q_i,v_j⟩$ instead.

For a first picture, take $v_1=(1,1{)}^T$ and $v_2=(1,0{)}^T$ in ${\mathbb{R}}^2$ with the ordinary dot product. ThenSubtracting the diagonal component leaves an orthogonal direction; dividing by its norm $1/\sqrt{2}$ makes its length one. If the second input were $(2,2{)}^T$ instead, its orthogonal projection would be the entire input and $u_2=0$, revealing linear dependence.

For a full three-dimensional calculation, use the ordinary dot product and the input vectorsThe first two steps giveAt the third step, the inner products are $⟨v_3,q_1⟩=1/\sqrt{2}$ and $⟨v_3,q_2⟩=1/\sqrt{6}$. ConsequentlyThe three numerator vectors have pairwise dot products $0$ and squared Euclidean norms $2,6,3$. Thus the displayed outputs are indeed an orthonormal basis, and the three nonzero remainders also confirm linear independence of the inputs. Keeping the removed coefficients produces the QR decomposition below.

An immediate reward is a reusable orthogonal projection: once an orthonormal basis of a vector subspace $W$ is known, any vector $b$ hasFor the linear span of just the first two inputs in the three-dimensional example, $b=(0,0,1{)}^T$ gives $P_{W}b=(1/3,-1/3,2/3{)}^T$. Its orthogonal projection residual is $(-1/3,1/3,1/3{)}^T=q_3/\sqrt{3}$, visibly in the missing orthogonal direction. This is the geometry behind the linear least-squares problem.

The formulas above describe exact arithmetic. With floating-point arithmetic, subtracting an orthogonal projection nearly equal to the input can leave a small remainder; dividing by its small norm can magnify existing roundoff error. The Classical Gram-Schmidt and Modified Gram-Schmidt arrangements agree in exact arithmetic but can give different rounded results. Exact linear dependence and a small computed remainder therefore require different judgments in a numerical implementation.^[ref]

For a linearly independent ordered family in a real inner product space, require each unit vector $w_j$ to lie in the span of the first $j$ inputs and be orthogonal to all earlier inputs. That orthogonal part of the successive span is one-dimensional. Subtracting the earlier orthonormal components gives its nonzero generator $r_j$. Every possible unit choice is therefore $\pm r_j/\Vert r_j\Vert$. Since $⟨r_j,u_j⟩=\Vert r_j{\Vert }^2>0$, the positive normalization fixes the sign uniquely. This distinguishes uniqueness of the normalized ordered construction from the many possible orthonormal bases of the whole space.

Classical Gram-Schmidt computes all orthogonal projection coefficients $⟨v_j,q_i⟩$ from the original input $v_j$, then subtracts their sum as in the displayed Gram-Schmidt process formula.

Modified Gram-Schmidt starts with $w=v_j$ and, for each retained orthonormal vector $q_i$ in order, replaces $w$ by $w-⟨w,q_i⟩q_i$. Each orthogonal projection coefficient uses the current remainder. In exact arithmetic, subtracting earlier orthogonal directions does not change the inner product with $q_i$, so the final $w$ equals the Classical Gram-Schmidt remainder $u_j$; then divide by $\Vert w\Vert$ if it is nonzero.

A reduced QR decomposition of an $m\times k$ matrix over the real numbers with full column rank, $m\ge k$, writes $A=QR$, where $Q$ is $m\times k$ with orthonormal columns and $R$ is a $k\times k$ upper triangular matrix with positive diagonal. Applying the Gram-Schmidt process to the columns $v_j$ of $A$ constructs $Q=(q_1\cdots q_k)$ and records the coefficientsThis last column identity is exactly the matrix multiplication $A=QR$. The upper triangular matrix records that the $j$th input uses only the first $j$ new directions, while its positive diagonal records the norms of the nonzero orthogonal projection residuals. Since $Q^{T}Q=I$, one also has $R=Q^{T}A$; $Q^T$ denotes the transpose and $I$ the identity matrix. Over the complex numbers, use the conjugate transpose $Q^{\ast }$ in these identities.

For the three-dimensional Gram-Schmidt process example, the exact matrices areFor example, the last column reconstructs asThe diagonal entries are approximately $1.414214,1.224745,1.154701$. Each is a norm of the new direction before it was scaled to length one.

For a linear least-squares problem, the orthogonal projection of $b$ onto the column space of $A$ is $Q{Q}^{T}b$. Hence the minimizing coefficients solve the linear system $Rx=Q^{T}b$ by back substitution. Let $A_2=(v_1{v}_2)$ consist of the first two columns of the example; its reduced factors are $Q_2=(q_1{q}_2)$ and the leading $2\times 2$ block $R_2$ of $R$. For $b=(0,0,1{)}^T$, the linear system $R_{2}x=Q_2^{T}b$ isThus $A_{2}x=(1/3,-1/3,2/3{)}^T$, agreeing with the orthogonal projection computed above, and the minimum Euclidean norm of $A_{2}x-b$ is $1/\sqrt{3}$. The QR decomposition separates building an orthonormal basis from recovering coefficients in the original basis.^[ref]

Let $A=QR$ with $Q^{T}Q=I$ and square invertible upper triangular $R$. The normal equations for linear least squares reduce from $R^{T}Rx=R^T{Q}^{T}b$ to the system in the title. It is solved by back substitution. The geometric identity $\Vert Ax-b{\Vert }_2^2=\Vert Rx-Q^{T}b{\Vert }_2^2+\Vert (I-Q{Q}^T)b{\Vert }_2^2$ proves existence and uniqueness of the minimizer. For tall $Q$, $Q{Q}^T$ is an orthogonal projection, not the identity on the entire observation space.

Let $K$ be block diagonal with $2\times 2$ blocks $J=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)$. If a real invertible matrix $A$ commutes with $K$, its QR decomposition $A=BC$, with $B$ orthogonal and $C$ upper triangular with positive diagonal, has both $BK=KB$ and $CK=KC$. Indeed $CK{C}^{-1}=B^{-1}KB$, so $E=K^{-1}CK{C}^{-1}$ is orthogonal and upper triangular in $2\times 2$ blocks. Orthogonality forces its off-diagonal blocks to vanish. Its diagonal blocks are $J^{-1}{C}_{i}J{C}_i^{-1}$. For $C_i=\left(\begin{array}{cc}a& b\\ 0& d\end{array}\right)$ with $a,d>0$, this block is $\left(\begin{array}{cc}d/a& -b/a\\ -b/a& (a^2+b^2)/(ad)\end{array}\right)$. Orthogonality of its columns forces $b=0,d=a$, and hence each block is the identity. Thus $E=I$, proving the claim.

If a real $m$-by-$n$ matrix $A$ has full column rank, with $m\ge n$, then $A^{T}A$ is positive-definite. Its upper-triangular Cholesky decomposition $A^{T}A=R^{T}R$, with positive diagonal, is unique. Defining $Q=A{R}^{-1}$ gives $Q^{T}Q=I$ and $A=QR$. Any other such factorization yields another Cholesky factor of $A^{T}A$, so its $R$ and then its $Q$ must agree. This proves existence and uniqueness; it does not assert that forming $A^{T}A$ is the most numerically stable algorithm.

For a full-column-rank $n\times k$ matrix $V$, extend the thin QR decomposition $V=Q_1{R}_1$ to an orthogonal matrix $Q=(Q_1{Q}_2)$. Thenso $Q^T$ maps the column space of $V$ onto the first $k$ coordinate directions. Successive Householder transformations construct the same reduction without first forming a full orthogonal basis.

A linear least-squares problem minimizes $\Vert Ax-b{\Vert }_2$ over $x$. If $A=QR$ has full column rank, orthogonality reduces it to the upper-triangular system formed from the leading block of $R$ and the leading entries of $Q^{T}b$.

For a symmetric bilinear form $B$, the same projection formulaproduces a $B$-orthogonal basis whenever the chosen pivots $B(u_j,u_j)$ are nonzero. Reordering the input vectors can avoid a zero pivot when a suitable anisotropic vector remains.