Should we always minimize squared deviations if we want to find the dependency of mean on features? It is known that for any distribution a minimization of squared deviations from a value $c$ gives mean of the distribution. In other words, if we generate many values $y_i$ using our distribution, and then try to find $c$ that minimizes the following expression:
$
\lim_{N \to \infty}  \frac{1}{N}\sum_{i}^N (y_i - c)^2
$
then we find that the minimum is achieved when $c = \mu$, where $\mu$ is the mean of the given distribution. And for any $N$ the above minimization gives the sample mean.

Now, let's assume tat we have a distribution which depends on a single real-valued feature $x$, and we have a data-set generated by this distribution: ($x_1$, $y_1$), ($x_2$, $y_2$), $\dots$, ($x_n$, $y_n$). Then we try to find a function that gives dependency of the mean of this distribution on the feature:
$
\mu = \mu(x, p_1, p_2, \dots, p_k)
$
where $p_i$ are model parameters that we want to find.
Obviously, that we can get the desired function (values of the parameters) if we choose them by minimization of squared deviations:
$
\frac{1}{N} \sum_{i}^{N} [y_i - \mu(x_i, p_1, p_2, \dots, p_k)]^2
$
In other words, we can find the correct answer (correct values of the parameters $p_i$) under condition that $N$ goes to infinity and that the real function, that we are searching, is covered by the space of functions that we consider.
To summarize, the minimization of squared deviations can give not only mean of a given constant distribution, it can also gives correct dependency of the mean on a feature.

In reality $N$ is never infinity and can be rather small. As a result minimization of squared deviation will not give exactly correct values of the model parameters $p_i$. For another sample of the same size we will give other estimations of the parameters $p_i$. In general case, we can say that our estimations of the parameters are somehow distributed.
If we, instead of squared deviations minimize some other measure, we will get other distributions of our model parameters $p_i$. So, different measure of accuracy give different distributions in the parameters space for a given fixed sample size.
These distributions are characterized by their "width" and possible a "shift" relative to the correct point (true values of the model parameters).
Can it be the case that a use of something different from squared deviation gives a better estimate of model parameters (for example more accurate (smaller width) and with smaller or no systematic from the correct model parameters?

A similar question (if not the same) is: If the predicted value of machine learning method is E(y | x), why bother with different cost functions for y | x?
The theoretical mean of a distribution minimizes the squared error, but that does not mean that the sample mean is always the best estimator with regards the squared error loss function. The sample mean has a statistical variation.
In an answer to the above question an example is givem that shows how the sample median is performing better than the sample mean when the errors are Laplace distributed. Below is a copy of the image:


*

*Another example question is: Why is the Median Less Sensitive to Extreme Values Compared to the Mean?
The median can be a better estimator in the case of distributions with outliers.


*Related is also: Could a mismatch between loss functions used for fitting vs. tuning parameter selection be justified?
The answer to that question explains how, when we wish to have an estimator that optimizes the mean squared error, then it doesn't mean that we need to use the squared error loss function in fitting/training the model.
Asside from using different estimators like the median or a maximum likelihood estimators, there is also the concept of biased estimators that can improve the expectation of the loss. Examples are regularisation (Ridge regression, lasso regression), Bayesian estimators, shrinking.
