When no solution exists, how do you satisfy as many constraints as possible? Suppose there exists no $x$ such that $Ax=b$. How would you go about choosing $x$ to maximize $|\{i : A_i x = b_i \}|$?
 A: My first hunch would be to create a Mixed Integer Programming model:
$$\begin{array}{ll}
\min & \sum_i y_i \\
    & -My_i \le  \sum_j a_{i,j} x_j-b_i \le My_i \\
    & y_i \in \{0,1\} \\
\end{array}
$$
where $M$ is a large enough constant. There are some other tools developed to help diagnose infeasible Linear Programming models, such as IIS (Irreducable Infeasible Set).
See also the book John W. Chinneck: Feasibility and Infeasibility in Optimization. 
In a different answer it is suggested to add a positive slack to each equation and then minimize the sum of the slacks. This is related but I think not the same. The first thing we should notice is that we need to allow positive and negative slacks. We then can either minimize the absolute values or the squared values. I.e.:
$$\begin{array}{ll}
\min & \sum_i |s_i| \\
    & Ax = b + s \\
\end{array}
$$
or
$$\begin{array}{ll}
\min & \sum_i s^2_i \\
    & Ax = b + s \\
\end{array}
$$
The latter is of course just ordinary least squares. Besides this technical issue there is a more fundamental difference. When minimizing slacks we minimize the sum of of the infeasibilities. This is very different from minimizing the number of infeasibilities. Minimizing the sum of the infeasibilities may lead to a large number of small infeasibilities. In my MIP formulation we explicitly minimize the number of the infeasibilities.
In a reaction we see two objections made to the MIP formulations:


*

*How to choose $M$. This is valid point. We need in advance an upper bound on the maximum infeasibility we expect. This is a major drawback.

*This method would not flag things exactly. The concept of exactness is somewhat difficult when solving a system of linear equations. At least in LP/MIP solvers there are a large number of tolerances involved. But even a straight equation solver will not find exact solutions unless very special techniques are used. I don't think the MIP approach adds much "inexactness". 

A: Rahul in the comments made an excellent point that your problem is to minimize the $L_0$ norm (which should be mentioned doesn't satisfy all properties of a norm). The L0 norm is the number of non-zero elements. 
$$ \text{minimize (over $x$) } \|Ax - b\|_0 $$
This is not a convex problem, and may be difficult to solve. An alternative is the $L_1$ norm, which is a well defined norm and can be minimized efficiently on a computer. The $L_1$ norm may be a reasonable approximation of the $L_0$ norm.
$$ \text{minimize (over $x$) } \|Ax - b\|_1 $$
Moving perhaps further afield from your original objective, another norm to minimize is the 2-norm which would give the solution in the least squares sense.
$$ \text{minimize (over $x$) } \|Ax - b\|_2 $$
Minimizing the 2 norm has the well known analytic solution:
$$ x = (A^{\top}A)^{-1}A^{\top} b $$
where $(A^{\top}A)^{-1}A^{\top}$ is the pseudo inverse for a skinny matrix.
A: Another option would be to define a positive penalty $p_i$ for every constraint $i$:
$$
a_ix_i+p_i =b_i
$$
and then to minimize 
$$
\sum_i \omega_ip_i
$$
subject to
$$
a_ix_i+p_i =b_i\\
p_i\ge 0,
$$
where $\omega_i$ is a positive weight associated to $p_i$. If the objective is 0, the solution is feasible. One advantage of this method is that you can control which constraint you want to satisfy in priority (by giving a large $\omega_i$). If you give the same $\omega_i$ to each constraint, you might violate every constraint, but by very little! up to you.
Note. As recommended in some of the comments below, it is a good idea to let $p_i$ be a variable in $\mathbb{R}$ and to minimize the squares, because if you impose $p_i\ge 0$, you will force $a_ix_i<b_i$, whereas it may be more interesting to have $a_ix_i>b_i$.
