From the BatchNorm paper, section 4.2.3, (https://arxiv.org/abs/1502.03167),
The ensemble prediction was based on the arithmetic average of class probabilities predicted by the constituent networks.
Is there a theoretical basis for doing this? Is the output value after averaging of individual probabilities, still a valid probability?
From the law of total probability we know that for disjoint events $H_n$, we can calculate:
$$P(A) = \sum_n P(A|H_n) * P(H_n)$$
Basically if $P(A|H_n), n:1,...,N$ are different networks emitting probabilities, and $H_n$ is a disjoint hypothesis space then the result is a probability.
When doing simple averaging they are assuming that $P(H_n) = \frac{1}{N}$ for all $n:1,..,N$; a discrete uniform distribution.
The biggest problem with these kind of averages is that nobody really checks if the hypotheses are in fact disjoint or whether it makes sense to assign equal probabilities to each or not. Hypotheses usually end up being very similar to each other. As a result, mathematically speaking the result is still a probability, but from a Bayesian averaging point of view, it is not a well thought prior.
As I already noticed in my comment, you can find a partial answer to your question and further details in my other answer and the references that were provided.
What you seem to be asking, is "how do we know that the average of probability forecasts is a valid probability?", at least this is how I understand it. Your question asks about taking averages of multiple probabilistic forecasts to take pooled forecast, so it is closely resorted to linear opinion pools (Stone, 1961).
First thing to notice is that a probability forecast, is in fact a conditional probability distribution. Taking arithmetic mean is a special case of taking a weighted sum $\sum_k w_k x_k$ with $\forall\, w_k > 0$ and $\sum_k w_k = 1$, where $w_1 = w_2 = \dots = w_n = n^{-1}$, so it is a convex combination. A weighted sum of probability distributions leads to a mixture distribution
$$ p(x) = \sum_k w_k \,p_k(x) $$
where $p_k$ are some probability density (or mass) functions.
As already said by Cowboy Trader, you can think of this in terms of basic laws of probability. Given the properties of weights $w_i$, we can think of them as of probabilities, the most meaningful interpretation would be considering them as prior probabilities for choosing those forecasts. In such case, their joint distribution is
$$ p(x, k) = p_k(x) \,w_k = p(x|k) \, p(k) $$
what follows from the definition of conditional probability. When we have joint distribution, we can calculate marginal distribution of it by the law of total probability
$$ p(x) = \sum_k p(x, k) = \sum_k p(x|k) \, p(k) $$
If you also want to ask "why do people use it?", then the answer is: because it just works.
Yes, there is theoretical basis, and No, we don't know why it works. Look up "forecast combination puzzle" in interweb, e.g. this presentation, p.20. Somehow, a simple average of multiple models appears to outperform single model forecast and weighted average forecasts in practice. There are many hypotheses of why this happens, but there is no consensus in forecasting literature. This could be because an optimal weight in a weighted average combination has too much noise, so in the end a simple average works better
Yes there is a theory for it, it's called ensemble learning. The method of bagging (bootstrap aggregating) relies on it. This is used for example in Random Forests.
The intuitive idea is that by averaging models that have a very low bias but a high variance, you can reduce that variance while still keeping the bias low. This is what happens with random forests where you usually use deep trees that can overfit (i.e. low bias-high variance), but averaging their prediction reduces this overfitting. This of course works best if the training sets of all the models are independent but in practice you use bagging.
In DL models the diversity in the ensemble comes from different hyperparameters: they highlight here different initialization, dropout levels, BN or not.
As for the second part of your question, I think Cowboy Trader answered it best. However, the ensembling also works with outputs that are not probabilities, like for example in the case of regression.
