  We describe a parallel iterative least squares solver named \texttt{LSRN}
that is based on random normal projection. \texttt{LSRN} computes the
min-length solution to $\min_{x \in \mathbb{R}^n} \|A x - b\|_2$, where $A \in
\mathbb{R}^{m \times n}$ with $m \gg n$ or $m \ll n$, and where $A$ may be
rank-deficient. Tikhonov regularization may also be included. Since $A$ is only
involved in matrix-matrix and matrix-vector multiplications, it can be a dense
or sparse matrix or a linear operator, and \texttt{LSRN} automatically speeds
up when $A$ is sparse or a fast linear operator. The preconditioning phase
consists of a random normal projection, which is embarrassingly parallel, and a
singular value decomposition of size $\lceil \gamma \min(m,n) \rceil \times
\min(m,n)$, where $\gamma$ is moderately larger than 1, e.g., $\gamma = 2$. We
prove that the preconditioned system is well-conditioned, with a strong
concentration result on the extreme singular values, and hence that the number
of iterations is fully predictable when we apply LSQR or the Chebyshev
semi-iterative method. As we demonstrate, the Chebyshev method is particularly
efficient for solving large problems on clusters with high communication cost.
Numerical results demonstrate that on a shared-memory machine, \texttt{LSRN}
outperforms LAPACK's DGELSD on large dense problems, and MATLAB's backslash
(SuiteSparseQR) on sparse problems. Further experiments demonstrate that
\texttt{LSRN} scales well on an Amazon Elastic Compute Cloud cluster.
