# Solving linear simultaneous equations statistically

I have lots of simultaneous equations, all containing a few variables from a set of a large number of variables. None of the equations are exact, it's possible that some are completely incorrect. Some equations may conflict with other equations, so it is unlikely that there is an exact solution. Is there a standard way to find the "best solution" to the system of simultaneous equations? For example a least squares solution or something similar?

Thanks

Suppose the system of equations $Ax = b$ is overdetermined, so there is no solution. What do we do?

One thing we can do is find $x$ so that the residual $b - Ax$ is as "small" as possible. To be specific, we can find $x$ to minimize the 2-norm of the residual. This $x$ is called the "least squares" solution to our system of equations. (We could also use the 1-norm, or some other norm, but in least squares we use the 2-norm.)

Here is some geometric intuition. For any $x$, $Ax$ is a point in the column space of $A$. $b$ is not in the column space of $A$, but we can choose $x$ so that $Ax$ is the vector in the column space which is as close as possible to $b$. Visually, you can see that the residual $b - Ax$ is then orthogonal to the column space of $A$. In particular, $b - Ax$ is orthogonal to each column of $A$. This tells us that \begin{equation} A^T(b - Ax) = 0 \end{equation} or in other words \begin{equation} A^T Ax = A^T b. \end{equation} This system of equations is called the "normal equations", and by solving it you can find a least squares solution to the system $Ax = b$.

You could solve the normal equations with whatever method you like, but numerical linear algebraists have carefully studied the problem of computing least squares solutions, and have come up with special methods (based on QR decomposition or SVD) that are efficient and insensitive to the effects of roundoff error. Trefethen's book Numerical Linear Algebra covers this topic nicely.

If you're working in Matlab, the command

x = A\b ;


gives you a least squares solution to $Ax = b$.

By the way, when you say that some of your equations might be completely incorrect, this suggests that it might be better for you to choose $x$ to minimize the 1-norm of the residual. The 1-norm doesn't mind if a few components of the residual are huge, whereas the 2-norm cannot stand it. The 1-norm won't try so hard to satisfy those equations which are completely incorrect. Finding $x$ to minimize $\|b - Ax\|_1$ can be formulated as a linear program. You could find an explanation in the book Convex Optimization by Boyd and Vandenberghe, which is free online. In Matlab, you could use cvx to minimize $\|b - Ax\|_1$ in just a few lines of code.