What is the point of having a dense layer in a neural network with no activation function? I am surprised by the default setting in keras:
keras.layers.Dense(units, activation=None, ...)

Why do we have the option of only using a dense layer (which is matrix multiplication) but without an activation function (non-linear transformation)? I think these two should always go together in a neural network. Is there another case where we can use a dense layer without an activation function?
 A: Another purpose of using linear layers is to reduce dimensionality (and the number of parameters). For example the Skip-gram and CBOW model for word embeddings.

The training task is to predict context words of a given word $p(w_o|w_i)$. A naive way is to count the occurrences of context words for each word and put them in a matrix $M$, then the probability is just $p(w_o|w_i) = f(M_{w_i})_{w_o}$, where $f$ is a normalization function.
The problems is often the number of words $n$ is huge and we can't afford an $n$ by $n$ matrix. So we can first use an $n$ by $d$ matrix $A$ to reduce the dimension (to say 128) and use another $d$ by $n$ matrix $B$ to turn it back, then number of parameters can be reduced to $2*d*n$. It's kind of like matrix decomposition in the sense that we use $BA$ to approximate $M$. 
The model can be implemented as two linear layers followed by a normalization function and can be trained using the cross-entropy loss $E[y_n\log f(BAw_i)]$, where for the Skip-gram model, $w_i$ is a one-hot vector, for the CBOW model, $w_i$ is a BOW vector.
A: Suppose you have a network for either a binary classification task. One way to implement this is using a final activation that yields predicted probabilities (non-negative, sum to 1). An inverse logistic function $f(x)=\frac{1}{1+\exp(-x)}$ is one way to do this in the binary case.
Using the standard binary cross-entropy loss $\mathcal{L}=-\sum_i y_i \log(\hat{y}_i)+(1-y_i)\log(1-\hat{y}_i)$ means that we're round-tripping exponential functions and logarithmic functions. This can cause a severe loss of precision due to accumulated numerical error.
On the other hand, working on the scale of $x$ instead of the probability scale eliminates the round-tripping. In other words, using an identity function the final layer, and a loss function that works on the scale of the linear predictor, achieves the same model without loss of precision.
We can extend the same reasoning to the case of a $k$-nary classification task with a softmax activation.
See: Numerical computation of cross entropy in practice
A: One such scenario is the output layer of a network performing regression, which should be naturally linear. This tutorial demonstrates this case.
Another case that comes to my mind are deep linear networks which are often being used in neural networks literature as a toy model for studying some phenomena that would be too complex with usual non-linear networks.
A: If you choose to use activation=None, you for example add a BatchNormalization layer before you actually use the activation. This is used often in convolutional neural networks, but is good for dense neural networks as well.
z = tf.keras.layers.Dense(20, activation=None)(z)
z = tf.keras.layers.BatchNormalization()(z)
z = tf.keras.layers.Activation("relu")(z)

A: You always want the most flexibility that is possible out of the library that you are using. For example, if you want to have a deep NN with skipped connections ( see this paper ) you need to apply your activation function after the operation F(x) + x is done. how can you implement that on a dense layer with the output F(x) that has no option to stop it from applying activation function before you can do the summation?
