QR decomposition (source code)

= QR decomposition
{c}
{title2=$A=QR$}
{wiki}

= QR factorization
{c}
{synonym}

A reduced <QR decomposition> of an $m\times k$ <matrix> over the <real numbers> with <full column rank>, $m\geq 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 <coefficients>
$$
r_{ij}=\begin{cases}
\langle v_j,q_i\rangle,&i<j,\\
\|u_j\|,&i=j,\\
0,&i>j.
\end{cases}
\qquad v_j=\sum_{i\leq j}r_{ij}q_i.
$$
This 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^TQ=I$, one also has $R=Q^TA$; $Q^T$ denotes the <transpose> and $I$ the <identity matrix>. Over the <complex numbers>, use the <conjugate transpose> $Q^*$ in these identities.

For the three-dimensional <Gram-Schmidt process> example, the exact <matrices> are
$$
A=\begin{pmatrix}1&1&0\\1&0&1\\0&1&1\end{pmatrix},\qquad
Q=\begin{pmatrix}
1/\sqrt2&1/\sqrt6&-1/\sqrt3\\
1/\sqrt2&-1/\sqrt6&1/\sqrt3\\
0&2/\sqrt6&1/\sqrt3
\end{pmatrix},
$$
$$
R=\begin{pmatrix}
\sqrt2&1/\sqrt2&1/\sqrt2\\
0&\sqrt{3/2}&1/\sqrt6\\
0&0&2/\sqrt3
\end{pmatrix}.
$$
For example, the last column reconstructs as
$$
\frac{q_1}{\sqrt2}+\frac{q_2}{\sqrt6}+\frac{2q_3}{\sqrt3}
=\frac12(1,1,0)^T+\frac16(1,-1,2)^T+\frac23(-1,1,1)^T
=(0,1,1)^T.
$$
The 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.

\Image[/gram-schmidt-qr-factorization.png]
{title=The computed <QR factorization>: input <matrix> $A$, <orthonormal> columns of $Q$, and the <upper triangular matrix> $R$ of <orthogonal projection> <coefficients>}

For a <linear least-squares problem>, the <orthogonal projection> of $b$ onto the <column space> of $A$ is $QQ^Tb$. Hence the minimizing <coefficients> solve the <linear system> $Rx=Q^Tb$ 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\times2$ block $R_2$ of $R$. For $b=(0,0,1)^T$, the <linear system> $R_2x=Q_2^Tb$ is
$$
\begin{pmatrix}\sqrt2&1/\sqrt2\\0&\sqrt{3/2}\end{pmatrix}
\begin{pmatrix}x_1\\x_2\end{pmatrix}
=\begin{pmatrix}0\\2/\sqrt6\end{pmatrix},
\qquad x_2=\frac23,\quad x_1=-\frac13.
$$
Thus $A_2x=(1/3,-1/3,2/3)^T$, agreeing with the <orthogonal projection> computed above, and the minimum <Euclidean norm> of $A_2x-b$ is $1/\sqrt3$. The <QR decomposition> separates building an <orthonormal basis> from recovering <coefficients> in the original <basis>.https://ocw.mit.edu/courses/18-065-matrix-methods-in-data-analysis-signal-processing-and-machine-learning-spring-2018/resources/lecture-9-four-ways-to-solve-least-squares-problems/[Strang, MIT 18.065, Lecture 9: Four Ways to Solve Least Squares Problems]{ref}