Why is using squared error the standard when absolute error is more relevant to most problems? I recognize that parts of this topic have been discussed on this forum. Some examples:


*

*Is minimizing squared error equivalent to minimizing absolute error? Why squared error is more popular than the latter?

*Why square the difference instead of taking the absolute value in standard deviation?
What I still don't understand is why OLS regression is still the default solution to the problem of linear regression. It seems to me that in the majority of practical situations, the costs associated with errors are linear or approximately linear. If I order 2 extra parts, I incur twice the unnecessary cost as compared to if I order 1 extra part. Therefore, the optimal solution that OLS produces will not correspond to an optimal solution in reality. I understand that in general, you can minimize whatever error metric makes the most sense in the situation. My question is not about whether it is possible or a good idea to use MAE in a specific case; it is about the convention. Why is MSE minimized in nearly all simple cases instead of MAE when the real cost is typically linear? 
The cases I have seen made for minimizing MSE include the following:


*

*MSE is continuously differentiable

*Squaring gives a greater penalty to large errors

*It makes sense because we assume errors to be normally distributed


However, if we can perform regression with absolute errors easily, and we are concerned primarily with prediction, I don't see how those points lead anyone to choose squared errors. If we are to choose a convention, isn't absolute error better? 
This post has received many excellent answer which have all been useful to me. Of those answers, and the answers elsewhere which the mods believe answer my question, none of them exactly address the real source of my confusion except for the answer by @richard-hardy.
 A: I think it's worth taking a step back and considering what the two losses imply.
Looking at it from a probabilistic point of view, the loss function is equivalent to the assumed log-likelihood function and thus should correspond to how we think our measurements are distributed around their unknown 'true' values.
As you say, in the case of OLS this is equivalent to assuming a Gaussian likelihood, where as an absolute error loss function is equivalent to a Laplacian likelihood. Gaussian likelihoods are far far more often a good match to real life as a consequence of the central limit theorem.
Our predictions are in general improved by making our assumed (and implicitly generative) model as close a match to reality as possible. In many (most?) cases this will improve the predictive accuracy by any sensible metric (including e.g. mean absolute error). It is far more often the case assuming a Gaussian likelihood will achieve this.
A: If errors are independent and follow the normal distribution (of any variance but consistent), then the sum of squared errors corresponds to their joint probability/likelihood.
$\Pi e^{-x_i^2}=e^{-\Sigma x_i^2}$
So under those conditions minimizing the sum of square errors is the same as maximizing the likelihood.

If a cost-minimizing prediction is needed (where the cost metric is different from MSE) the general/accurate approach would be to explicitly minimize the expected cost over the entire distribution of models weighted by their likelihoods (or probabilies if you have prior knowledge). This completely decouples the problem of minimizing expected cost from the problem of estimation in the presence of noise.
Suppose you are measuring a constant quantity in the presence of Gaussian noise. Even if your cost metric for future outcomes is MAE, you would rather predict with the mean (minimizing past MSE) than the median (minimizing past MAE), if indeed you know the quantity is constant and the measurement noise is Gaussian.
Example
Consider the following spread of hits produced by a gun that was mechanically fixed in place. You place a circle of a given size somewhere on the target. If the next shot lands entirely inside your circle, you win, else you lose. The cost function is of the form $f_C(x,y)=sign((x-x_C)^2+(y-y_C)^2-R^2)$.

If you minimize $\sum_i f_C(x_i,y_i)$, you would place the circle in the blue position, containing entirely the maximum number of past shots.
But if you knew that the gun is fixed in place and the error is Gaussian, you would place the circle in the green position, centered on the data's mean/centroid (minimizing MSE), as you are optimizing future expected payoff, not average past payoff.
A: TLDR; when nothing is known about actual cost of error to the user of the model, MSE is a better default option compared to MAE because, in my opinion, it is easier to manipulate analytically and is more likely to match the actual cost of error.
It's a great question. I like that you start with desire to make your loss function match actual costs. This is how it's supposed to be done ideally in my opinion. However, it is impractical to derive the cost function from actual costs every time you build a model, so we tend to gravitate to using one of the loss functions available in software. Least squares is one of the most popular functions mainly due to mathematical convenience. It is easier to deal with it analytically. Also, in some cases least squares produces unbiased point forecast, that is $E[y]-\hat y=0$, which is often considered desirable for sentimental reasons.
Having said this, I must argue that it is not obvious to me that absolute value loss is more realistic. Consider, drug overdoses - they are much costlier than underdoses in some situations: not getting high enough vs dying. 
Within your parts example, consider this: what if you underestimated the cost of parts to be \$1, and entered into a forward agreement to deliver one million parts one month later at \$1.1 knowing that you will have $1M one month from today. You are going to make 10% profit!
Then comes the day and parts are actually $1.2 a piece. So, you are not only going to incur loss of \$100K, but will also lack funds to deliver 1M parts. So, you are forced to default and go into bankruptcy which is very expensive. On the other hand if you overestimated the cost of parts, then you wold forego some profit but wouldn't end up in dire situation of insolvency or liquidity crisis.
This is a very common situation in business where the losses are asymmetric and highly nonlinear with rapidly escalating costs in one direction of forecast error but not the other. 
Hence, I'd argue that absolute loss, which is symmetric and has linear losses on forecasting error, is not realistic in most business situations. Also, although symmetric, the squared loss is at least non linear.
Yet the differences between absolute and squared loss functions don't end here. For instance, it can be shown that the optimal point forecast in absolute loss is the median while for the squared loss it is mean. 
I think that the following loss function is more suitable to business forecasting in many cases where over forecasting error $e=y-\hat y$ can become very costly very quickly: $$\mathcal L(e,\hat y)=|\ln\left(1+\frac e {\hat y}\right)|$$
Here, if you are forecasting a non negative quantity $y$, then over forecasting is potentially devastating. Imagine you are bank forecasting the deposit volume, and the actual deposit volume turned out to be much lower than you hoped for. This can have severe consequences. This type of asymmetric loss function will lead to a biased optimal point forecast, i.e. $E[y]-\hat y\ne 0$, but that's exactly what you want: you want to err on the side of under forecasting in this kind of business problem.
A: The first 5 answers fail to distinguish between estimation loss1 and prediction loss2, something that is crucial in answering the question. A priori, there is no reason that the two should coincide. I will discuss both types of loss in the context of point prediction using linear regression. The discussion can be extended to models other than linear regression and tasks other than point prediction, but the essence remains the same.
Setup
Suppose you are facing a prediction problem where the model is
$$
y=X\beta+\varepsilon
$$
with $\varepsilon\sim D(0,\sigma)$, $D$ being some probability distribution with location $0$ and scale $\sigma$. You aim to predict $y_0$ given $x_0$, and your point prediction will be $\hat y_0$, a function of $x_0$, the data sample, the model and the penalty (the negative of reward) function defined on the prediction error. The penalty function you are facing is $L_P(y-\hat y)$. It has a minimum at zero (the value $L_P(0)$ can be set to zero without loss of generality) and is nondecreasing to both sides of zero; this is a typical characterization of a sensible prediction loss function. You can freely choose an estimation loss function $L_E(\cdot)$ and a point prediction function $y_hat_0$. What are your optimal choices for each? This will depend on the error distribution $D$ and the prediction loss function $L_P(\cdot)$.
Estimation loss
Estimation loss specifies how parameter estimates of a model are obtained from sample data. In our example of linear regression, it concern the estimation of $\beta$ and $\sigma$. You can estimate them by minimizing the sum of squared residuals (OLS) between the actual $y$ and the corresponding fitted values, sum of absolute residuals (quantile regression at the median) or another function. The choice of the estimation loss can be determined by the distribution of model errors. The most accurate estimator in some technical sense* will be achieved by the estimation loss that makes the parameter estimator the maximum likelihood (ML) estimator. If the model errors are distributed normally ($D$ is normal), this will be OLS; if they are distributed according to a Laplace distribution ($D$ is Laplace), this will be quantile regression at the mean; etc.
*To simplify, given a ML estimator, you may expect more accurate parameter estimates from your model than provided by alternative estimators.
Prediction loss
Prediction loss specifies how prediction errors are penalized. You do not choose it, it is given. (Usually, it is the client that specifies it. If the client is not capable of doing that mathematically, the analyst should strive to do that by listening carefully to the client's arguments.) If the prediction error causes the client's loss (e.g. financial loss) to grow quadratically and symmetrically about zero, you are facing square prediction loss. If the client's loss grows linearly and symmetrically about zero, you are facing absolute prediction loss. There are plenty of other possibilities for types of prediction loss you may be facing, too.
Prediction
Given the parameter estimates of the model and the values of the regressors of the point of interest, $x_0$, you should choose the point prediction $\hat y_0$ based on the prediction loss. For square loss, you will choose the estimated mean of $y_0$, as the true mean minimizes square loss on average (where the average is taken across random samples of $y_0$ subject to $x=x_0$). For absolute loss, you will choose the estimated median. For other loss function, you will choose other features of the distribution of $y_0$ that you have modelled.
Back to your question
Why do people frequently choose square error rather than absolute error, or correspondingly square loss rather than absolute loss, as estimation loss? Because normal errors ($D$ being normal) are common in applications, arguably more so than Laplace errors ($D$ being Laplace). They also make the regression estimators analytically tractable. They are not much easier to compute, however. Computational complexity of OLS (corresponding to ML estimation under normal errors) vs. quantile regression at the median (corresponding to ML estimation under Laplace errors) are not vastly different. Thus there are some sound arguments for the choice of OLS over quantile regression at the median, or square error over absolute error.
Why do people choose square error, or correspondingly square loss, as prediction loss? Perhaps for simplicity. As some of the previous answers might have mentioned, you have to choose some baseline for a textbook exposition; you cannot discuss all possible cases in detail. However, the case for preferring square loss over absolute loss as prediction loss is less convincing than in the case of estimation loss. Actual prediction loss is likely to be asymmetric (as discussed in some previous answers) and not more likely to grow quadratically than linearly with prediction error. Of course, in practice you should follow the client's specification of prediction loss. Meanwhile, in casual examples and discussions where there is no concrete client around, I do not see a strong argument for preferring square error over absolute error.
1 Also known as estimation cost, fitting loss, fitting cost, training loss, training cost.
2 Also known as prediction cost, evaluation loss, evaluation cost.
A: I think the reason is more sociological that statistical.
Short version: We do it this way because we always have.
Longer version:
Historically, we could not do many of the things we now take for granted. Many things are computer intensive and Ronald Fisher was born before Alan Turing. 
So, people did OLS regression - a lot. And people read those regressions in all sorts of substantive fields and statistics courses in those fields taught ANOVA/regression and not more modern methods.
Additionally, editors of journals learned those methods and not others, and many will reject articles with modern methods because e.g. "they won't be understood". 
Many practitioners reject modern methods too; I used to be a sort of data analysis geek at a hospital. Doctors would come to ask my advice and, if it wasn't "do OLS regression" or "do logistic regression" they would reject my advice.
I got my PhD in psychometrics and many of my professors in other branches of psychology did not know any modern methods (one said: "just report the p value, that's what matters"). 
A: Suppose one rolls one die (numered 1-6), and wants to compute its average deviation from the average value of 3.5.  Two rolls would differ by 0.5, two by 1.5, and two by 2.5, for an average deviation of 1.5.  If one takes the average of the squares of the values, one would have one deviation of 0.25, one of 2.25, and one of 6.25, for an average of 2.916 (35/12).
Now suppose instead of rolling one die, one rolls two.  The average deviation would be 1.94 (35/18), and the average square of the deviation would be 5.833 (70/12).
If instead of rolling two dice, one wanted to estimate the expected deviation based upon what it was with one die, doubling the linear average single-die deviation (i.e. 1.5) would yield a value of 3, which is much larger than the actual linear average deviation of 1.94.  On the other hand, doubling the average square of the deviation when using a single die (2.916) would yield precisely the average square of the deviation when using two dice.
In general, the square root of the average of the squares is a more useful number than the average of the squares itself, but if one wants to compute the square root of the average of a bunch of squares, it's easier to keep the values to be added as squares, than to take the square roots whenever reporting them and then have to square them before they can be added or averaged.
A: In my opinion, it boils to that the squared error guarantees a unique solution, easier to work with and hence much more intuition.  By only two main assumptions (and linearity of the error term), a quadratic loss function guarantees that the estimated coefficient is the unique minimized.  Least-absolute deviations does not have this property.  There is always a potential for an infinite number of solutions.  Assuming that $\exists\theta_o\in\Theta$ such that $E(y|x)=m(x,\theta_o)$ and $E((m(x,\theta)-m(x,\theta_o)^2)>0$ for all $\theta\neq\theta_o$, then $\theta_o$ is the unique minimizer for non-linear least squares.
Proof: Let $y=m(x,\theta_o)+u$ and $E(u|x)=0$.  Then $$E_{\theta_o}((y-m(x,\theta))^2)=E_{\theta_o}((y-m(x,\theta_o)+m(x,\theta_0)-m(x,\theta))^2)$$
$$=E_{\theta_o}(u^2)+E_{\theta_o}((m(x,\theta_o)-m(x,\theta))^2)+2E_{\theta_o}(u(m(x,\theta_o)-m(x,\theta))).$$
By the law of iterated expectations, the third term is zero.  Therefore
$$E_{\theta_o}((y-m(x,\theta))^2)=u^2+E_{\theta_o}((m(x,\theta_o)-m(x,\theta))^2)$$
is uniquely minimized at $\theta_o$.  
Another nice property is the total law of variance
$$Var(Y)=Var_X(E_Y(Y|X))+E_X(Var_Y(Y|X)),$$
which can be read as the variance of the dependent variable is the variance of the fitted value plus the residual's variance.
On a more technical note, the asymptotic formulas are much easier for a quadratic loss function.  Importantly, the formulas don't depend on the probability density of the error term.  Unfortunately, that is not true for least-absolute deviations.  Therefore most practitioners end up having to assume independence of the error term (the formula has the conditional density of the error term at 0 conditioned on $x$, which is impossible to estimate($f_{u|x}(0)$)) to estimate $f_u(0)$.  
And the least rigorous point is that people have an easy time understanding what a mean or expected value is, and the quadratic loss solves for the conditional expectation.  Least-absolute deviations soles for the median, which is just harder to interpret.  Another reason quantile regressions aren't very popular.
