Is it possible to regularize a covariance matrix? I have many "parallel" time-series (about 100). They have relatively short history. I calculate a covariance matrix between these time-series.
Now, I believe that the observed covariance matrix is very inaccurate / noisy. For example some time-series might look to be very correlated just because of noise / chance / coincidence.
My question is if there is a method for smoothing / regularization (in machine learning sense) of the covariance matrix.
 A: You are trying to model high dimensional data with few observations.  As @whuber notes, you can add a diagonal matrix to your sample covariance matrix.  This and similar approaches reduce estimated correlation across the series, although the observed data will eventually dominate estimation process (however, you claim few observations).
Many covariance estimation frameworks, including Bayesian and maximum entropy frameworks, focus on the precision matrix (the inverse covariance matrix). In Bayesian analysis, the precision matrix is often characterized by inverse Wishart distribution conjugate priors. For Gaussian random vectors, zero entries in the precision matrix imply conditional independence of corresponding pair of variables (see graphical model literature).
In a Bayesian framework, look at Alvarez, Ignacio, Jarad Niemi, and Matt Simpson. "Bayesian inference for a covariance matrix." arXiv (2014).
Alvarez et al provide experiments with different methods, and report that  in certain circumstances a prior analogous to what @whuber suggests may significantly bias the results towards low correlation.  I don't think this is a problem, because I too prefer the data to demonstrate the correlation before I entertain complex models. Over-regularizing (or under-data-collecting) is expected to be problematic when model building.
In a maximum entropy framework, the identity matrix is the maximum entropy starting point for the correlation matrix and for the inverse correlation matrix. Pairwise or multiple variable composite constraints can be imposed.  For pairwise constraints on covariance, see Dempster, Arthur P. "Covariance selection." Biometrics (1972).
Be sure to read past the summary. The summary is misleading, "setting elements of the inverse of the covariance matrix to zero" is not shown in Dempster, rather non-zeros are iteratively added to the precision matrix in section 3. The demonstrated algorithm is an iterative maximum entropy approach, where pairwise constraints are iteratively added until the model is adequate. This is explained further in the appendix of my article Kevin R. Keane, Jason J. Corso "The Wrong Tool for Inference - A Critical View of Gaussian Graphical Models." ICPRAM (2018).
I impose sets of multivariable constraints in Keane, Kevin R. "Portfolio Variance Constraints." ICPRAM (2019).  The covariance model for stock log-price-return time series is obtained by imposing constraints requiring that exchange traded funds' ("ETFs", linear combinations of variables) observed variance match model variance.  For your scenario, an $N \times N$ covariance model may be constructed by identifying the maximum entropy covariance matrix subject to $N$ individual variance constraints and a few $J$ ETF variance constraints.  The estimation burden is thus only $N+J$ parameters rather than $\frac{N(N+1)}{2}$ parameters.
