How to efficiently compute an "inverse" of a tall matrix

Question

The following question came to mind when looking at this question and thinking about the approach I took here.

Suppose that we are given an $m \times k$ matrix $M$ with $k<m$ (potentially $k$ is much smaller than $m$) and linearly independent columns. How can one (efficiently) find a square, invertible matrix $P$ such that $$ PX = \pmatrix{I\\ 0}? \tag{*} $$ What if we also require knowledge of the full form of $P^{-1}$? That is, what if we also need to know the matrix $Y$ in $$ P^{-1} = \pmatrix{X & Y}? \tag{**} $$

I was primarily thinking about this problem for matrices with real-number or complex-number entries, but approaches for other fields are welcome.

One approach to compute such a $P$ (which is the first that I had thought of but I suspect is significantly less efficient than is optimal) is to extend the columns of $X$ into a complete basis, which is to say find an $m \times (m-k)$ matrix $Y$ such that $[X\ \ Y]$ is invertible. From there, one could simply take $P = [X \ \ Y]^{-1}$.

Notably, finding such a matrix $Y$ is itself a noteworthy task; one approach is to row-reduce $[X \mid I_m]$ and extract pivot columns. Another "probabilistic" approach is to just select $Y$ randomly since if $Y$ is selected using an absolutely continuous probability measure, $[X\ \ Y]$ will be invertible with probability $1$.

As it turns out (and I'm only noticing this now), row reducing $[X \mid I_m]$ will produce a matrix of the form $[R \mid P]$ where $R$ denotes the right-hand side of $(*)$ and $P$ satisfies $(*)$, i.e. $PX = R$. That already brings the complexity down from $O(m^3)$ to $O(k^3)$. It's not clear to me if $P^{-1}$ could quickly be constructed here.

I have some other approaches that come to mind that I'll be sharing in my own answer, but I am interested in seeing any ideas and, if anyone recognizes it, knowing whether this problem has been studied.

Answer 1

Another approach: begin by solving $P_1 X = I_k$, i.e. finding a "left inverse" to $X$. For that, one can take the Moore-Penrose pseudoinverse $P_1 = X^+$.

From there, we can take $$ P = \pmatrix{P_1\\P_2}, $$ where the rows of $P_2$ must be orthogonal to the columns of $X$. That is, one must (independently) solve the problem $P_2 X = 0_{(m-k) \times k}$. For this problem, one solution is to select a spanning set of rows from $I - XX^+$. Another approach is to use a full $QR$-decomposition $X = QR$ (where $Q$ is $m \times m$) and take the last $m-k$ columns of $Q$ to be the rows of $P_2$.

One could also rely entirely on the $QR$ decomposition. First, compute a decomposition $$ X = \pmatrix{Q_1 & Q_2} \pmatrix{R_1\\0} = Q_1R_1. $$ From there, we can solve $P_1 (Q_1R_1) = I$ by taking $P_1 = R_1^{-1}Q_1^\top$, where the upper-triangular structure of $k$ makes its inversion $O(k^2)$ instead of $O(k^3)$. As in the previously discussed approach, we can take $P_2 = Q_2^\top$. We can also recover the corresponding $Y$ as follows: $$ P = \pmatrix{R_1^{-1}&0\\0&I} Q^\top \implies\\ P^{-1} = Q \pmatrix{R_1&0\\0&I} = \pmatrix{X & Q_2}. $$