First, consider the problem $\Sigma x = b$, where
$$
\Sigma = \pmatrix{\sigma_1\\& \ddots\\&&\sigma_r\\ &&&0\\&&&&\ddots\\&&&&&0}
$$
Note that $b$ is only in the range of $\Sigma$ if its entries $b_{r+1},\dots,b_n$ are all zero. Furthermore, you should be able to convince yourself (geometrically or otherwise) that the least squares solution must be
$$
x = (b_1/\sigma_1,\dots,b_r/\sigma_r,0,\dots,0)^T = \Sigma^+ b
$$
From there, note that
$$
U\Sigma V^T x = b \implies\\
\Sigma (V^T x ) = U^T b
$$
By the above argument, the least squares solution for $(V^T x)$ is given by
$V^T x = \Sigma^+ U^T b$. Noting that $\|V^T x\| = \|x\|$, we can use this to conclude that $x = (V \Sigma ^+ U^T)b$ must be the least squares solution (for $x$).
I hope you find this explanation sufficient.
For the full rank least squares problem, where $A \in \mathbb{K}^{m \times n},m>n=\mathrm{rank}(A)$ ($\mathbb{K}$ is the base field), the solution is $(A^T A)^{-1} A^T b$. This is a very bad way to approach the problem numerically for condition number reasons: you roughly square the condition number, so a relatively tractable problem with $\kappa=10^8$ becomes a hopelessly intractable problem with $\kappa=10^{16}$ (where we think about tractability in double precision floating point). The condition number also enters into convergence rates for certain iterative methods, so such methods often perform poorly for the normal equations.
The SVD pseudoinverse is exactly the same as the normal equations pseudoinverse i.e. $(A^T A)^{-1} A^T$. You simply compute it using the SVD and simplify. There is indeed a simplification; the end result is
$$(A^T A)^{-1} A^T = V (\Sigma^T \Sigma)^{-1} \Sigma^T V^T.$$
This means that if I know the matrix of right singular vectors $V$, then I can transform the problem of finding the pseudoinverse of $A$ to the (trivial) problem of finding the pseudoinverse of $\Sigma$.
The above is for the full rank problem. For the rank deficient problem with $m>n>\mathrm{rank}(A)$, the LS solution is not unique; in particular, $A^T A$ is not invertible. The usual choice is to choose the solution of minimal Euclidean norm (I don't really know exactly why people do this, but you do need some criterion). It turns out that the SVD pseudoinverse gives you this minimal norm solution. Note that the SVD pseudoinverse still makes sense here, although it does not take the form I wrote above since $\Sigma^T \Sigma$ is no longer invertible either. But you still obtain it in basically the same way (invert the nonzero singular values, leave the zeros alone).
One nice thing about considering the rank-deficient problem is that even in the full rank case, if $A$ has some singular value "gap", one can forget about the singular values below this gap and obtain a good approximate solution to the full rank least squares problem. The SVD is the ideal method for elucidating this.
The homogeneous problem is sort of unrelated to least squares, it is really an eigenvector problem which should be understood using different methods entirely.
Finally a fourth comment, not directly related to your three questions: in reasonably small problems, there isn't much reason to do the SVD. You still should not use the normal equations, but the QR decomposition will do the job just as well and it will terminate in an amount of time that you can know in advance.
Best Answer
If $A=\sum\limits_{i=1}^r\sigma_i u_iv_i^T$ is the SVD of $A$ (where $\sigma_i\neq 0$, $i=1,\ldots,r$), any vector $b$ can be expressed as a linear combination of $u_i$'s plus a term which is orthogonal to them, that is, $$ b=\sum_{i=1}^r\beta_iu_i+c, $$ where $\beta_i=u_i^Tb$ and $u_i^Tc=0$ for $i=1,\ldots,r$. Hence the residual $b-Ax$ is given by $$ b-Ax=c+\sum_{i=1}^r\left[\beta_iu_i-\sigma_i (v_i^Tx) u_i\right]=c+\sum_{i=1}^r\left[u_i^Tb-\sigma_i (v_i^Tx)\right]u_i. $$ Since the summands are mutually orthogonal, we can use the Pythagorean theorem to compute the norm of $b-Ax$ as $$ \|b-Ax\|_2^2=\|c\|_2^2+\sum_{i=1}^r\left[u_i^Tb-\sigma_i (v_i^Tx)\right]^2. $$ The first term is constant and, therefore, the residual norm is minimal if and only if $$\tag{1}v_i^Tx=\frac{u_i^Tb}{\sigma_i}, \quad i=1,\ldots,r.$$ Now any $x$ can be similarly to $b$ expressed as a linear combination of $v_i$'s plus a term orthogonal to these vectors: $$ x=\sum_{i=1}^r\xi_iv_i+y, $$ where $\xi_i=v_i^Tx$ and $v_i^Ty=0$ for $i=1,\ldots,r$. Hence, using (1), any $x$ of the form $$\tag{2} x=\sum_{i=1}^r\frac{u_i^Tb}{\sigma_i}v_i+y $$ minimizes the Euclidean norm of the residual $b-Ax$. There can be an infinite number of vectors $x$ minimizing the residual norm if the orthogonal complement of $\mathrm{span}(v_i)_{i=1}^r$ is nontrivial (that is, $A$ has a nontrivial kernel). Since $$ \|x\|_2^2=\sum_{i=1}^r\left(\frac{u_i^Tb}{\sigma_i}\right)^2+\|y\|_2^2, $$ we can see that setting $y=0$ in (2) leads to the solution, which, in addition to minimizing the residual, has the minimal norm.