Algorithm for Constructing a Projection Matrix onto the Null Space

linear algebramatricesprojectionprojection-matricesvector-spaces

I am using MATLAB to write functions that generate projection matrices onto the column space, col($A$), and null space, null($A$), of an input matrix $A$.

I have the column space projection working but when I apply a slightly altered algorithm to write the function generating the null space projection matrix, my null space projector results do not fulfill the general conditions for a projection matrix, as I understand them.

Here is how I algorithmically proceed for (successfully) generating projection matrices onto col(A):

  1. Take input $A$ which is an ($m\times n)$ matrix.

  2. Run B = span_basis(A). This function creates a matrix B whose columns are basis vectors spanning col(A).

  3. Run [Q,~] = qr(B). This grabs matrix Q from B's QR-factorization.

  4. Re-define Q to Q = Q(all rows, columns 1 thru rank(A)). This makes sure that Q's columns are now only the vectors constituting the ortho-normal basis of col(A).

  5. I then use nested for-loops to generate the entries of the projection matrix $P$ onto col(A). The loops run on indices $i$ and $j$ from 1 to size(Q). Each entry $P(i,j)$ is the dot product between rows of $Q$:

    $P(i,j) = Q(i,:) \cdot Q(j,:)$

Once I've done these things, I have a projection matrix $P$ onto the column space of $A$. I then test

(i) $P = P^T$

(ii) $P^2=P$

(iii) rank(A) = rank(P)

(iv) $PA = A$

Those all work out. Thus, $P$ is a projection matrix.

However, problems arise when I repeat these same steps – slightly altered – to get the projection onto null(A):

  1. Take input $A$ which is an ($m\times n)$ matrix. Make $A^T$ which is $(n \times m)$.

  2. Run B = span_basis($A^T$). This creates a matrix B whose columns are basis vectors spanning col($A^T$).

  3. Run [Q,~] = qr(B). This grabs matrix Q from B's QR-factorization.

  4. Re-define Q to Q = Q(all rows, columns 1 thru rank($A^T$)). This makes sure that Q's columns are the vectors constituting the ortho-normal basis of col($A^T$).

  5. I then use nested for-loops to generate the entries of the projection matrix $P$ onto col($A^T$). The loops run on indices $i$ and $j$ from 1 to size(Q). Each $P(i,j)$ is the dot product between rows of $Q$:

    $P(i,j) = Q(i,:) \cdot Q(j,:)$

  6. I then finally get the projection onto null(A) as $\tilde{P} = \mathbb{I} -P $.

I then check rank($A^T$) = rank($\tilde{P}$), $\tilde{P}^2=\tilde{P}$ and $\tilde{P}^T = \tilde{P}$.

However, it is my understanding that, in addition to these conditions, I should be getting $\tilde{P}A^T=A^T$. However, this is not the case. Instead, when I grab the absolute value of the maximum difference between entries of $\tilde{P}A^T$ and $A^T$, I always get something around 1 (usually 0.99…), i.e.,

max(abs($\tilde{P}A^T-A^T$), [], "all") = 0.99…

That means that $\tilde{P}A^T$ and $A^T$ are not equal, as – at least to my understanding – they should be.

Where am I going wrong?

Thank you.

Best Answer

Your algorithm is fine. Steps 1-4 is equivalent to running Gram-Schmidt on the columns of $A$, weeding out the linearly dependent vectors. The resulting matrix $Q$ has columns that form an orthonormal basis whose span is the same as $A$. Thus, projecting onto $\operatorname{colspace} Q$ is equivalent to projecting onto $\operatorname{colspace} A$. Step 5 simply computes $QQ^\top$, which is the projection matrix $Q(Q^\top Q)^{-1}Q^\top$, since the columns of $Q$ are orthonormal, and hence $Q^\top Q = I$.

When you modify your algorithm, you are simply performing the same steps on $A^\top$. The resulting matrix $P$ will be the projector onto $\operatorname{col} (A^\top) = (\operatorname{null} A)^\perp$. To get the projector onto the orthogonal complement $\operatorname{null} A$, you take $\tilde{P} = I - P$.

As such, $\tilde{P}^2 = \tilde{P} = \tilde{P}^\top$, as with all orthogonal projections. I'm not sure how you got $\operatorname{rank} \tilde{P} = \operatorname{rank} A$; you should be getting $$\operatorname{rank} \tilde{P} = \dim \operatorname{null} A = n - \operatorname{rank} A.$$ Perhaps you computed $\operatorname{rank} P$ instead?

Correspondingly, we would also expect $P$, the projector onto $\operatorname{col}(A^\top)$, to satisfy $PA^\top = A^\top$, but not for $\tilde{P}$. In fact, we would expect $\tilde{P}A^\top = 0$; all the columns of $A^\top$ are orthogonal to $\operatorname{null} A$, so projecting each of them should produce $0$. Or, more directly, $$PA^\top = A^\top \implies IA^\top - PA^\top = 0 \implies (I - P)A^\top = 0 \implies \tilde{P}A^\top = 0.$$ I would double check these things: check that $\tilde{P} A^\top$ is indeed $0$, and that $\operatorname{rank} \tilde{P} = n - \operatorname{rank} A$.

Related Question