Statistically you would apply Bessel's correction to address the bias you point out. However, that misses the point that the variance-covariance matrix is non-stationary, suffers from the curse of dimensionality, and that the noisy mean return estimates have significantly more impact than a biased covariance matrix on portfolio weights.
The best ways to build a covariance matrix are to:
- build a multi-factor risk model (see BARRA, Axioma, Northfield, or Finanalytica literature for examples). There is quite a bit of nuance in building these properly (alpha and risk interactions, adjustments for out-of-sample bias, factor identification, errors-in-variables bias, etc.)
- cleanse the sample covariance matrix using random matrix theory
- estimate the matrix via some type of exponential weighting
- use a shrinkage estimate such as Ledoit and Wolf which blends the sample covariance with a prior (usually a constant covariance, constant correlation, or identity matrix)
- Or some combination of the above
I would recommend taking a shrinkage approach as in Ledoit's aptly titled "Honey I Shrunk the Covariance Matrix" with a constant covariance prior since it is quite easy to implement and generates good results. There are over 200+ citations from the original approach which cover extensions such as how much weight to assign to the prior and the sample covariance matrix.
Note that estimation error in expected returns has been estimated to be about 10x more important than the estimation error in variances and 20x times as important as estimation in covariances (see Ziemba 2003). So you may want to do a decent job on the risk side of the utility function, and then concentrate your efforts on addressing the noisiness in expected returns through robust Bayesian optimization procedures.
