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.com Sep 18 '16 at 5:26

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

  • $\begingroup$ 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. $\endgroup$ – iheap Sep 17 '16 at 23:29
  • $\begingroup$ I guess you can simply delete this one and create it in stats.stackexchage.com $\endgroup$ – leonbloy Sep 18 '16 at 1:34
  • $\begingroup$ @leonbloy You can vote to close the question as off-topic (This question belongs on another site in the Stack Exchange network). $\endgroup$ – 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.

  • $\begingroup$ Possible error in formula, shouldn't it be a minus before $\Sigma_Y$ in the last formula? $\endgroup$ – tomka Sep 18 '16 at 10:40
  • $\begingroup$ @tomka Changed the typo in the first formula (sum of variances). I think now it leads to the last formula (unchanged) correctly. $\endgroup$ – villybyun Sep 18 '16 at 12:09
  • $\begingroup$ @Greenparker, right. I am deleting my answer! $\endgroup$ – villybyun Sep 18 '16 at 12:57
  • $\begingroup$ Deleted the content of the first answer due to wrong logic and re-posted another suggestion. $\endgroup$ – villybyun Sep 18 '16 at 13:49
  • $\begingroup$ 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? $\endgroup$ – 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.

variable A(n,n) symmetric
subject to
CovZsample - CovY + A -mineig * eye(n) == semidefinite(n)

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.

  • $\begingroup$ 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. $\endgroup$ – iheap Sep 18 '16 at 19:55
  • $\begingroup$ 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. $\endgroup$ – Mark L. Stone Sep 18 '16 at 19:59
  • $\begingroup$ Hi Mark. It seems like what I proposed coincides with the actual optimal correction. Check out this paper by Higham. $\endgroup$ – iheap Sep 19 '16 at 18:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.