For least squares problems of minimizing
$\|\boldsymbol{b}-A\boldsymbol{x}\|_2$ where $A$ is a large sparse
$m\times n$ ($m\ge n$) matrix, the common method is to apply the
conjugate gradient method to the normal equation
$A^{\mathrm{T}} A\boldsymbol{x}=A^{\mathrm{T}} \boldsymbol{b}$.
However, the condition number of $A^{\mathrm{T}} A$ is square of that of $A$,
and convergence becomes problematic for severely ill-conditioned
problems even with preconditioning.
In this paper, we propose two methods for applying the GMRES method to
the least squares problem by using an $n\times m$ matrix $B$.
We give the necessary and sufficient condition that $B$ should satisfy
in order that the proposed methods give a least squares solution.
Then, for implementations for $B$, we propose an incomplete QR
decomposition IMGS($l$).
Numerical experiments showed that the simplest case $l=0$ gives the
best results, and converges faster than previous methods for severely
ill-conditioned problems. The preconditioner IMGS(0) is equivalent to
the case $B=(\diag (A^{\mathrm{T}} A))^{-1} A^{\mathrm{T}}$,
so $(\diag (A^{\mathrm{T}} A))^{-1} A^{\mathrm{T}}$ was
the best preconditioner among
IMGS($l$) and Jennings' IMGS($\tau$).
On the other hand, CG-IMGS(0) was the fastest
for well-conditioned problems.