Finding a solution of a non-square linear system with SVD According to this solution, it seems that it's possible to find one solution of a non-square linear system that has infinitely many solutions with the following method.

*

*Write the linear system as $A X = B$.
(A is non invertible, since A is non-square)


*Consider the singular value decomposition of $A$:
$$A = U\  D\  ^t V$$
where $D$ is a $(m, n)$ diagonal matrix with non-negative real numbers on the diagonal.


*Let $\Delta$ be a $(m, n)$ diagonal matrix, with:
$\Delta_{i,i} = 0$ if $D_{i,i} = 0$
$\Delta_{i,i} = \frac{1}{D_{i,i}}$ if $D_{i,i} \neq 0$


*One solution of the initial linear system is given by:
$$ X= V\ \Delta\ ^tU\ B$$
Is this true in general? (it seems to be true at least for this example). How can you prove this?

Note: It seems that  $A^+ = V\ \Delta\ ^tU$ is the pseudoinverse of the matrix $A$, see this article.
PS: Here are some ideas but I'm looking for a real proof.

Note2: It seems that this problem is in fact made of 2 steps:

*

*Prove that if $A = U \Sigma ^tV$ is the SVD of $A$, then $A^+ = V \Sigma^+\, ^tU$ respects the conditions defining the pseudo inverse, i.e. $A A^+ A = A$, $A^+ A A^+ = A^+$, and the fact that $A^+ A$ and $A A^+$ are symmetrical matrices.


*Prove that if $(1) A X = B$ has infinitely many solutions, then $X = A^+ B$ is a solution of (1), i.e. we have to prove that $A A^+ B = B$.
 A: If we have the SVD $A=UDV^*$, with $U$, $V$ unitary, $D$ diagonal, then the matrix $V\Delta U^*$ (with $\Delta$ as in your post) is equal to the Moore-Penrose inverse of $A$: For instance
$$
A(V\Delta U^*)A = UDV^*(V\Delta U^*)UDV^*=UD\Delta DV^*=UDV^*=A,
$$
since $D\Delta D=D$. All other properties can be verified similarly.
If $Ax=b$ is solvable, i.e., there is $x_0$ such that $Ax_0=b$, then $A^+b$ is a solution:
$$
A(A^+b) = AA^+Ax_0 = Ax_0=b
$$
by the properties of the Moore-Penrose inverse.
A: Let’s follow your links shall and consider all matrices over a field $K$ of real or complex numbers. 
Let a $m\times n$ matrix $A$ has a singular-value decomposition $A=UDV^*$. Then 
the pseudoinverse $A^+$ of the matrix $A$ is $VD^+U^*$, where $D^+$ is the pseudoinverse of $D$, which is formed by replacing every non-zero diagonal entry by its reciprocal and transposing the resulting matrix (see here). That is, $D^+=\ ^t{\Delta}$. 
Since we have assumed that a linear system $Ax=b$ has (infinitely many) solutions, they are all given by $x=A^+b+[I-A^+A]w$ for arbitrary vector $w$, (see here). That is $x=VD^+U^*b+[I- VD^+U^*A]w$. In particular, if $w=0$ then $x=VD^+U^*b$. 
A: Additional notes to @daw's answer:
If $AX=B$ has infinitely many solutions, the set of solutions is exactly:
$$\Big\{ A^+ B + [I_n - A^+ A]w,\ \rm{ with }\ \ w \in \mathbb{R}^n \Big\}$$
Proof:
As mentioned by @daw, since we assume that there are solutions, there exists a $X_0$ such that $A X_0 = B$.
Now let $X_1 = A^+ B + [I_n - A^+ A]w$, for any vector $w$. Then: 
$$\begin{align}A X_1 &= A A^+ B + A [I_n - A^+ A]w \\ &= A A^+ A X_0 + A w - A A^+ A w \\ &= A X_0 + Aw - Aw \\ &=B  \end{align}$$
so $X_1$ is indeed solution of $A X = B$.
Conversly, let's assume $X_2$ is a solution of $AX=B$. Then:
$$X_2-A^+B=X_2-A^+AX_2 = (I_n - A^+ A) X_2$$ so, with $w=X_2$, we have:
$$X_2=A^+B+ (I_n - A^+ A) w.$$
