  We develop the first stochastic incremental method for calculating the
Moore-Penrose pseudoinverse of a real matrix. By leveraging three alternative
characterizations of pseudoinverse matrices, we design three methods for
calculating the pseudoinverse: two general purpose methods and one specialized
to symmetric matrices. The two general purpose methods are proven to converge
linearly to the pseudoinverse of any given matrix. For calculating the
pseudoinverse of full rank matrices we present two additional specialized
methods which enjoy a faster convergence rate than the general purpose methods.
We also indicate how to develop randomized methods for calculating approximate
range space projections, a much needed tool in inexact Newton type methods or
quadratic solvers when linear constraints are present. Finally, we present
numerical experiments of our general purpose methods for calculating
pseudoinverses and show that our methods greatly outperform the Newton-Schulz
method on large dimensional matrices.
