Why GLMs predict the mean and not the mode? Why does a GLM predict the mean and not the mode of a signal? Doesn't this contradict the very basis behind the GLM, i.e., maximum likelihood?  The equations to solve for the model parameters in a GLM are based on maximization of the likelihood as described by the probability distribution of the modeled signal. This probability distribution is maximum for the mode not for the mean (the normal distribution is an exception: both mode and mean are the same). Therefore, a GLM should predict the mode, not the mean of a signal! (For some background for this question see here.)
 A: Thanks for all the comments and answers. Although in none of them is 100% the answer to my question, all of them helped me to see through the apparent contradiction. Thus, I decided to formulate the answer myself, I think this is a summary of all ideas involved in the comments and answers:
Maximization of likelihood through the data PDF $f(y; \theta, \phi)$ in GLMs is not related to the mode of $f$ (but to its mean) because of 2 reasons:


*

*When you maximize $f(y; \theta, \phi)$ you do not consider $f$ as a function of $y$, but as a function of $\boldsymbol\beta$ (the parameters of the linear model). More specifically, when you differentiate $f$ to obtain a system of equations leading to determine $\boldsymbol\beta$, you do not do it with respect to $y$; you do it with respect to $\boldsymbol\beta$. Thus, the maximization process gives you the $\boldsymbol\beta$ that maximizes $f$. An optimal $\boldsymbol\beta$, and not an optimal $y$ (which, indeed, would be the mode), is the output of the maximization process.

*Additionally, in the maximization process, the mean, $\boldsymbol\mu$, is a function of $\boldsymbol\beta$. Therefore, through the maximization process we also obtain the optimal $\boldsymbol\mu$.
A: There are two things to argue here:


*

*The facts that a glm attempts the predict $y$ as the mean of a conditional distribution, and estimates its parameters $\beta$ by maximum likelihood are consistent.

*Estimating the parameters by maximum likelihood is not determining the mode of any distribution.  At least not in the classical formulation of a glm.


Lets take the simplest non-trivial glm as a working example, the logistic model.  In logistic regression we have a response $y$ which is 0, 1 valued.  We postulate that $y$ is bernoulli distributed conditional on our data
$$ y \mid X \sim Bernoulli(p(X)) $$
And we attempt to estimate the mean of this conditional distribution (which in this case is just $p$) by linking it to a linear function of $X$
$$ \log\left(\frac{p}{1-p}\right) = X \beta $$
Pausing and reflecting, we see in this case that it is natural to want to know $p$, which is a  mean of a conditional distribution.
In the glm setup, $p$ is not estimated directly, it is $\beta$ that the estimation procedure targets.  To get at $\beta$ we use maximum likelihood.  The probability of observing a datapoint $y$ from the conditional bernoulli distribution, given the value of $X$ observed, and a specific set of parameters $\beta$ ,is
$$ P \left( y \mid X, \beta \right) = p^y (1-p)^{1-y} $$
where $p$ is a function of $\beta$ and $X$ through the linking relationship.
Notice that it is $y$ that is sampled from a probability distribution here, not beta.
To apply maximum likelihood, we flip this around into a function of $\beta$, considering both $X$ and $y$ as fixed and observed:
$$ L(\beta) =  p^y (1-p)^{1-y} $$
But, $L$ is not a density function, it is a likelihood.  When you maximize the likelihood you are not estimating the mode of a distribution because there simply is no distribution to, well, mode-ize.
You can produce a density from $L$ by providing a prior distribution on the parameters $\beta$ and using Bayes's rule, but in the classical glm formulation, this is not done.
A: The goal of maximum likelihood fitting is to determine the parameters of some distribution that best fit the data - and more generally, how said parameters may vary with covariates. In the case of GLMs, we want to determine the parameters $\theta$ of some exponential family distribution, and how they are a function of some covariates $X$. 
For any probability distribution in the overdispersed exponential family, the mean $\mu$ is guaranteed to be related to the canonical exponential family parameter $\mathbf{\theta}$ through the canonical link function, $\theta = g(\mu)$. We can even determine a general formula for $g$, and typically $g$ is invertible as well. If we simply set $\mu = g^{-1}(\theta)$ and $\theta = X\beta$, we automatically get a model for how $\mu$ and $\theta$ vary with $X$, no matter what distribution we are dealing with, and that model can be easily and reliably fit to data by convex optimization. Matt's answer shows how it works for the Bernoulli distribution, but the real magic is that it works for every distribution in the family.
The mode does not enjoy these properties. In fact, as Cliff AB points out, the mode may not even have a bijective relationship with the distribution parameter, so inference from the mode is of very limited power. Take the Bernoulli distribution, for example. Its mode is either 0 or 1, and knowing the mode only tells you whether $p$, the probability of 1, is greater or less than 1/2. In contrast, the mean tells you exactly what $p$ is.
Now, to clarify some confusion in the question: maximum likelihood is not about finding the mode of a distribution, because the likelihood is not the same function as the distribution. The likelihood involves your model distribution in its formula, but that's where the similarities end. The likelihood function $L(\theta)$ takes a parameter value $\theta$ as input, and tells you how "likely" your entire dataset is, given the model distribution has that $\theta$. The model distribution $f_\theta(y)$ depends on $\theta$, but as a function, it takes a value $y$ as input and tells you how often a random sample from that distribution will equal $y$. The maximum of $L(\theta)$ and the mode of $f_\theta(y)$ are not the same thing.
Maybe it helps to see the likelihood's formula. In the case of IID data $y_1,y_2,\ldots,y_n$, we have
$$L(\theta) = \prod_{i=1}^n f_\theta(y_i)$$
The values of $y_i$ are all fixed - they are the values from your data. Maximum likelihood is finding the $\theta$ that maximizes $L(\theta)$. Finding the mode of the distribution would be finding the $y$ that maximizes $f_\theta(y)$, which is not what we want: $y$ is fixed in the likelihood, not a variable.
So finding the maximum of the likelihood function is not, in general, the same as finding the mode of the model distribution. (It is the mode of another distribution, if you ask an objective Bayesian, but that's a very different story!)
