If $A^TA$ is not positive definite (or it's not invertible, or equivalently some of its eigen values are zero), then we can either use the Moore–Penrose inverse, or we can use Tikhonov regularization.
Tikhonov regularization calculates the regular inverse of $(A^TA+\epsilon I)$ for some small $\epsilon$, and uses this inverse instead.
Moore–Penrose inverse is basically based on the SVD decomposition of $A^TA$ say $V\Lambda V^T$, where $V$ is a unitary matrix and $\Lambda$ is a diagonal matrix. Then $(A^TA)^{-1}=V\Lambda^{-1} V^T$, if an element $\lambda_i$ (a singular value) on the diagonal of $\Lambda$ is not zero then we replace it with $1/\lambda_i$ in $\Lambda^{-1}$, but if it's zero, we set the corresponding element to zero in $\Lambda^{-1}$.
Now we can answer your question, about why people use iterative solutions. There can be several reasons, one is when you have other constraints on the variable $x$, then $(A^TA)^{-1}A^Tb$ is not the solution. But the main reason is that for large scale $A$, calculating the SVD decomposition for example can be expensive, and at the same time it maybe not memory efficient. Some times you don't have access to full $A$, and can only observe a sparse representation of it, and so on.
Computing the singular value decomposition (SVD) of symmetric, rank-$1$ matrix $\rm A$,
$$\mathrm A = \begin{bmatrix} 1 & -1\\ -1 & 1\end{bmatrix} = \begin{bmatrix} 1\\ -1\end{bmatrix} \begin{bmatrix} 1\\ -1\end{bmatrix}^\top = \left( \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1\\ -1 & 1\end{bmatrix} \right) \begin{bmatrix} 2 & 0\\ 0 & 0\end{bmatrix} \left( \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1\\ -1 & 1\end{bmatrix} \right)^\top$$
Hence, the pseudoinverse of $\rm A$ is
$$\mathrm A^+ = \left( \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1\\ -1 & 1\end{bmatrix} \right) \begin{bmatrix} \frac12 & 0\\ 0 & 0\end{bmatrix} \left( \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1\\ -1 & 1\end{bmatrix} \right)^\top = \color{blue}{\frac14 \mathrm A}$$
Best Answer
The Moore-Penrose inverse of $A\in \mathbb{C}^{m\times n}$, denoted by $A^{+}$, is the unique matrix $X$ satisfying the following four Penrose equations. \begin{eqnarray} (i)~ AXA = A,~ (ii)~ XAX = X, ~ (iii)~ (AX)^* = AX,~ (iv)~ (XA)^* = XA \end{eqnarray}
There are various methods to find out Moore-Penrose inverse of a matrix. For example Rank factorization method, singular value decomposition method, QR decomposition etc. If $Ax = b$ is inconsistent then $x = A^{+}b$ represents the least square solution of minimum 2 norm.
In general $x = A^{+}b$ represents the least square solution to system $Ax = b$. For more detailed information refers to generalized inverse by Ben-Israel http://www.amazon.com/Generalized-inverses-applications-Adi-Ben-Israel/dp/0882759914