Solution (source code)

= Solution

Use the <singular value decomposition> and the <orthogonal invariance of the Euclidean norm>. Put $z=V^Tx$ and $c=U^Tb$. Since precisely $\sigma_1,\ldots,\sigma_r$ are positive,
$$
\|Ax-b\|_2^2=\sum_{j=1}^r(\sigma_jz_j-c_j)^2+\sum_{j=r+1}^m c_j^2.
$$
The second sum is independent of $x$ and cannot be removed. The first sum vanishes exactly when $z_j=c_j/\sigma_j$ for $1\leq j\leq r$, while $z_{r+1},\ldots,z_n$ are arbitrary. Thus every <least-squares solution> is
$$
\boxed{x=\sum_{j=1}^r\frac{u_j^Tb}{\sigma_j}v_j+\sum_{j=r+1}^n\alpha_jv_j,\qquad\alpha_j\in\mathbb R.}
$$
The unavoidable <linear-system residual> is the component of $b$ <orthogonal> to the <column space> of $A$:
$$
\boxed{\min_x\|Ax-b\|_2=\left(\sum_{j=r+1}^m|u_j^Tb|^2\right)^{1/2}.}
$$
The <minimum-norm least-squares solution> sets all $\alpha_j=0$ and equals $A^+b$, where $A^+$ is the <Moore-Penrose inverse>. The squared minimum, rather than the minimum itself, is the second sum in the first display.