Variance of the product of a random matrix and a random vector If $X$ and $Y$ are independent random variables, then the variance
of the product $XY$ is given by
$
\mathbb{V}\left(XY\right)=\left\{ \mathbb{E}\left(X\right)\right\} ^{2}\mathbb{V}\left(Y\right)+\left\{ \mathbb{E}\left(Y\right)\right\} ^{2}\mathbb{V}\left(X\right)+\mathbb{V}\left(X\right)\mathbb{V}\left(Y\right)
$
If $\mathbf{X}$ and $\mathbf{y}$ are independent matrix and vector
of $m\times m$ and $m\times1$ dimension respectively, then what would be the variance
of the product $\mathbf{X}\mathbf{y}$?
My Attempt
$
\mathbb{V}\left(\mathbf{X}\mathbf{y}\right)=\mathbb{E}\left(\mathbf{X}\right)\mathbb{V}\left(\mathbf{y}\right)\left\{ \mathbb{E}\left(\mathbf{X}\right)\right\} ^{\prime}+\left\{ \mathbb{E}\left(\mathbf{y}\right)\otimes\mathbf{I}_{m}\right\} ^{\prime}\mathbb{V}\left\{ \textrm{vec}\left(\mathbf{X}\right)\right\} \left\{ \mathbb{E}\left(\mathbf{y}\right)\otimes\mathbf{I}_{m}\right\} +\mathbb{V}\left\{ \textrm{vec}\left(\mathbf{X}\right)\right\} \left\{ \mathbb{V}\left(\mathbf{y}\right)\otimes\mathbf{I}_{m}\right\} 
$
I know this is not right, at least the last term is wrong. I'd highly appreciate if you give me the right identity or point out any reference. Thanks in advance for your help and time.
 A: The co-variance matrix of $W x$ is
$$
\large{\mathrm{V}[Wx] = \mathrm{diag}(S m) + M \Sigma M^T},
$$
where $\mathrm{diag}(S m)$ is the diagonal matrix with the vector $Sm$ on its diagonal.
This is a generalization of jbowman's answer when we don't assume the entries of the vector are independent, but instead have covariance matrix $\Sigma$.
More precisely:

*

*We let $W$ be a random matrix with iid. entries, mean $\mathrm{E}[W]=M$ and scalar variances $\mathrm{V}[W_{i,j}] = S_{i,j}$.

*We let $x$ be a random vector with co-variance matrix $\mathrm{V}[x]=\Sigma$ and raw second moments $\mathrm{E}[x_i^2] = m_i$.

The the means $W$ work on the covariance matrix in the term $M\Sigma M^T$; and the variances of $W$ work on the second moment of $x$ in the term $\mathrm{diag}(Sm)$.

Proof.
I couldn't find that formula in the Matrix Cookbook, but it follows relatively simply from the law of total variance: $\mathrm{V}(y) = \mathrm{E}(\mathrm{V}(y \mid x)) + \mathrm{V}(\mathrm{E}(y \mid x))$.
For the second term
\begin{align}
    \mathrm{V}(\mathrm{E}(y \mid x))
    &=
    \mathrm{V}(\mathrm{E}(Wx \mid x))
    \\&=
    \mathrm{V}(\mathrm{E}(W) x)
    \\&=
    \mathrm{V}(M x)
    \\&=
    M \Sigma M^T
    .
\end{align}
For the term,
$\mathrm{V}(Wx \mid x)$, note that we may assume $M=\mathrm{E}(W)=0$ since
$\mathrm{V}((W'+M)x \mid x)
= \mathrm{V}(W'x \mid x) + \mathrm{V}(Mx \mid x)$, where the second term is constant (given $x$) and thus has variance 0.
Assuming $x$ is constant and $W$ thus mean 0, $Wx$ is a vector of independent entries, and so the co-variance matrix is just a diagonal.
We can compute those diagonal entries:
\begin{align}
    (\mathrm{V} Wx)_{ii}
    &=
    \mathrm{E}(W_ix)^2 - E[W_ix]^2
    \\&=
    \sum_{k,\ell}
    \mathrm{E}(W_{i,k}W_{i,\ell}) x_k x_\ell - 0
    \\&=
    \label{eq:reduce-sum}
    \sum_k
    \mathrm{E}(W_{i,k}^2) x_k^2
    \\&=
    \sum_k
    S_{i,k} x_k^2
    \\&=
    \langle S_i, x^2\rangle,
\end{align}
where again we used that $W$ is here assumed to be mean 0.
And so $\mathrm{V}(Wx)$ is a diagonal matrix with $S (x^2)$ on the diagonal.
Putting it all together we get
\begin{align}
    \mathrm{V}(y)
    &= \mathrm{E}(\mathrm{V}(y \mid x)) + \mathrm{V}(\mathrm{E}(y \mid x))
    \\&= \mathrm{diag}(S \mathrm E(x^2)) + M\Sigma M^T.
\end{align}
A: I'll assume that the elements of $\mathbf{y}$ are i.i.d. and likewise for the elements of $\mathbf{X}$.  This is important, though, so be forewarned!


*

*The diagonal elements of the covariance matrix equal the sum of $m$ products of i.i.d. random variates, so the variance will equal $m \mathbb{V}(x_{ij}y_j)$, which variance you have above in your first row.

*The off-diagonal elements all equal zero, as the rows of $\mathbf{X}$ are independent.  To see this, without loss of generality assume $\mathbb{E}x_{ij} = \mathbb{E}y_i = 0 \space \forall\thinspace i,j$. Define $\mathbf{x}_i$ as the $i^{\text{th}}$ row of $\mathbf{X}$, transposed to be a column vector.  Then:
$\text{Cov}(\mathbf{x_i^\text{T}y},\mathbf{x_j^\text{T}y}) = \mathbb{E}(\mathbf{x_i^\text{T}y})^\text{T}(\mathbf{x_j^\text{T}y}) = \mathbb{E}\mathbf{y}^{\text{T}}\mathbf{x}_i\mathbf{x}_j^\text{T}\mathbf{y}=\mathbb{E}_y\mathbb{E}_x \mathbf{y}^{\text{T}}\mathbf{x}_i\mathbf{x}_j^\text{T}\mathbf{y}$
Note that $\mathbf{x}_i\mathbf{x}_j^\text{T}$ is a matrix, the $(p,q)^\text{th}$ element of which equals $x_{ip}x_{jq}$.   When $i \ne j$, the expectation with respect to $x$ of $\mathbf{y}^{\text{T}}\mathbf{x}_i\mathbf{x}_j^\text{T}\mathbf{y}$ equals 0 for any $\mathbf{y}$, as each element is just the expectation of the product of two independent r.v.s with mean 0 times $y_py_q$.  Consequently, the entire expectation equals 0.  
A: From [1]
The joint-covariance matrix of the product of a real random matrix $X$ of dimension $m\times m$ and a real random matrix $y$ of dimension $m\times 1$ is a real matrix of dimension $m\times m$. The element   on the $k^\textrm{th}$ row and $l^\textrm{th}$ column of the joint-covariance matrix, denoted as $\operatorname {E} \left[(\mathbf{X} \,\mathbf{y}- \operatorname {E} \left[\mathbf{X} \,\mathbf{y}  \right]   )(\mathbf{X} \,\mathbf{y}- \operatorname {E} \left[\mathbf{X} \,\mathbf{y}  \right]  )^{\top }\right]_{k,l}$, is given as
$$\sum\limits_{i=1}^m\sum\limits_{j=1}^m 
\Bigl(
\operatorname {cov}_X( X_{ki},    X_{lj})
 + \operatorname {E}_X \left[ X_{ki}     \right]
\operatorname {E}_X \left[ X_{lj}    \right]
 \Bigr)\Bigl(
 \operatorname {cov}_Y( y_{i},    y_{j}
 )
 + \operatorname {E}_Y \left[ y_{i}     \right]
\operatorname {E}_Y \left[ y_{j}    \right]
 \Bigr) -\operatorname {E}_X \left[ X_{ki}     \right]
\operatorname {E}_X \left[ X_{lj}    \right] \operatorname {E}_Y \left[ y_{i}     \right]
\operatorname {E}_Y \left[ y_{j}    \right]
$$
Bibliography
[1]
Proof Verification: Joint variance of the product of a random matrix with a random vector
