How to verify solution to pre-conditioned linear systems solver?

You've started with a singular linear system of equations $Ax=b$. As a practical matter, it's unlikely that $b$ lies exactly in the range of $A$, so at best you can find a least squares solution that minimizes

$min | Ax - b |_{2}$

Because the system is singular, the null space of $A$ is non-empty, and there will be an infinite number of solutions to this least squares problem. In particular, if $x_{LS}$ is some particular least squares solution and $v$ is any vector in the null space of $A$, then any $x$ of the form $x=x_{LS}+alpha v$ will be a least squares solution.

Which of these infinitely many least squares solutions is the "correct" one?

Perhaps you don't care and will accept any least squares solution. In that case, any solution that satisfies the normal equations $A^{T}Ax=A^{T}b$ will be correct.

Perhaps you want to simultaneously minimize the norm of $x$ to get a minimum norm least squares solution? You can regularize the least squares problem by minimizing

$min | Ax-b |_{2}^{2} + lambda^{2} | x |_{2}^{2}$

where $lambda$ is a small parameter. This least squares problem is strictly convex, so it has a unique optimal solution.

What happens in floating point arithmetic with limited precision? In that case, it's likely that $A$ will have a very large but still finite condition number.

You can throw a preconditioner $M$ at $A$, in hopes that $M^{-1}A$ will be better conditioned than $A$. This might work in the sense at $M^{-1}A$ is numerically better conditioned than $A$, but since we know that $A$ is theoretically singular, this is a mirage! $M^{-1}A$ should still be singular.

So what is the solution to the numerically well-conditioned system $M^{-1}Ax=M^{-1}b$ doing? This is similar to regularization in the sense that $M^{-1}A$ is better conditioned than $A$, but you don't know much about the direction in which the preconditioner has pushed the solution. If $M^{-1}Ax-M^{-1}b$ is small, and $M^{-1}$ is well behaved, then perhaps $Ax-b$ is small and you've got a least squares solution. On the other hand, if $M^{-1}$ is very badly scaled, it may be that the residual in the original problem $Ax-b$ is quite large.

There has been some research on preconditioners that are designed to achieve a particular regularization. See for example:

Calvetti, Daniela. 2007. “Preconditioned Iterative Methods for Linear Discrete Ill-Posed Problems from a Bayesian Inversion Perspective.” Journal of Computational and Applied Mathematics, Special Issue: Applied Computational Inverse Problems, 198 (2):378–95. https://doi.org/10.1016/j.cam.2005.10.038.