Your problem might well have infinitely many solutions. It might also have no solutions! You should consider whether or not you would like to find a least squares solution rather than an exact solution, and whether you want to pick out some particular solution if there are many solutions. Both of these could be dealt with in the polynomial optimization approach that I've suggested.

Asymptotically, for large values of $n$ or $m$, you're in deep trouble. For small values of $m$ and $n$ there are some interesting and useful approaches that can work very well in practice. How big are your $m$ and $n$?

Solving rank deficient and discrete ill-posed least squares problems is a very well studied subject. There's an inevitable tradeoff between stability of the computed solution with respect to round-off as well as error in the right hand side (your $N$ matrix) and biasing the solution away from a true least squares solution. By adjusting the regularization term, you can control in what direction the solution is biased by the regularization. By adjusting a regularization parameter you can control the trade-off between stability and bias. See for example Per Hansen's book.

What do you know about $M$? In particular do you know whether it is simply rank deficient or whether it might be ill conditioned? In terms of the eigenvalues of $M$ are the nonzero eigenvalues well separated from the 0 eigenvalues, or could you have small nonzero eigenvalues?

If you want to use this sparse QR "factorization" to solve $\min \| Ax-b \|$, the procedure is to first apply the Givens rotations to reduce $A$ to $R$, possibly reordering the columns of $A$ to minimize fill-in in the $R$ matrix. As you do these Givens rotations you also apply them to $b$ to obtain $Q^{T}b$. Then you solve the triangular system $Rx=Q^{T}b$ to obtain the least squares solution. There is lots of available software for performing these computations.

The typical case in which sparse QR factorization is used is where $A$ is of size $m$ by $n$, $m \gg n$, and $A$ has full column rank. In this situation, $Q$ would be $m$ by $m$ and typically dense with nonzeros, but the effect of multiplication by $Q$ can be achieved by using sequences of Givens rotations. So $Q$ isn't stored explicitly. It is also possible to order the columns of $A$ to reduce the nonzero fill-in in the $R$ matrix.

The MATLAB toolbox is finding a local minimum solution- these will occur at points where several of the ellipsoids intersect. However, there's no way to be sure that the local minimum it finds will be a global minimum. The comment about "smallest sphere touching the intersection" is inconsistent with the original statement of the problem.

This particular lower bound might be quite useful within a branch and bound algorithm for the problem.

I don't know of any code for this specialized problem. You could certainly give it to a more general purpose branch and bound code for non-convex (MI)NLP problems like BARON. Using such a solver (or a custom program written by you), it should be possible to get reasonably good solutions with bounds (e.g. "Here's a solution with objective value 21.72, and our best bound on the optimal value is 21.45.") within a few minutes of computation. The process will be much faster if you have a reasonably good (in objective value) feasible solution to start with.