Moore Penrose Pseudo-Inverse Fast Algorithm in R I want to apply Moore Penrose Pseudo-Inverse on my matrix, which is a 20,000 * 20,000 symmetric matrix with rank 19,999. I found ginv() function from the MASS package and MPInv() functioon from gnm package but they are not as fast as I expected.
Is there any other more efficient algorithmm to do Moore Penrose Pseudo-Inverse in R?
 A: The main problem you will run into here is not with the speed of the computation, but the ability of R to store matrices of sufficient size for what you want to do.  A numeric $N \times N$ matrix is stored as a vector with an additional set of dimensions.  The numeric vector takes up about $N^2 \times 8 \text{ bytes}$, so a numeric matrix of the size you want is beyond the storage limits in R.  This problem gets worse if you want to compute the singular value decomposition of the matrix, because you then have several matrices of this size.
In any case, setting this problem aside, the base package of R has a function svd to compute the singular value decomposition of a matrix.  It should be possible to use this function to compute the Moore-Penrose pseudo-inverse of a fairly large matrix.  Below I create a function to do this and test it on an $N \times N$ random matrix where I have set $N=5000$.  The computation time on my laptop computer is $$, although it may be faster or slower on other machines.
#Create function for computing Moore-Penrose pseudo-inverse via SVD
MPI <- function(matrix) {
  SVD <- svd(matrix);
  DDD <- rep(0, length(SVD$d));
  for (i in 1:length(DDD)) { DDD[i] <- ifelse(SVD$d[i] == 0, 0, 1/SVD$d[i]);  }
  SVD$v %*% diag(DDD) %*% t(SVD$u); }

#Compute the Moore-Penrose pseudo-inverse for an N x N random matrix
set.seed(75134903);
N   <- 5000;
MAT <- matrix(rnorm(N^2), nrow = N, ncol = N);
TT1 <- Sys.time()
INV <- MPI(MAT);
TT2 <- Sys.time()

#Check the result and computation time
identical(zapsmall(MAT %*% INV), diag(N));
[1] TRUE

#Check the computation time
TT2-TT1;

Note that the computation involves some rounding error, so you will have to ensure that it is sufficiently accurate for your purposes.  In the above simulation the computation is sufficiently accurate that the product of the matrix and the computed Moore-Penrose pseudo-inverse is close enough to the identity matrix that the zapsmall function makes them identical.
If you can find a way to compute and store the relevant matrices within the memory then it should be possible to obtain the Moore-Penrose pseudo-inverse of the matrix you are dealing with.  In my view, the main challenge will be dealing with memory issues rather than computation speed.
