Pooled Covariance Matrix with very different amount of samples per class I have a dataset with 10 classes, and want to estimate the covariance. It turns out that due to numerical stabilitiy, it is much better to use a pooled covariance matrix. Suppose I have $N$ samples per class. Then 
$$S_{\mathrm{pooled}} = \frac{1}{10 * N} \sum_{i=1}^{10} \sum_{j=1}^N (x_j - \mathrm{mean}(x_j)) (x_j - \mathrm{mean}(x_j))'$$
I would also like to perform LDA (for dimensionality reduction and later classification) on the dataset, and the computation of the Within-Class Scatter Matrix $S_W$ is almost the same upto the scaling factor.
I do have a very unbalanced classes - quite different number of samples per class. For example for class 1, I have only $100$ samples whereas for class 4 I have $5000$ samples!
Let $N_1, N_2, ... N_{10}$ denote the amount of samples for each class. According to https://en.wikipedia.org/wiki/Pooled_variance, one computes as
$$S_{\mathrm{pooled}} = \frac{1}{N_1 + N_2 + ... N_{10} - 10} \sum_{i=1}^{10} \sum_{j=1}^{N_i} (x_j - \mathrm{mean}(x_j)) (x_j - \mathrm{mean}(x_j))'$$
For my data set this gives not a good estimation of the covariance matrix. What did work was to artifically equalize the amount of samples per class by re-using data, and then estimating the covariance by the first equation. For example for class 1, I artifically increased the number of samples by re-using the 100 samples 50 times to get an amount of 5000 samples for this class, i.e. adding the same data again and again. This seems to remove an apparent bias.
It works quite well, but I have no mathematical explanation or intuition why???
 A: In this kind of analysis you need to be clearer on what covariance matrix you are actually trying to estimate.  (You refer to estimating "the covariance", but the covariance of what with what?)  I am going to assume that "the covariance" you are trying to estimate is the covariance matrix for the vector of deviations of each output vector from the true mean for that category.  This can be framed by writing your model as a standard multivariate regression model, with error vectors representing the deviation of each output from the true mean for the category.  This framework will be used to elucidate the derivation of the pooled covariance estimator, and the conditions under which this is a reasonable estimator.

Where does the pooled covariance estimator come from: Consider a multivariate linear regression model for a response vector $\mathbf{Y}_i$ where we have a single categorical variable $z =1,...,m$ (with $m$ categories) and homoscedastic errors.  We can write the multivariate regression model as:
$$\mathbf{Y}_i = \sum_{k=1}^{m} \boldsymbol{\beta}_k \cdot \mathbb{I}(z_i=k) + \boldsymbol{\varepsilon}_i \quad \quad \quad
\boldsymbol{\varepsilon}_1, ... , \boldsymbol{\varepsilon}_n \sim \text{IID N}(\mathbf{0}, \boldsymbol{\Sigma}).$$
It is important to stress here that this model assumes homoscedastic errors, which means that we assume that the true covariance matrix of the error terms in each of the categories is the same.  We will soon see that this is a crucial assumption for the use of the pooled covariance estimator.
Without loss of generality, suppose that we have $N_1,...,N_m$ values in each of the respective categories, leading to corresponding response means $\bar{\mathbf{y}}_1, ..., \bar{\mathbf{y}}_m$.  It is simple to show that OLS estimation leads to the following coefficient estimators and corresponding residuals:
$$\hat{\boldsymbol{\beta}}_k = \bar{\mathbf{y}}_k \quad \quad \quad \quad \quad
\mathbf{r}_i = \mathbf{y}_i - \bar{\mathbf{y}}_{z_i}.$$
For the next step we will define $\mathbf{y}_{k,i}$ to be the $i$th response vector in the $k$th category.  Using this notation, the unbiased estimator of the covariance matrix of the error vector is:
$$\begin{equation} \begin{aligned}
\hat{\mathbf{\Sigma}} 
&= \frac{1}{N-m} \sum_{i=1}^n \mathbf{r}_i \mathbf{r}_i^\text{T} \\[6pt]
&= \frac{1}{N-m} \sum_{i=1}^n (\mathbf{y}_i - \bar{\mathbf{y}}_{z_i}) (\mathbf{y}_i - \bar{\mathbf{y}}_{z_i})^\text{T} \\[6pt]
&= \frac{1}{N-m} \sum_{k=1}^m \sum_{i=1}^{N_i} (\mathbf{y}_{k,i} - \bar{\mathbf{y}}_k) (\mathbf{y}_{k,i} - \bar{\mathbf{y}}_{z_i})^\text{T} \\[6pt]
&= \frac{\sum_{k=1}^m \sum_{i=1}^{N_i} (\mathbf{y}_{k,i} - \bar{\mathbf{y}}_k) (\mathbf{y}_{k,i} - \bar{\mathbf{y}}_{z_i})^\text{T}}{\sum_{k=1}^m (N_k -1) } \\[6pt]
&= \hat{\mathbf{\Sigma}}_\text{Pooled}. \\[6pt]
\end{aligned} \end{equation}$$
We can see here that, in the case of OLS estimation, the corresponding estimator of the covariance of the error vector is equivalent to the pooled covariance estimator.  That is, this form of estimator arises as a natural unbiased estimator in the case where the output data is treated as coming from a standard multivariate regression based on an observed categorical variable, with homoscedastic error terms.

When is this a good/bad estimator: Since the pooled covariance estimator arises from standard OLS estimation under a model with homoscedastic errors, it should not be surprising that it is generally a good estimator in cases where this model form is appropriate, and a bad estimator in cases where this model form is inappropriate.  In particular, the legitimacy of the pooled covariance estimator hinges on the assumption of homoscedastic error terms, which is equivalent to assuming that the true covariance matrix of the error terms in each of the categories is the same.
If the true covariance matrices for the error terms in difference categories are substantially different then it should not be surprising that the above estimator will be a poor estimator of the assumed common covariance matrix.  Firstly, it means that you are attempting to estimate a single covariance matrix that does not in fact exist (i.e., there is not really a common covariance matrix across categories to begin with), and secondly, you are using OLS estimation that is generally poor in cases of serious heteroskedasticity.
In such cases it is best to revise your model to allow for heteroskedasticity and then use estimation techniques that are appropriate to this.  That might include some parametric specification of a heteroskedastic weighting structure (and corresponding estimation using WLS) or it might mean using basic sample estimates of the covariance across each of the categories, and not attempting to aggregate them.

Issues with your problem: In your problem you have $m=10$ categories, which is not many, but you have noted that they are very unbalanced.  In your description of your problem you way that you are estimating "the covariance" and that you have compared the pooled covariance estimator to "the covariance" and found it performed poorly.  Unfortunately, it is unclear exactly what covariance you are seeking to estimate, and it what covariance matrix you are comparing the pooled estimator to.
To improve your analysis, I would suggest that you begin by making a clearer specification of exactly what covariance you are trying to estimate, and what (if any) assumptions about the variance across groups is reasonable with this data.  I hope the above exposition assists you in understanding the derivation of the pooled covariance estimator and the assumptions behind this.
