Direct pseudoinverse on a different decomposition

 We hold the following truths to be self-evident (for real-matrices):

$$\displaystyle A = U \Sigma V^T$$

$$\displaystyle A^+=V \Sigma^{-1} U^T$$

$A^+$ is (up to quibbling) unique, and we see that it has a condition number $\kappa(A^+)=\kappa(\Sigma^{-1})$. When calculating a pseudoinverse from an SVD, we multiply all these components together, but all our instability is inherently baked into the central $\Sigma^{-1}$ term. Now, we might speculate, suppose we calculated $\Sigma^{-0.5}$, our condition number for this is the square root of that of $Sigma^{-1}$. With $U$ and $V$ being unitary we can* see that $V \Sigma^{-0.5}$ and $\Sigma^{-0.5} U^T$ have these same square root condition numbers and have only undergone a square root condition number multiplication operation. Multiplying these products together we then get a much more numerically stable calculation. A lot of handwaving went on there but feel free to verify in python that:

$$\displaystyle \textbf{EVAL: }\|A(V \Sigma^{-0.5})(\Sigma^{-0.5} U^T)A - A\| \leq \|A((V \Sigma^{-1}) U^T)A-A\|$$

To reiterate the above statement is a literal equality, but numerically there is a difference. Because pseudoinverses are often calculated using SVDs, it might make more sense to directly calculate $(V \Sigma^{-0.5})$ and $(\Sigma^{-0.5} U^T)$. Let us denote $L = (V \Sigma^{-0.5})$ and $R = (\Sigma^{-0.5} U^T)$.

We can calculate $L$ and $R$ directly by noting that $RAL=I$ and $L^TL = R R^T = \Sigma^{-1}$. We also are aided via $\|row(R,k)\| = \|col(L,k)\|$.

Modifying the standard householder based SVD algorithm, we can actually constrain it to avoid calculating the central bidiagonal matrix, instead creating an identity matrix via successive updates to $R$ and $L$ that diminish the deviation of $RAL-I$. Even if done in a fashion where $R$ and $L$ don't have their diagonalizable property or an equal distribution of singular values, as long as $RAL=I$ we find that $LR$ is still such that $LRx = A^+ x$.

There are some nice group properties to go along with this. When A is nonsingular and square we see that $RAL = LRA = ALR = I$, thus there are exact inverses for R, A, and L. Now the trick would be seeing if we can form a ring out of this. I have a sneaking suspicion that we can, and I have a copy of Hungerford that probably has the right lemma for it somewhere. I'll examine the algebraic properties in a later post.

*(I shouldn't need to prove this but if I must, just note that $V\Sigma^{-0.5} = V\Sigma^{-0.5}I$, which is of course a standard SVD format that we know carries the condition numbers in the central term.)