Construction of QR decomposition for a singular matrix

Question

Consider a matrix $A \in \mathbb{C}^{n \times n}$. Suppose it's non-singular. Then it's columns $a_1, a_2, ..., a_n$ form a basis in $\mathbb{C}^n$.

Let's apply Gram-Schmidt process to them. We will obtain some basis $q_1, q_2, ..., q_n$. The transition matrix $S$ from the basis $a$ to the basis $q$ is upper-triangular(by construction). Now, consider the matrix $Q=\big[ q_1|q_2|...|q_n \big]$. It's columns are orthonormal $\implies$ $Q$ is unitary.

And now we know that the matrix $A$ can be represented as: $$ AS=Q \\ A=QS^{-1} $$

$S$ is invertible since it is a transition matrix. $S$ is upper-triangular $\implies$ $S^{-1}$ is upper-triangular as well.

E.g. we have shown that each non-degenerate matrix $A$ can be represented as a product of unitary and upper-triangular matrices.

But what if $A$ is singular? I've tried the following:

Let's select a basis of column space of $A$: $a_1, ..., a_r$ (they are columns of $A$). And also let's reorder columns of $A$ s.t. these linearly independent vectors $a_1, ..., a_r$ are the first $r$ vectors of the matrix. So, we've done permutation of columns of $A$. This operation can be represented by right-multiplication of A with some matrix $P$: $A \rightarrow AP$.

Now we apply Gram-Schmidt to the basis of column space(vectors $a_1,...,a_r$), and after extend the resulting basis to an ONB of $\mathbb{C}^n$:

$$q_1,q_2, ..., q_n$$

We have constructed matrix $Q=\big[ q_1 | q_2 | ... | q_n \big]$ that is unitary(since it's columns form an ONB). The matrix $S$ s.t. $AP=QS$ is upper-triangular. Moreover, it looks like: $$ \begin{bmatrix} s_{11} & s_{12} & s_{13} & ... & s_{1r} & ... & s_{1n} \\ 0 & s_{22} & s_{23} & ... & s_{2r} & ... & s_{2n} \\ 0 & 0 & s_{33} & ... & s_{3r} & ... & s_{3n} \\ ... & ... & ... & ... & .... & ... & ... \\ 0 & 0 & 0 & ... & s_{rr} & ... & s_{rn} \\ 0 & 0 & 0 & ... & 0 & ... & 0\\ 0 & 0 & 0 & ... & 0 & ... & 0\\ ... & ... & ... & ... & .... & ... & ... \\ 0 & 0 & 0 & ... & 0 & ... & 0\\ \end{bmatrix} $$

But how to show now that the matrix $SP^{-1}$ is upper-triangular? For me, it seems like it's not. Can you see a mistake in my thoughts? Maybe, you have other suggestions for construction of $QR$-decomposition in case of singular matrix? I'm stuck with this issue for a week now and I have no ideas what to do. Thanks!

Answer 1

conceptually the simplest approach using Gram Schmidt, is to lean on what you already know. Suppose $\text{rank}\big(A\big)=r$.

1.) run Gram Schmidt on $A$ as you normally do, and each time you run into a linearly dependent column, discard it. The end result is

$$A' = Q'R'=\bigg[\begin{array}{c|c|c|c} \mathbf q_1' & \mathbf q_2' &\cdots & \mathbf q_{r}'\end{array}\bigg]R'$$

where $A'$ is $n\times r$ and $R$ is $r\times r$ and upper triangular. This is equivalent to running QR facotrization on a tall skinny matrix with full column rank.

2.) Now extend $Q'$ to an orthonormal basis (i.e. 'find' $n-r$ linearly independent vectors not in the span of $A'$ and run more Gram Schmidt. Common places to look -- left nullspace of $A$ or use a random number generator, ref here: In $A = QR$ factorization do we need basis vectors from the left null space?)

3.) Now one at a time re-insert the deleted columns to recreate $A$. E.g. suppose $\mathbf a_1$ was non-zero but $\mathbf a_2$ was linearly dependent (and hence deleted when constructing $A'$). So re-insert $\mathbf a_2$ on the LHS

$$A'=\bigg[\begin{array}{c|c|c|c} \mathbf a_1' & \mathbf a_2' &\cdots & \mathbf a_{r}'\end{array}\bigg]=\bigg[\begin{array}{c|c|c|c} \mathbf a_1 & \mathbf a_2' &\cdots & \mathbf a_{r}'\end{array}\bigg]\rightarrow \bigg[\begin{array}{c|c|c|c|c} \mathbf a_1 & \mathbf a_2& \mathbf a_2' &\cdots & \mathbf a_{r}'\end{array}\bigg]$$

Now to make dimensions match on the RHS, first insert some $\mathbf q_j$ that you found when doing step (2) in the appropriate column (i.e. 2nd column)
$$\bigg[\begin{array}{c|c|c|c} \mathbf q_1' & \mathbf q_2' &\cdots & \mathbf q_{r}'\end{array}\bigg]\rightarrow \bigg[\begin{array}{c|c|c|c} \mathbf q_1'& \mathbf q_j & \mathbf q_2' &\cdots & \mathbf q_{r}'\end{array}\bigg]$$

and you need to insert a row and column into $R'$ to make dimensions match (here $k=1$)
$$R'=\begin{bmatrix} R_{k\times k}' & * \\ \mathbf 0 & R_{r-k\times r-k}'\end{bmatrix}\rightarrow \begin{bmatrix} R_{k+1\times k+1}' & * \\ \mathbf 0 & R_{r-k\times r-k}'\end{bmatrix}$$ and the RHS is still upper triangular because its two diagonal blocks are upper triangular, i.e. one of them is unchanged (bottom right) and the other (top right) had a row and column inserted on the bottom and right respectively, but this is still upper triangular because $\mathbf a_2$ is written a linear combination of $\mathbf a_i$ with lower indices (i.e. it is linearly dependent).

Thus the structure holds one step at a time and after finitely many steps you recover $A=QR$

remark 1:
QR factorization is not unique for non-invertible $A$

remark 2:
a much cleaner way of getting to this result involves using Householder matrices. The proof is essentially is the same as when using Householder matrices for QR factorization of invertible A, except you skip the step at iteration $j$ if the block submatrix you are looking at has a zero in column j (which occurs iff $\mathbf a_j = \sum_{k=1}^{j-1}x_k \mathbf a_k$).

Answer 2

This will be a similar approach as @user8675309, but slightly different.

In your setup, $SP^{-1}$ is likely not upper triangular. There is a simple fix. Just do not reorder the columns. Leave as they are.

Consider a reduced row echelon form (rref) of the matrix $A$. It encodes the linear relations between columns of $A$. The columns of $A$ including the pivot positions, say $a_1, \ldots, a_r$ from the left to right, are the basis of columns space of $A$.

Construction of $Q$

To construct $Q$, first extend $a_1, \ldots, a_r$ to the basis of $\mathbb{C}^n$. Then apply Gram-Schmidt to obtain an orthonormal basis $q_1, \ldots, q_n$ of $\mathbb{C}^n$. Take $Q= [q_1 | \cdots | \ q_n]$.

The matrix $R$

To obtain $R$, consider columns of rref of $A$ between two pivot positions $j$ th and $j+1$ th (necessarily $j<r$). We write from the $(i-1)$th column to the $(i+k+1)$th column of the rref of $A$. $$ \begin{bmatrix} 0 & a_{1,i} & \cdots & a_{1,i+k} & 0\\ \vdots & \vdots & \cdots & \vdots & \vdots\\ 0 & a_{j-1,i} & \cdots & a_{j-1, i+k} & 0\\ 1 & a_{j,i} & \cdots & a_{j,i+k} & 0\\ 0 & 0 & \cdots & 0 & 1 \end{bmatrix} $$

From the rref of $A$, we see that from the $i$ th column to the $i+k$ th column of $A$ are linear combinations of $a_1, \ldots a_j$. The $i-1$ th column of $A$ is $a_j$, and $i+k+1$ th column of $A$ is $a_{j+1}$. Since $\mathrm{span} (a_1,\ldots, a_j) = \mathrm{span} (q_1, \ldots, q_j)$, the $i-1$ to $i+k$ th columns of $A$ are linear combinations of $q_1, \ldots, q_j$, and $i+k+1$ th column is a linear combination of $q_1, \ldots, q_{j+1}$. Let $A^{(j)}$ be $j$-th column of $A$. Then we have $$ \begin{align} A^{(i-1)}&=r_{1,i-1} q_1+ \cdots + r_{j,i-1} q_j\\ A^{(i)}&=r_{1,i} q_1+\cdots + r_{j,i} q_j \\ &\vdots \ \ \vdots \ \\ A^{(i+k)}&=r_{1,i+k} q_1 + \cdots r_{j,i+k} q_j\\ A^{(i+k+1)}&=r_{1,i+k+1} q_1+\cdots + r_{j,i+k+1} q_j + r_{j+1, i+k+1} q_{j+1} \end{align} $$ for some $r_{j,j'}\in\mathbb{C}$. That is, $$ \begin{bmatrix} | & | & \cdots & |\\ A^{(i-1)} & A^{(i)} & \cdots & A^{(i+k+1)} \\ | & | & \cdots & | \end{bmatrix} $$ $$ =\begin{bmatrix} | & | & \cdots & |\\ q_{1} & q_2 & \cdots & q_{j+1}\\ | & | & \cdots & | \end{bmatrix} \begin{bmatrix} r_{1,i-1} & r_{1,i} & \cdots & r_{1,i+k} & r_{1,i+k+1}\\ \vdots & \vdots & \cdots & \vdots & \vdots\\ r_{j-1,i-1} & r_{j-1,i} & \cdots & r_{j-1, i+k} & r_{j-1,i+k+1}\\ r_{j,i-1} & r_{j,i} & \cdots & r_{j,i+k} & r_{j,i+k+1}\\ 0 & 0 & \cdots & 0 & r_{j+1,i+k+1} \end{bmatrix}. $$

Thus, $R$ must be upper triangular. Moreover, $R$ has the same pivot positions as the rref of $A$. We may have nonzero entries above pivots of $R$ though.

The computer algebra system will obtain $R$ by computing $\overline{Q}^T A$.

This approach is implemented in the online calculator.

I used the matrix $$ A =\begin{bmatrix} 1 & 3 & 2 & 10 & 2\\ 2 & 4 & 2 & 14 & 5 \\ -1 & -2 & -1 & -7 & 3 \\ 0 & 1 & 1 & 3 & 2\\ 4 & 1 & -3 & 7 & 5 \end{bmatrix} $$ The matrix has

$A^{(1)}$, $A^{(2)}$, $A^{(5)}$ as basis of column space, and

$A^{(3)}=A^{(2)}-A^{(1)}$

$A^{(4)}=A^{(1)}+3A^{(2)}$.