Gram-Schmidt Orthonormalization and QR Decomposition of a Matrix

Theorem 4.

In every finite dimensional unitary (Euclidean) vector space an orthonormal basis exists.

Proof of this theorem, which provides an effective method of construction of such a basis, poceeds in two parallel stages. One starts with an arbitrary basis, applies the Gram-Schmidt orthogonalization method to construct an orthogonal basis, and then by normalisation obtains an orthonormal basis.

Construction of the Orthonormal Basis

Let \(\,\dim\,V=n.\ \) Proof proceeds according to the following diagram:

\[\begin{split}\begin{array}{cccccc} & \mathcal{B}=(v_1,v_2,\dots,v_n) & \longrightarrow & \mathcal{P}=(y_1,y_2,\dots,y_n) & \longrightarrow & \mathcal{Q}=(u_1,u_2,\dots,u_n) \\ \text{Basis:} & \text{initial} & & \text{orthogonal} & & \text{orthonormal} \end{array}\end{split}\]

We choose vectors for bases \(\ \mathcal{P}\ \) and \(\ \mathcal{Q}\ \) as follows:

\(\triangleright\quad y_1\,=\,v_1,\quad u_1\,=\,\displaystyle\frac{y_1}{\|y_1\|}\,;\)

\(\triangleright\quad y_2\,=\,\alpha_1\,y_1+\,v_2\ =\ v_2\,-\ \displaystyle\frac{\langle\,y_1,v_2\rangle}{\|y_1\|^2}\ \ y_1\ =\ v_2-\,\langle\,u_1,v_2\rangle\ u_1\,,\qquad u_2\,=\,\displaystyle\frac{y_2}{\|y_2\|}\ ,\)

where \(\ \alpha_1\ \) has been determined from the condition \(\ \langle\,y_1,y_2\rangle=0\,;\)

\(\triangleright\quad y_3\,=\,\beta_1\,y_1+\,\beta_2\,y_2+\,v_3\ =\ v_3\,-\ \displaystyle\frac{\langle y_1,v_3\rangle}{\|y_1\|^2}\ \ y_1\,-\ \displaystyle\frac{\langle y_2,v_3\rangle}{\|y_2\|^2}\ \ y_2\ =\ v_3\,-\,\langle u_1,v_3\rangle\ u_1 - \langle u_2,v_3\rangle\ u_2\,,\)

where \(\ \beta_1,\,\beta_2\ \) have been determined from the conditions \(\,\langle\,y_1,y_3\rangle = \langle\,y_2,y_3\rangle = 0\,;\qquad u_3\,=\,\displaystyle\frac{y_3}{\|y_3\|}\,;\)

Assume that in this way we have constructed an orthogonal set \(\ (y_1,y_2,\dots,y_{j-1}),\ \) where, as follows from the construction,

(1)\[y_k\,\in\,L(v_1,v_2,\dots,v_k)\qquad k=1,2,\dots,j-1\,,\]

that is, every \(\,k\)-th vector of the set \(\,(y_1,y_2,\dots,y_{j-1})\ \) is a linear combination of only first \(\,k\ \) vectors from the basis \(\,\mathcal{B}.\ \) Next vectors of bases \(\ \mathcal{P}\ \) and \(\ \mathcal{Q}\ \,\) are \(\\\)

\(\begin{array}{rcl} \triangleright\quad y_j & = & \lambda_1\,y_1\,+\;\lambda_2\,y_2\,+\;\dots\,+\;\lambda_{j-1}\,y_{j-1}\,+\;v_j\ \ = \\ \\ & = & v_j\ -\ \, \displaystyle\frac{\langle\,y_1,v_j\rangle}{\|y_1\|^2}\ \ y_1\ \,-\ \, \displaystyle\frac{\langle\,y_2,v_j\rangle}{\|y_2\|^2}\ \ y_2\ \,-\ \, \ldots\ -\ \, \displaystyle\frac{\langle\,y_{j-1},v_j\rangle}{\|y_{j-1}\|^2}\ \ y_{j-1}\ \ = \\ \\ & = & v_j\,-\,\langle\,u_1,v_j\rangle\ u_1 - \langle\,u_2,v_j\rangle\ u_2 \ -\ \ldots\ -\ \langle\,u_{j-1},v_j\rangle\ u_{j-1}\,, \qquad u_j\ =\ \displaystyle\frac{y_j}{\|y_j\|}\ , \end{array}\)

where the coefficients \(\ \lambda_1,\,\lambda_2,\,\dots,\,\lambda_{j-1}\ \,\) have been determined from the orthogonality conditions

\[\langle\,y_1,y_j\rangle\ \ =\ \ \langle\,y_2,y_j\rangle\ \ =\ \ldots\ \ =\ \ \langle\,y_{j-1},y_j\rangle\ \ =\ \ 0\,.\]

The vector \(\,y_j\ \) is a combination of linearly independent vectors \(\,v_1,\,v_2,\,\dots,\,v_j.\ \) The property (1) implies that it is a non-trivial combination, because the coefficient of the vector \(\,v_j\ \) equals \(\,+1.\,\) This means that all the vectors \(\,y_j\ \) are non-zero: \(\,y_j\neq\theta,\ j=1,2,\dots,n.\ \) Hence, the set \(\,\mathcal{P}\ \) is orthogonal, and if we divide each vector \(\,y_j\ \) by its norm, we obtain a unit vector \(\,u_j\,.\)

Gram-Schmidt orthogonalization may be applied to any (not necessarily linearly independent set) of vectors. As a result we obtain a set \(\,\mathcal{Y},\ \) whose vectors are pairwise orthogonal, but \(\,\mathcal{Y}\ \) itself is not necessarily orthogonal as it may contain zero vectors.

After the last step of the procedure, where \(\,j=n,\ \) we obtain an orthogonal basis \(\ \mathcal{P}\ \\\) and an orthonormal basis \(\ \mathcal{Q}\ \) of the space \(\,V.\)

Example: Construction of Legendre polynomials.

Consider a Euclidean space \(\,P\ \) of real polynomials with an inner product

\[\langle\,p,q\,\rangle\ \,=\ \,\int_{-1}^{+1}\ p(x)\,q(x)\,dx\,,\qquad p,q\in P\,.\]

Taking \(\,\mathcal{B}\,=\,(1,\,x,\,x^2,\,x^3,\,\dots)\ \) as the initial basis, we will construct an orthogonal basis \(\,\mathcal{P}\,=\,(p_0,\,p_1,\,p_2,\,p_3,\,\dots)\ \,\) and \(\,\) an orthonormal basis \(\,\mathcal{Q}\,=\,(q_0,\,q_1,\,q_2,\,q_3,\,\dots)\,.\\\)

0. \(\ p_0(x)\,=\,1\,;\qquad \|\,p_0\|^2\,=\,\langle\,p_0,p_0\,\rangle\ =\ \int_{-1}^{+1}\ 1\cdot 1\ \ dx\ =\ 2\,;\qquad \|\,p_0\|\,=\,\sqrt{2}\,.\)

   \(\blacktriangleright\quad q_0(x)\ \,=\ \, \displaystyle\frac{1}{\|\,p_0\|}\ \ p_0(x)\ \,=\ \, \frac{1}{\sqrt{2}}\ \ .\\\)

1. \(\ p_1(x)\ =\ \alpha_0\cdot 1+x\,,\ \ \) where \(\ \ \alpha_0\ \) is determined from the condition \(\ \ \langle\,p_0,p_1\rangle\ =\ 0\,.\)

   \(\ \langle\,p_0,p_1\rangle\ =\ \int_{-1}^{+1}\ 1\cdot(\alpha_0+x)\ dx\ \ =\ \ \int_{-1}^{+1}\ (\alpha_0+x)\,dx\ =\ 2\,\alpha_0\ =\ 0\,:\quad\alpha_0=0\,.\)

   \(\ p_1(x)\,=\,x\,;\qquad \|p_1\|^2\,=\,\int_{-1}^{+1}\ x^2\;dx\ =\ \left.\frac{1}{3}\ x^3\,\right|_{-1}^{+1}\ =\ \frac{2}{3}\ ;\qquad \|\,p_1\|\ =\ \sqrt{\frac{2}{3}}\ .\)

   \(\blacktriangleright\quad q_1(x)\ \,=\ \, \displaystyle\frac{1}{\|\,p_1\|}\ \ p_1(x)\ \,=\ \, \sqrt{\,\frac{3}{2}}\ \ x\ .\\\)

2. \(\ p_2(x)\ =\ \beta_0\cdot 1+\beta_1\cdot x+x^2\,,\ \ \) where \(\ \ \beta_0,\,\beta_1:\ \ \langle\,p_0,p_2\rangle\ =\ \langle\,p_1,p_2\rangle\ =\ 0\,.\)

   \(\ \langle\,p_0,p_2\rangle\ =\ \int_{-1}^{+1}\ (\beta_0+\beta_1\,x+x^2)\,dx\ =\ 2\,\beta_0+\frac{2}{3}\ =\ 0\,:\quad\,\beta_0\,=\ -\ \frac{1}{3}\,;\)

   \(\ \langle\,p_1,p_2\rangle\ =\ \int_{-1}^{+1}\ x\,(\beta_0+\beta_1\,x+x^2)\,dx\ =\ \frac{2}{3}\,\beta_1\ =\ 0\,:\qquad\beta_1\,=\,0\,.\)

   \(\ p_2(x)\,=\,x^2-\;\frac{1}{3}\ ;\)

   \(\ \|\,p_2\|^2\,=\,\int_{-1}^{+1}\ (x^2-\,\frac{1}{3})^2\,dx\ =\ \int_{-1}^{+1}\ (x^4-\frac{2}{3}\,x^2+\frac{1}{9})\,dx\ =\ \frac{8}{45}\ ;\quad \|\,p_2\|\ =\ \frac{2}{3}\ \sqrt{\frac{2}{5}}\,.\)

   \(\blacktriangleright\quad q_2(x)\ \,=\ \, \displaystyle\frac{1}{\|\,p_2\|}\ \ p_2(x)\ \ =\ \ \frac{3}{2}\ \ \sqrt{\,\frac{5}{2}}\ \ \left(x^2-\;\frac{1}{3}\right)\ \ =\ \ \frac{1}{2}\ \ \sqrt{\,\frac{5}{2}}\ \ (3\,x^2-\,1)\,.\\\)

3. \(\ p_3(x)\ =\ \widetilde{\gamma_0}\cdot 1+\gamma_1\cdot x+\gamma_2\cdot(x^2-\;\frac{1}{3})+x^3\ =\ \gamma_0+\,\gamma_1\,x\,+\,\gamma_2\,x^2\,+\,x^3\,,\)

   where \(\quad\gamma_0,\,\gamma_1,\,\gamma_2:\quad \langle\,p_0,p_3\rangle\,=\,\langle\,p_1,p_3\rangle\,=\,\langle\,p_2,p_3\rangle\,=\,0\,.\)

   \(\ \langle\,p_0,p_3\rangle\,=\,2\,\gamma_0+\frac{2}{3}\,\gamma_2\,=\,0\,;\quad \langle\,p_1,p_3\rangle\,=\,\frac{2}{3}\,\gamma_1+\frac{2}{5}\,=\,0\,;\quad \langle\,p_2,p_3\rangle\,=\,\frac{8}{45}\,\gamma_2\,=\,0\,;\)

   hence \(\quad\gamma_0=\gamma_2=0\,,\quad\gamma_1=-\ \frac{3}{5}\,.\)

   \(\ p_3(x)\,=\,x^3-\,\frac{3}{5}\,x\,;\qquad\|\,p_3\|\ =\ \frac{2}{5}\ \sqrt{\frac{2}{7}}\,.\)

   \(\blacktriangleright\quad q_3(x)\ \,=\ \, \displaystyle\frac{1}{\|\,p_3\|}\ \ p_3(x)\ \ =\ \ \frac{5}{2}\ \ \sqrt{\,\frac{7}{2}}\ \ \left(x^3-\,\frac{3}{5}\,x\right)\ =\ \frac{1}{2}\ \ \sqrt{\,\frac{7}{2}}\ \ (5\,x^3-3\,x)\,.\)

\[\mathcal{Q}\ \ =\ \ \left(\ \ \frac{1}{\sqrt{2}}\ ,\quad \sqrt{\,\frac{3}{2}}\ \ x\ ,\quad \frac{1}{2}\ \ \sqrt{\,\frac{5}{2}}\ \ (3\,x^2-\,1)\ ,\quad \frac{1}{2}\ \ \sqrt{\,\frac{7}{2}}\ \ (5\,x^3-3\,x)\ ,\ \ \dots\ \right)\]

Orthogonal Matrices and the QR Decomposition

Considerations of this section concern real domain, that is, real matrices \(\\\) and \(\,\) Euclidean spaces. Later we will translate the notions and theorems introduced below for complex domain.

Definition.

Matrix \(\ \boldsymbol{B}\,=\,[\,\boldsymbol{b}_1\,|\,\boldsymbol{b}_2\,|\,\dots\,|\, \boldsymbol{b}_n\,]\,=\,[\,\beta_{ij}\,]_{n\times n}\in M_n(R)\ \,\) is \(\,\) orthogonal, \(\,\) if

(2)\[\boldsymbol{B}^T\boldsymbol{B}\,=\,\boldsymbol{I}_n\,.\]

Properties of orthogonal matrices.

1. Taking determinant on both sides of the equation (2), we obtain

   \[\det\,(\boldsymbol{B}^T\boldsymbol{B})=\det\boldsymbol{B}^T\cdot\,\det\boldsymbol{B}= (\det\boldsymbol{B})^2\quad=\quad\det\boldsymbol{I}_n=1\,,\]

   and thus \(\,\det\boldsymbol{B}=\pm 1.\ \) Hence, the orthogonal matrix is nondegenerate, and thus invertible: \(\ \boldsymbol{B}^{-1}=\,\boldsymbol{B}^T\,.\ \) After multilplying this equality on the left by \(\ \boldsymbol{B}\ \) we obtain \(\ \boldsymbol{B}\boldsymbol{B}^T=\boldsymbol{I}_n\,.\ \) This means that orthogonal matrices satisfy equalities \(\ \ \boldsymbol{B}^T\boldsymbol{B}\,=\, \boldsymbol{B}\boldsymbol{B}^T=\boldsymbol{I}_n\,.\)

2. The condition \(\ \boldsymbol{B}\boldsymbol{B}^T=\boldsymbol{I}_n\ \) may be written as \(\ (\boldsymbol{B}^T)^T\boldsymbol{B}^T=\boldsymbol{I}_n\,,\ \) which implies that if \(\ \boldsymbol{B}\in M_n(R)\ \) is an orthogonal matrix, then orthogonal are also the transposed matrix \(\ \boldsymbol{B}^T\ \) and the inverse matrix \(\ \boldsymbol{B}^{-1}\,.\)

3. Let \(\ \boldsymbol{B}_1,\boldsymbol{B}_2\in M_n(R)\ \) be orthogonal matrices: \(\ \ \boldsymbol{B}_1^T\,\boldsymbol{B}_1=\boldsymbol{B}_2^T\,\boldsymbol{B}_2= \boldsymbol{I}_n\,.\ \) Then, by properties of the transpose of a matrix, we obtain

   \[(\boldsymbol{B}_1\boldsymbol{B}_2)^T(\boldsymbol{B}_1\boldsymbol{B}_2)\ =\ \boldsymbol{B}_2^T\,(\boldsymbol{B}_1^T\boldsymbol{B}_1)\,\boldsymbol{B}_2\ =\ \boldsymbol{B}_2^T\,\boldsymbol{I}_n\,\boldsymbol{B}_2\ =\ \boldsymbol{B}_2^T\,\boldsymbol{B}_2\ =\ \boldsymbol{I}_n\,.\]

   Hence, product of orthogonal matrices is also an orthogonal matrix. Since the identity matrix \(\ \boldsymbol{I}_n\ \) is orthogonal, we deduce

   Corollary 1.

   A set of orthogonal matrices of size \(\,n\ \) together with matrix multiplication comprises a (nonabelian) group.

4. The condition (2) may be rewritten in terms of matrix’ entries as

   \[\sum_{k\,=\,1}^n\ \beta_{ik}^T\;\beta_{kj}\,=\ \sum_{k\,=\,1}^n\ \beta_{ki}\;\beta_{kj}\,=\ \langle\,\boldsymbol{b}_i,\boldsymbol{b}_j\rangle\ =\ \delta_{ij}\,,\qquad i,j=1,2,\dots,n\,.\]

   This means that columns of the matrix \(\,\boldsymbol{B},\,\) interpreted as vectors of the space \(\,R^n\,,\ \,\) comprise an orthonormal system. Since the matrix \(\ \boldsymbol{B}^T\ \) is also orthogonal, this is also the case for rows of the matrix \(\ \boldsymbol{B}.\)

   Corollary 2.

   Matrix \(\ \boldsymbol{B}\in M_n(R)\ \) is orthogonal if \(\,\) and \(\,\) only if its columns \(\,\) (and also rows) \(\,\) comprise an orthonormal system in the space \(\,R^n\).

Consider a nondegenerate matrix \(\ \boldsymbol{A}\,=\, [\,\boldsymbol{a}_1\,|\,\boldsymbol{a}_2\,|\,\dots\,|\,\boldsymbol{a}_n\,]\in M_n(R).\ \) Its columns are linearly independent and comprise a basis \(\,\mathcal{B}=(\boldsymbol{a}_1\,,\,\boldsymbol{a}_2\,,\,\dots,\,\boldsymbol{a}_n)\,\) of the space \(\,R^n\). We will apply to this basis Gram-Schmidt orthogonalization in order to obtain an orthogonal basis \(\,\mathcal{P}=(\boldsymbol{y}_1\,,\,\boldsymbol{y}_2\,,\,\dots,\,\boldsymbol{y}_n)\ \ \) and \(\,\) an orthonormal basis \(\,\mathcal{Q}=(\boldsymbol{u}_1\,,\,\boldsymbol{u}_2\,,\,\dots,\,\boldsymbol{u}_n).\ \) Vectors of the bases \(\,\mathcal{P}\ \) and \(\ \mathcal{Q}\ \) are related by identities

\[\boldsymbol{u}_i\ =\ \,\frac{1}{\|\,\boldsymbol{y}_i\|}\ \ \boldsymbol{y}_i\,, \qquad i=1,2,\dots,n.\]

Our goal is to write the matrix \(\ \boldsymbol{A}\ \) as a product of an orthogonal matrix \(\\ \boldsymbol{Q}\,=\, [\,\boldsymbol{u}_1\,|\,\boldsymbol{u}_2\,|\,\dots\,|\,\boldsymbol{u}_n\,]\ \,\) and \(\,\) certain uppertriangular matrix \(\ \boldsymbol{R} : \ \boldsymbol{A}=\boldsymbol{Q}\boldsymbol{R}\,.\) Such a decomposition is called a QR decomposition.

The Gram-Schmidt process applied to the basis \(\,\mathcal{B}\ \,\) gives (\(j=2,\dots,n\)) :

\[\begin{split}\begin{array}{rcl} \boldsymbol{y}_1 & = & \boldsymbol{a}_1\,, \\ \boldsymbol{y}_j & = & \lambda_1\,\boldsymbol{y}_1\;+\ \lambda_2\,\boldsymbol{y}_2\;+\ \ldots\ +\ \lambda_{j-1}\,\boldsymbol{y}_{j-1}\;+\ \boldsymbol{a}_j\;\ = \\ & = & \boldsymbol{a}_j \;-\ \langle\,\boldsymbol{u}_1,\boldsymbol{a}_j\rangle\ \,\boldsymbol{u}_1 \;-\ \langle\,\boldsymbol{u}_2,\boldsymbol{a}_j\rangle\ \,\boldsymbol{u}_2 \;-\ \ldots \;-\ \langle\,\boldsymbol{u}_{j-1},\boldsymbol{a}_j\rangle\ \, \boldsymbol{u}_{j-1}\,; \end{array}\end{split}\]

(\(\,\lambda_i\ \) have been computed from the orthogonality conditions \(\,\langle\,\boldsymbol{y}_i,\boldsymbol{y}_j\rangle=0\,,\ \ i=1,2,\dots,j-1.\))

Hence, vectors from the initial basis \(\,\mathcal{B}\,\) may be represented by vectors from an orthonormal basis \(\,\mathcal{Q}\,\) in the following way:

\[\begin{split}\begin{array}{rcl} \boldsymbol{a}_j & = & \langle\,\boldsymbol{u}_1,\boldsymbol{a}_j\rangle\ \,\boldsymbol{u}_1 \;+\ \langle\,\boldsymbol{u}_2,\boldsymbol{a}_j\rangle\ \,\boldsymbol{u}_2 \;+\ \ldots \;+\ \langle\,\boldsymbol{u}_{j-1},\boldsymbol{a}_j\rangle\ \, \boldsymbol{u}_{j-1} \;+\ \,\boldsymbol{y}_j\ \,= \\ & = & \langle\,\boldsymbol{u}_1,\boldsymbol{a}_j\rangle\ \,\boldsymbol{u}_1 \;+\ \langle\,\boldsymbol{u}_2,\boldsymbol{a}_j\rangle\ \,\boldsymbol{u}_2 \;+\ \ldots \;+\ \langle\,\boldsymbol{u}_{j-1},\boldsymbol{a}_j\rangle\ \, \boldsymbol{u}_{j-1} \;+\ \,\|\,y_j\|\ \,\boldsymbol{u}_j\ \,= \\ & = & \langle\,\boldsymbol{u}_1,\boldsymbol{a}_j\rangle\ \,\boldsymbol{u}_1 \;+\ \langle\,\boldsymbol{u}_2,\boldsymbol{a}_j\rangle\ \,\boldsymbol{u}_2 \;+\ \ldots \;+\ \langle\,\boldsymbol{u}_{j-1},\boldsymbol{a}_j\rangle\ \, \boldsymbol{u}_{j-1} \;+\ \langle\,\boldsymbol{u}_j,\boldsymbol{a}_j\rangle\ \, \boldsymbol{u}_j\,. \end{array}\end{split}\]

These relations may be written in terms of an (uppertriangular) matrix \(\ \boldsymbol{R}\,=\,[\,\rho_{ij}\,]_{n\times n}\,:\)

\[ \begin{align}\begin{aligned}\begin{split}\rho_{ij}\ \ :\,=\ \ \left\{\ \begin{array}{ccc} \langle\,\boldsymbol{u}_i,\boldsymbol{a}_j\rangle & \text{dla} & i\leq j \\ 0 & \text{dla} & i > j \end{array} \right.,\quad i,j=1,2,\dots,n\,.\end{split}\\\boldsymbol{a}_j\ \;=\ \ \sum_{i\,=\,1}^j\ \langle\,\boldsymbol{u}_i,\boldsymbol{a}_j\rangle\ \boldsymbol{u}_i\ \ =\ \ \sum_{i\,=\,1}^n\ \rho_{ij}\;\boldsymbol{u}_i\,,\qquad j=1,2,\dots,n.\end{aligned}\end{align} \]

The last equality states that the \(\,j\)-th column of the matrix \(\,\boldsymbol{A}\ \) is a linear combination of columns of the matrix \(\,\boldsymbol{Q},\ \) with coefficients from the \(\,j\)-th column of the matrix \(\,\boldsymbol{R}.\ \) According to the column rule of matrix multiplication this is nothing else than the identity \(\ \boldsymbol{A}=\boldsymbol{Q}\boldsymbol{R}\ \) we were looking for.

Application of the QR decomposition. \(\ \)

Consider a Cramer system of equations over a field \(\,R\ \) with the coefficient matrix \(\,\boldsymbol{A}\in M_n(R)\ \) and a column vector of constants \(\,\boldsymbol{b}\in R^n\,.\ \) If \(\ \boldsymbol{A}=\boldsymbol{Q}\boldsymbol{R},\ \) then this system may be written as:

\begin{eqnarray*} \boldsymbol{A}\,\boldsymbol{x} & = & \boldsymbol{b}\,, \\ (\boldsymbol{Q}\boldsymbol{R})\,\boldsymbol{x} & = & \boldsymbol{b}\,, \\ \boldsymbol{Q}(\boldsymbol{R}\boldsymbol{x}) & = & \boldsymbol{b}\,. \end{eqnarray*}

Orthogonality of matrix \(\,\boldsymbol{Q}\,\) allows to replace an expensive operation of taking an inverse by a transpose: \(\ \ \boldsymbol{Q}^{-1}=\;\boldsymbol{Q}^T,\ \ \) and so

\[\boldsymbol{R}\,\boldsymbol{x}\ =\ \boldsymbol{Q}^T\,\boldsymbol{b}\,.\]

A system of equations with triangular matrix \(\,\boldsymbol{R}\ \) may be solved quickly by the “back-substitution” method.

For an illustration, we perform the QR decomposition for the matrix

\[\begin{split}\boldsymbol{A}\ =\ \left[\begin{array}{rrr} -2 & 8 & 19 \\ -2 & 11 & -14 \\ 1 & -7 & -8 \end{array}\right]\,.\end{split}\]

The Gram-Schmidt orthogonalization applied to columns of the matrix \(\,\boldsymbol{A}\ \) gives a matrix \(\,\boldsymbol{P}\ \) whose columns comprise an orthogonal system and further the target orthogonal matrix \(\,\boldsymbol{Q}.\ \\\) If one knows \(\,\boldsymbol{Q},\ \) the matrix \(\,\boldsymbol{R}\ \) may be easily computed from \(\ \boldsymbol{R}=\boldsymbol{Q}^T\boldsymbol{A}\,.\)

sage: A = matrix(QQ,[[-2,  8,  19],
                     [-2, 11, -14],
                     [ 1, -7,  -8]])

sage: P,Q = copy(A),copy(A)

sage: P[:,0] = A[:,0]
sage: Q[:,0] = P[:,0]/P[:,0].norm()

sage: P[:,1] = A[:,1] - Q.column(0)*A.column(1)*Q[:,0]
sage: Q[:,1] = P[:,1]/P[:,1].norm()

sage: P[:,2] = A[:,2] - Q.column(0)*A.column(2)*Q[:,0]\
                      - Q.column(1)*A.column(2)*Q[:,1]

sage: Q[:,2] = P[:,2]/P[:,2].norm()

sage: R = Q.T*A

sage: table([['$A$','','$Q$','','$R$'],[A,'=',Q,'*',R]])

\[\begin{split}\begin{array}{ccccc} A & & Q & & R \\ \\ \left(\begin{array}{rrr} -2 & 8 & 19 \\ -2 & 11 & -14 \\ 1 & -7 & -8 \end{array}\right) & = & \left(\begin{array}{rrr} -\textstyle\frac{2}{3} & -\textstyle\frac{2}{3} & \textstyle\frac{1}{3} \\ -\textstyle\frac{2}{3} & \textstyle\frac{1}{3} & -\textstyle\frac{2}{3} \\ \textstyle\frac{1}{3} & -\textstyle\frac{2}{3} & -\textstyle\frac{2}{3} \end{array}\right) & * & \left(\begin{array}{rrr} 3 & -15 & -6 \\ 0 & 3 & -12 \\ 0 & 0 & 21 \end{array}\right) \end{array}\end{split}\]

Sage contains a function (method) QR(), which performs the QR decomposition for a given matrix \(\,\boldsymbol{A}\ \) and returns a pair of matrices \(\,(\boldsymbol{Q},\boldsymbol{R})\,.\ \) It requires a ring which contains square roots and rational numbers (for example a field of algebraic numbers QQbar). Numerical computations should be done over a field RDF of real numbers with double precision.

sage: B = A.change_ring(RDF)
sage: (Q,R) = B.QR()
sage: show((Q.round(2),R))

\[\begin{split}\left(\left( \begin{array}{rrr} -0.67 & 0.67 & -0.33 \\ -0.67 & -0.33 & 0.67 \\ 0.33 & 0.67 & 0.67 \end{array}\right), \left(\begin{array}{rrr} 3.0 & -15.0 & -6.0 \\ 0.0 & -3.0 & 12.0 \\ -0.0 & -0.0 & -21.0 \end{array}\right)\right)\end{split}\]