This came up twice now recently, so I'm writing it out as a reference.
Suppose you have a system:
A x = b
where A is n by n, x is n by 1 and b is n by 1. However, rank(A) = m < n. This means there is not just one unique solution, but a family of solutions. Now, let's assume we can build a basis for the null space of A. That is we can write an n by n-m matrix N such that
A N u = 0
for any n-m by 1 vector u. Then we can also write the matrix M = N⊥ which spans the space orthogonal to N. Then we can write any n by 1 vector in these bases:
x = N u + M v
Substituting this into our original system, we see that the choice of u doesn't matter
A (N u + M v) = b
A M v = b
Now we have a system of n equations and m unknowns. This is over-constrained. Let's choose the minimum in the least squares sense.
minimize |A M v - b|^2
Differentiating with respect to v and setting to zero we get
M' A' A M v = M' A' b
which we can solve for v.
Then we can recover x as
x = M v
We can show that indeed this satisfies the original equation by substituting this back in above:
A M v = b
A M ( M' A' A M ) \ M' A' b = b
Let Z = A M, then write this as
Z (Z' Z) \ Z b = b
But (Z' Z) \ Z is simply the pseudo-inverse of Z so we know that Z (Z' Z) \ Z = I and we have that b = b.