# 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!

## migrated from math.stackexchange.comSep 18 '16 at 5:26

This question came from our site for people studying math at any level and professionals in related fields.

• Thanks @leonbloy. How can I get this question moved to the stats site? I assume opening duplicate questions is not the right approach, but I'm not clear on how to proceed. – iheap Sep 17 '16 at 23:29
• I guess you can simply delete this one and create it in stats.stackexchage.com – leonbloy Sep 18 '16 at 1:34
• @leonbloy You can vote to close the question as off-topic (This question belongs on another site in the Stack Exchange network). – Math1000 Sep 18 '16 at 3:49

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.

• Possible error in formula, shouldn't it be a minus before $\Sigma_Y$ in the last formula? – tomka Sep 18 '16 at 10:40
• @tomka Changed the typo in the first formula (sum of variances). I think now it leads to the last formula (unchanged) correctly. – villybyun Sep 18 '16 at 12:09
• @Greenparker, right. I am deleting my answer! – villybyun Sep 18 '16 at 12:57
• Deleted the content of the first answer due to wrong logic and re-posted another suggestion. – villybyun Sep 18 '16 at 13:49
• Thanks! I attempted to derive a simple closed-form expression for $p(\Sigma_X|Z)$, but failed. The main difficulty I had was dealing with $(\Sigma_X + \Sigma_Y)^{-1}$. Do you guys think it is possible to derive a simple closed-form expression for the posterior, or do I have resort to numerics? – iheap Sep 18 '16 at 18:53

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.

• Thanks Mark. We can also pick $A = - \Sigma_j \min \left(0, \lambda_j \right) v^{\phantom{T}}_j \! v^T_j$, where $v_j$ is the $j$-th eigenvector of $\bar{\Sigma}_X$. I wonder how this will compare with picking $A$ according to your criteria. – iheap Sep 18 '16 at 19:55
• Try it and see. Or alternatively, add a multiple of the identity matrix with diagonal elements equal to (-smallest eigenvalue + mineig). My method will have a smaller adjustment (generally) to the diagonal, due to making adjustments elsewhere as well. – Mark L. Stone Sep 18 '16 at 19:59
• Hi Mark. It seems like what I proposed coincides with the actual optimal correction. Check out this paper by Higham. – iheap Sep 19 '16 at 18:30