$\newcommand{\cov}{\operatorname{cov}}\newcommand{\corr}{\operatorname{corr}}$I have two dependent multivariate random variables, or random vectors, $X$ and $Y$, their respective covariance matrices $\Sigma_X = \cov(X,X)$, $\Sigma_Y = \cov(Y,Y)$ and a non-zero correlation matrix $P_{XY} = \corr(X, Y)$.
I would like to compute the covariance matrix of the sum of the two random vectors:
\begin{align}
Z {}={}& X + Y
\\
\Sigma_Z {}={}& \Sigma_X + \Sigma_Y + 2 \cov(X, Y)
\end{align}
Intuitively, I would assume the cross-covariance can be computed as follows:
$$\cov(X, Y) = P_{XY} \sqrt{\Sigma_X} \sqrt{\Sigma_Y}$$
But I can see that this matrix is not necessarily symmetric, and therefore the resulting $\Sigma_Z$ wouldn't be either. Am I wrong? Is it at all possible to compute $\Sigma_Z$ from the information I have? If not, is there some way to approximate it?
EDIT: I think it isn't correct to write:
$$ \Sigma_Z = \Sigma_X + \Sigma_Y + 2 \cov(X, Y) $$
as $\cov(X,Y) \ne \cov(Y,X)$, but instead $\cov(X,Y) = \cov(Y,X)^T$.
Rewriting the equation for $\Sigma_Z$ to:
\begin{align}
\Sigma_Z {}={}& \Sigma_X + \Sigma_Y + \cov(X, Y) + \cov(Y, X)
\\
{}={}& \Sigma_X + \Sigma_Y + \cov(X, Y) + \cov(X, Y)^T
\end{align}
I can now prove that indeed $\Sigma_Z$ is guaranteed to be symmetric. As antkam points out in the chat below, it isn't necessary that $\cov(X, Y) + \cov(X, Y)^T$ be positive semi-definite for $\Sigma_Z$ to be positive semi-definite.
Best Answer
Your original intuition does not work, i.e., $cov(X, Y) \neq P_{XY} \sqrt{\Sigma_X} \sqrt{\Sigma_Y}$, because the definition of $corr(X_i, Y_j)$ involves scalar divisions, which cannot be reversed like this. For one thing, you have three matrices, so how do you pick the order of multiplication?
Instead, I would go back to the per-entry definitions:
$$[cov(X,Y)]_{ij} = cov(X_i, Y_j) = corr(X_i, Y_j) \ \sigma(X_i) \ \sigma(Y_j) = [P_{XY}]_{ij} \sqrt{ [\Sigma_X]_{ii}} \sqrt{[\Sigma_Y]_{jj}}$$
As you can see, the non-diagonal entries of $\Sigma_X, \Sigma_Y$ don't enter the calculations at all. So I'd think you cannot write this in terms of standard matrix multiplications.
However, if you have a $Diag$ function that keeps the diagonal entries but zeroes out all the non-diagonal entries, then you can write:
$$Cov(X,Y) = \sqrt{Diag(\Sigma_X)} \ P_{XY} \ \sqrt{Diag(\Sigma_Y)}$$
The multiplication is in that order because left-multiplying by a diagonal matrix scales the rows and right-multiplying scales the columns. (I know some typical software packages have a similar function although in the case of python np.diag you actually have to call it twice, so be careful!)