Estimating the covariance of $X$ using i.i.d. samples of $X + Y$, where $Y$ has known statistics Let $X$ be a vector-valued random variable whose mean and covariance we would like to estimate. However, we only have access to i.i.d. samples of $Z = X + Y$, where $Y \perp\!\!\!\perp X$ is some other vector-valued random variable with zero mean and (known) covariance $\Sigma_Y$.
The first thing that comes to mind is to use
$$ \bar{\Sigma}_X = \frac{1}{N - 1} \sum^{N}_{j = 1} \left(Z_j - \bar{\mu}\right) \left(Z_j - \bar{\mu}\right)^T - \Sigma_Y \;\; \text{ where } \;\;\bar{\mu} = \frac{1}{N} \sum^{N}_{j = 1} Z_j \;.$$
However, this estimate is not always positive semi-definite. Therefore, I believe we should be able to better than this. Any ideas, or pointers to relevant results, are appreciated. Thanks!
 A: Deleted the original answer, since it was wrong. Here is a new suggestion:
If $\Sigma_X$ acquired by $\hat{\Sigma}_Z - \Sigma_Y$ is indeed not positive semidefinite, it indicates that the estimator $\hat{\Sigma}_Z$ is not a good estimator for $\Sigma_Z$ probably due to small sample size. Since by design, $\Sigma_Z = \Sigma_X + \Sigma_Y$, should be positive definite. However, since $\hat{\Sigma}_Z$ still may have some information in the covariance structure, what if we take a Bayesian way? We set a prior for $\Sigma_X$ to use $Z=X+Y$ information and use observed $Z$ for the likelihood function.
\begin{align}
p(\Sigma_X) =\text{Inverse-Wishart}(\Psi, \nu)\\
p(Z|\Sigma_X) = \mathcal{N}(\mu_Y, \Sigma_X + \Sigma_Y)
\end{align}
Then, we will get posterior distribution for $p(\Sigma_X|Z)$ and the point estimate of the distribution will always be positive semi-definite.
A: It's easy to see how this could happen. Take the 1-D case, and as an extreme. consider the case $Var(X) = 0$. Then if sample variance of $Y$ is less than assumed true variance of $Y$, your variance estimator of $X$ will be negative (i.e., not positive semi-definite) even if all your assumptions are correct.
Here is an approach which is an alternative to the answer provided by @villybyun . I offer no opinion on its relative merits.
Adjust {sample covariance of $Z$ - covariance of $Y$} by adding the minimal Frobenius norm matrix, such that the adjusted matrix has minimum eigenvalue $\ge$ specified eigenvalue, call it mineig. The choden value of mineig could be $0$ or some positive number. 
This can readily be formulated and solved (providing the matrix is not excessively large) as a convex SemiDefinite program (SDP). Here is how it could be coded in CVX http://cvxr.com , which is quite popular in the Stanford EE dept.
cvx_begin
variable A(n,n) symmetric
minimize(norm(A,'fro'))
subject to
CovZsample - CovY + A -mineig * eye(n) == semidefinite(n)
cvx_end

The constraint can be expressed more succinctly using CVX's lambda_min function (minimum eigenvalue) as
lambda_min(CovZsample - CovY + A) >= mineig

Note that declaring A to be symmetric in the variable declaration, as shown above, is not necessary, because symmetry (of the "argument") will be enforced by either version of the constraint.
Your adjusted covariance estimate for $X$ is CovZsample - CovY + A, using the optimal; value of $A$ found by solving the SDP.
