3

What is the step by step numerical approach to calculate the pseudo-inverse of a matrix with M rows and N columns, using LU decomposition?

So far, I have found this, but it uses singular value decomposition.

Community
  • 1
sgarizvi
  • 133
  • I see you have linked to the question for which I have given an answer that exactly answers this question now. My answer there uses LU to calculate $A^\dagger=(A^\top A)^{-1}A^\top$ without explicitly calculating the cross-product matrix. It does LU, but calculates the columns as needed, rather than calculating the matrix, then updating the matrix, then reading the needed values from that matrix. It avoids all that extra and uneccessary work. – adam W Mar 01 '13 at 17:47

1 Answers1

4

Given an $m\times n$ matrix $\mathbf A$, there are a number of non-SVD methods for computing the Moore-Penrose inverse. Most of them require an accurate determination of "numerical rank"; to drive the point home, would you say the matrix

$$\begin{pmatrix}1&&\\&1&\\&&\varepsilon^{2/3}\end{pmatrix}$$

(where $\varepsilon$ is machine epsilon) has rank $3$ or rank $2$? It is well-known that no method based on Gaussian elimination is foolproof with respect to rank determination, and thus SVD methods are preferable.

Of course, if your matrix has full rank, the classical formula applies:

$$\mathbf A^\dagger=(\mathbf A^\top\mathbf A)^{-1}\mathbf A^\top$$

but of course the formation of the cross-product matrix is wrought with danger on its own.

Having cautioned you on why persisting on the use of Gaussian elimination is unsound, let me mention a few papers (which you could have found on your own by searching with, say, Google Scholar). As a start you will want to look at this gentle introduction. From here, you will want to look at paper 1, paper 2, paper 3, and paper 4, among others. (You can find more using the search terms "Moore Penrose inverse" or "generalized inverse" along with elimination.)

As a bonus, you might also be interested in the Schulz iteration (see e.g. paper 5),

$$\mathbf Z_{k+1}=\mathbf Z_k(2\mathbf I-\mathbf A\mathbf Z_k)$$

which is nothing more than the Newton-Raphson method applied to matrix inversion.

Delete Me Please
  • 56