How can we compute Pseudoinverse for any Matrix

Question

If we have a system of linear equations $Ax=b$, then in case $A$ is invertible it easy to say that the solution is $x=A^{-1}b$. In all other cases, even if $A$ is invertible, the solution if it exist will be in the form $x=A^{+}b+(I-A^{+}A)w$, where $A^{+}$is called pseudoinvers of $A$, and $w$ is a vector of free parameters. My question here is how to compute $A^{+}$ for any matrix $A$? what is the way for doing that? It is easy to compute $A^{+}$ if the column are linearly independent (so that $m>=n$), $A^{+}=(A^{T}A)^{-1}A$, and also if the rows are linearly independent (so that $m<=n)$, $A^{+}=A^{T}(AA^{T})^{-1}$, but I dont know how to compute $A^{+}$ if $A$ is sequare non invertible matrix or if the columns or rows are not linearly independents. The above does not works. If anyone have any idea about computing pseudoinverse or if there is easier method for finding the solution $x$ for system of linear equations $Ax=b$, please help me in all cases, and please keep in mind that I want to find the general form for the solution if it exist not just for a particular $b$.

Answer 1

Let $A = U \Sigma V^T$ be an SVD of $A$, i.e., $U$ and $V$ are orthogonal and $\Sigma = \mathop{\rm diag}(\sigma_1, \sigma_2, \dots, \sigma_r, 0, \dots, 0)$ is real diagonal with $\sigma_k > 0$ for all $k=1,\dots,r$. Then

$$A^+ = V \Sigma^+ U^T,$$

with

$$\Sigma^+ = \mathop{\rm diag}(\sigma_1^{-1}, \sigma_2^{-1}, \dots, \sigma_r^{-1}, 0, \dots, 0).$$

Answer 2

The reason that the SVD is commonly used for the Pseudoinverse is because of the unitary operations. The SVD uses unitary operations only, which gives that the final results have the same error as any number in the environment, the machine error. This means that the small numbers may as well be treated as zero. This is generally useful. But it comes at the cost of being an iterative algorithm; there is no direct formula for the SVD.

However, if you have some reason to believe that the matrix is well conditioned (for its dimension whether it is square or not), or if you are using exact arithmetic, there is a direct (left inverse) formula: $$A^+ = A^\dagger =(A^\top A)^{-1}A^\top$$ I have written an algorithm here which avoids the extra matrix multiplications that the formula implies. I consider the algorithm as the Gaussian elimination equivalent for the non-square matrix. It operates on each column of the original matrix in turn, and a zero appears when the current column is already within the span of the previous columns. By comparison, the Gaussian elimination reaches this state when a column of all zeros appears. The difference in my algorithm though is that the column may just be ignored while noting that the dimensionality is reduced by one along with it.

For keeping track of what kind of errors could be in your problem, it is imperative to use something like the SVD. Thus if the results ever seem off when using my algorithm, it is likely related to such errors, as would happen when the matrix is ill conditioned.