Why do we use a Least Squares fit? I've been wondering for a while now if there's any deep mathematical or statistical significance to finding the line that minimizes the square of the errors between the line and the data points. 
If we use a less common method like LAD, where we just consider the absolute deviation, then outliers make less difference to the final model, while if we take the cube of the error (or any other power higher than 2), then outliers are far more significant than with the least squares model.
I suppose what I'm really asking is mathematically, is raising the error to the power of 2 really that special. Is it say more "accurate" in some sense than raising the error to the power of 1.95 or 2.05??? 
Thanks!
 A: Read section 5.14 Why Least Squares? of Meyer's Matrix Analysis and Applied Linear Algebra.
A: I was researching this myself and I found the following answet to be most convincing, atleast in justifying why square but not the absolute values:
If the random scatter follows a Gaussian distribution, it is far more likely to have two medium size deviations (say 5 units each) than to have one small deviation (1 unit) and one large (9 units). A procedure that minimized the sum of the absolute value of the distances would have no preference over a curve that was 5 units away from two points and one that was 1 unit away from one point and 9 units from another. The sum of the distances (more precisely, the sum of the absolute value of the distances) is 10 units in each case. A procedure that minimizes the sum of the squares of the distances prefers to be 5 units away from two points (sum-of-squares = 50) rather than 1 unit away from one point and 9 units away from another (sum-of-squares = 82). If the scatter is Gaussian (or nearly so), the curve determined by minimizing the sum-of-squares is most likely to be correct.
Courtesy: Graphpad.com
A: Consider $n (x_1, x_2,\ldots,x_n)$ measurements which have a normal distribution. Then the probability that only one of the measurements will occur is given by $P(X=x_i) = e^{(x_i - x_T)^2/h^2}$ - where $x_T$ is the true value of the variable $x$ to be measured. (I am ignoring the root $\pi$ factor here because it occurs for all the errors and is not central to our argument.). Now the probability that all of the measurements will occur in an experiment is given by the product: $P(X=x_1)\cdot P(X=x_2)\cdot P(X=x_3)\cdots P(X=x_n)$- (assuming that the measurements are independent of each other - i.e the error from one measurement won't be carried over to the other - which would be the case if the experiment was well designed). The resulting probability is given by:
$$P(x)= e^{(d_1^2 + d_2^2 + \cdots + d_n^2)/h^2} \text{ where }d_i = x_i - x_T$$
Now our aim is to find the value of $x$ for which the above probability is maximum. This would be the true value of the measurement. The above probability is maximum when, the exponent is minimum (remember the negative sign of the exponent) and the exponent is nothing but the sum of the squares of the deviations of the measurements from the true value. This is the idea behind sum of least squares.
We could use calculus:
$$\frac{dP(x)}{dx} = 0$$
If you solve the resulting equation you find that the required value of x is 
$(x_1 + x_2 + x_3 +\cdots+x_n) / n$ - the arithemetic mean of the measurements. That's why we use the AM so much. In a reasonably well designed experiment, where the probability of small errors is large and of large errors is small and positive and negative errors occur with the same probablity (which is when you use a normal distribution), the AM is the most probable true value.
But remember that, the least squares method is only applicable when the measurements can be assumed to have a normal distribution. In other cases, the least squares method does not give the most probable true value.
A: The answers provided to this point have considerable limitations and no answer in favor of squares is clear and convincing.  Some are just dogmatic statements (squares are recommended).  Others fail to relate the question to human life. In real life, we aspire to reduce errors, not their square.  If squares are said to better predict the population, that is tautological if the standard of "prediction" uses squares rather than differences in its definition.  What we need to know is whether squares predict differences better than differences predict differences.  When two people predict the time outcomes of a race or the scores of a football match, who in their nutty mind would square their errors to compare accuracy?  Appealing to "N dimensional space" does not address this.  If we need to predict city water needs in multiple locations, and minimize errors, do you want to minimize the total water shortage, or for some odd reason treat differently the places that were poorly estimated (by squaring the errors)? that is to arbitrarily treat previous extreme observations differently that other cases in the study; if that is what you want to do, do so openly.  Is there seriously real evidence that squares predict population differences better than differences predict population differences?  That is what we need to know and that is what would answer the question.  Gorand suggests that is not what Fisher estimated.  Gorard notes that several fields are using least absolute regression instead of least squares.  See Gorard, S. (2013) 'The possible advantages of the mean absolute deviation 'effect' size.', Social research update., 65 (Winter 2013). pp. 1-4.
and Revisiting a 90-year-old debate: the advantages of the mean deviation.  Stephen Gorard
A: Carl Gauss (the most famous person to live on earth in the 19th century, except for people who did not work in the physical and mathematical sciences) showed that least squares estimates coincide with maximum-likelihood estimates when one assumes independent normally distributed errors with $0$ mean and equal variances.
POSTSCRIPT four years later:
Here are a couple of other points about raising errors to the power $2$ instead of $1.95$ or $2.05$ or whatever.

*

*The variance is the mean squared deviation from the average.  The variance of the sum of ten-thousand random variables is the sum of their variances.  That doesn't work for other powers of the absolute value of the deviation.  That means if you roll a die $6000$ times, so that the expected number of $1$s you get is $1000$, then you also know that the variance of the number of $1$s is $6000\times\frac 1 6\times\frac 5 6$, so if you want the probability that the number of $1$s is between $990$ and $1020$, you can approximate the distribution by the normal distribution with the same mean and the same variance.  You couldn't do that if you didn't know the variance, and you couldn't know the variances without additivity of variances, and if the exponent is anything besides $2$, then you don't have that.  (Oddly, you do have additivity with the $3$rd powers of the deviations, but not with the $3$rd powers of the absolute values of the deviations.)


*Suppose the errors are not necessarily independent but are uncorrelated, and are not necessarily identically distributed but have identical variances and expected value $0$.  You have $Y_i = \alpha + \beta x_i + \text{error}_i$.  The $Y$s and $x$s are observed; the $x$s are treated as non-random (hence the lower-case letter) the coefficients $\alpha$ and $\beta$ are two be estimated.  The least-squares estimate of $\beta$ is $$\widehat\beta = \frac{\sum_i (x_i-\bar x)(Y_i-\bar Y)}{\sum_i(x_i-\bar x)^2} \tag 1$$ where $\bar x$ and $\bar Y$ are the respective averages of the observed $x$ and $Y$ values.  Notice that $(1)$ is  linear in the vector of observed $Y$ values.  Then among all unbiased estimators of $\beta$ that are linear in the vector of $Y$ values, the one with the smallest variance is the least-squares estimator.  And similarly for $\widehat\alpha$.  That is the Gauss–Markov theorem.
A: One very practical reason for using squared errors is that we are going to want to minimize the error, and minimizing a quadratic function is easy - you just differentiate it and set the derivatives to zero, which results in a linear equation - which we have centuries of tricks to help us solve.
I'll walk through a simple example: finding the best line through the origin that fits the data points $(y_i,x_i)$ for $i=1,\dots,n$. Our model for the data is
$$y_i = ax_i + \varepsilon_i$$
where $\varepsilon_i$ is the error of the approximation for the $i$th data point. Let's raise each error to the power $2$ and then add them all up:
$$E = \sum_{i=1}^n \varepsilon_i^2 = \sum_{i=1}^n (y_i - ax_i)^2$$
Minimizing this error with respect to $a$, we differentiate and set the derivative equal to zero:
$$\frac{\partial E}{\partial a} = -2\sum_{i=1}^n x_i (y_i - ax_i) = 0$$
which rearranges to
$$\sum_{i=1}^n x_iy_i = a\sum_{i=1}^n x_i^2$$
so we get the standard least-squares estimator of the slope,
$$a = \frac{\sum_{i=1}^n x_iy_i}{\sum_{i=1}^n x_i^2}$$
If we had raised the errors to any power other than 2 before summing them, the resulting equation would be much harder to solve. Minimizing anything other than the squared error is generally only achievable with a numerical method.
A: Least squares fitting has the desirable property that if you have two different output values for the same input value, and you replace them with two copies of their mean, the least squares fit is unaffected.  For example, the best fit line is the same for the following two sets of data:
0 1
0 5
1 5
2 6

and
0 3
0 3
1 5
2 6

If you use minimum-distance fitting, this is no longer the case.
A: Honestly, it's just really easy to compute.  
While there are other methods that may give better answers in certain situations,  Least Squares with a much simpler algorithm.  It just involves constructing a matrix, where each element is a sum.  Then we can use Row Reduction Echelon Form to find the coefficients.  
This simplicity means that it can be used on cheaper hardware (like your calculator) and it means it's pretty fast on modern cpu's.
On top of that, it's used to find the coefficients of polynomials and rational functions,  two function forms that are extremely fast and easy to compute, and which we can make more accurate by just adding another term.
Honestly, I don't know if it is the theoretical optimal methods, but I do know that practically its one of, if not the easiest to use.
A: Here's a visualization of some of the alternative error functions one could use. 
Clockwise: Absolute error, Square error, Absolute cubed error, Absolute 4th power error, Absolute exponential error. 
Note how the non-square errors exaggerate different error combinations. 

Interactive version here: https://www.desmos.com/calculator/eugeucw7te
A: 16 months after asking the question, I've run into a different and quite physicsy answer, which I hope is useful to someone.
Suppose that in order to determine $m$ parameters of some model:
$$\text{Output} = f(\text{Inputs, parameters})$$
we have conducted $N>m$ experiments. We want to use the information from these experiments to best choose the $m$ parameters (so that the output of the model and the actual experimental value are as close together as possible). 
Now comes the physicsy part: lets construct an N-dimensional phase space so that the $N$ (experimentally determined) outputs from our $N$ experiments are represented by a single point in this space (the Cartesian coordinates of this point are the outputs from each experiment). Call this the 'data-point'.
Secondly, if we choose an arbitrary set of parameters for our model, we can use the inputs for each experiment to construct a 'predicted output' for each experiment (by parsing the inputs through our model). There will be $N$ predicted outputs (one for each experiment) and these form a second point in the phase space, say the 'prediction-point'. As we varying the parameters this point moves about in an $m$-dimensional subspace of the phase space. And this is the important point:

The sum of the squares of the error terms (SSE) is the square of the distance between these two points in the $N$-dimensional phase space, just by Pythagoras' Theorem.

So minimising the sum of squares error is equivalent to minimising the distance between the data-point and the prediction-point in the $N$-dimensional phase space - a very natural way of calibrating our model. 
Finally, from this Gauss' result makes some sense - if the data point is allowed to vary normally with 0 mean and equal variances, the error will be spherically-symmetric around the data-point, and so the closer our prediction-point is to the data-point in space, the better, and minimising this distance should  give the maximum-likelihood estimator.
A: The choice of least squares is often due to familiarity with the method, a herd instinct. Yet there are compelling mathematical reasons to use the 2-norm.
The first reason is the Gauss-Markov theorem. (Meyers, Matrix Analysis and Applied Linear Algebra, 2000, S 5.14). Essentially, the least squares estimators are unbiased and they have minimal variance.
Next, we have continuous differentiation in this norm. We can use derivatives to find the extrema and we know the extrema are minima (Proof that least squares minimizers are a convex set.)
