I have a situation where observed random variables $X_i$ are the sum of two independent (but unobserved) variables, $$X_i = S_i + N_i,$$ (e.g. what you observe is a random signal plus random noise).
I have samples of the $X_i$ and I also know the distribution of $N_i$ (it's multivariate Gaussian with mean 0 and known, but non-diagonal, covariance). How do I estimate the covariance among the $S_i$?
My naive attempt:
Since $S_i$ and $N_i$ are independent, $\mathrm{Cov}(X) = \mathrm{Cov}(S) + \mathrm{Cov}(N)$. I can take my samples of $X_i$, calculate the sample covariance, and subtract off the known covariance matrix of $N$, yielding an estimate of the covariance of $S$. The problem is that with a finite sample size, the sample covariance of $X$ is not a perfect measurement of $\mathrm{Cov}(X)$, and when I do the subtraction I end up with an estimate of $\mathrm{Cov}(S)$ that is not positive-definite.
This seems like it's probably a straightforward statistics problem that I don't know the terminology to search for. Is there a simple, reliable way to "subtract off" a known component from a measured sample covariance that yields a valid (i.e. positive definite) covariance estimate?
(I would be happy to assume $S_i$ is multivariate Gaussian if it helps. Unfortunately, the dimension of the variables is about $1500$ so it seems like a tall order to parameterize the covariance matrix of $S$ with $1500^2 \approx 2$ million free parameters and try to fit them simultaneously.)