title: Large-scale pseudoinverse
body: <p>I would like to compute the <a href="http://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse" rel="nofollow noreferrer">Moore–Penrose pseudoinverse</a> of an enormous matrix.  Ideally, I would like to do it on a matrix that has 23 million rows and 1000 columns, but if necessary I can reduce the number of rows to 4 million by only running on one part of my experiment.</p>
  
  <p>Obviously, loading the matrix in to memory and running SVD on it is not going to work.  <a href="http://en.wikipedia.org/wiki/Block_matrix_pseudoinverse" rel="nofollow noreferrer">Wikipedia</a> points to <a href="http://en.wikipedia.org/wiki/Krylov_subspace" rel="nofollow noreferrer">Krylov subspace</a> methods and mention the <a href="http://en.wikipedia.org/wiki/Arnoldi_iteration" rel="nofollow noreferrer">Arnoldi</a>, <a href="http://en.wikipedia.org/wiki/Lanczos_iteration" rel="nofollow noreferrer">Lanczos</a>, <a href="http://en.wikipedia.org/wiki/Conjugate_gradient" rel="nofollow noreferrer">Conjugate gradient</a>, <a href="http://en.wikipedia.org/wiki/GMRES" rel="nofollow noreferrer">GMRES</a> (generalized minimum residual), BiCGSTAB (biconjugate gradient stabilized), QMR (quasi minimal residual), TFQMR (transpose-free QMR), and MINRES (minimal residual) methods as being among the best Krylov subspace methods.  But I don't know where to go from here.  Is computing the pseudoinverse of such a huge matrix even feasible?  If so, using which algorithms or software libraries?  I have a large computing cluster available, so parallel approaches are welcome.</p>
  
  <p><a href="https://stackoverflow.com/questions/3169371/large-scale-regression-in-r-with-a-sparse-feature-matrix#3172815">This answer</a> points to the R package <a href="http://cran.r-project.org/web/packages/biglm/index.html" rel="nofollow noreferrer">biglm</a>.  Would that work?  Has anyone used it?  I normally work in Python, but don't mind using other languages and tools for this particular task.</p>