Variance of product of 2 independent random vector This should be a fairly straight forward but I still couldn't quite get it.
Let $X$ be $n$ by $1$ random vector, $Y$ be another $n$ by $1$ and that $X$ and $Y$ are independent.
What is $Var(X^tY)=:V_{XY}$?
I tried generalizing this formula from the univariate case
$Var(xy)=E(x)^2Var(y)+E(y)^2Var(x)+Var(x)Var(y)$
except I hit a bit of a roadblock for the 3rd term above. That is, I get
$$
V_{XY}=X^tV_YX+Y^tV_XY\\+E[(X^t-\mu_X^t)(Y-\mu_Y)(Y^t-\mu_Y^t)(X-\mu_X)]
$$
So... what is it?
 A: This answer supposes that $X^TY$ (where $X$ and $Y$ are $n\times 1$ vectors) is a $1\times 1$ vector or scalar $\sum_i X_iY_i$ and so we need to consider the
variance of a single random variable that is this sum of products.  Since $X$
and $Y$ are independent random vectors, we note that $X_1, Y_1$ are independent
random variables as are $X_2, Y_2$. Also, $(X_1, X_2)$ is independent of 
$(Y_1, Y_2)$.
Since $\operatorname{var}(Z_1+Z_2) = \operatorname{var}(Z_1) + \operatorname{var}(Z_2) + 2\operatorname{cov}(Z_1,Z_2)$, we get that
$$\begin{align}
\operatorname{var}(X_1Y_1+X_2Y_2) &= \operatorname{var}(X_1Y_1)+
\operatorname{var}(X_2Y_2) + 2\operatorname{cov}(X_1Y_1, X_2Y_2)\\
\end{align}$$
The first two terms are the variances of the products of independent
random variables, and the OP knows the formula for handling these.
Turning to the covariance, we have that
$$\begin{align}
\operatorname{cov}(X_1Y_1, X_2Y_2) &=
E[X_1Y_1X_2Y_2] - E[X_1Y_1]E[X_2Y_2]\\
&= E[X_1X_2]E[Y_1Y_2] - E[X_1]E[Y_1]E[X_2]E[Y_2]\\
&= \left(\operatorname{cov}(X_1, X_2)+E[X_1]E[X_2]\right)\left(\operatorname{cov}(Y_1, Y_2)+E[Y_1]E[Y_2]\right)\\
& \qquad \qquad - E[X_1]E[Y_1]E[X_2]E[Y_2]\\
\end{align}$$
which simplifies to
$$
\operatorname{cov}(X_1Y_1, X_2Y_2) = 
\operatorname{cov}(X_1, X_2)E[Y_1]E[Y_2]+ \operatorname{cov}(Y_1, Y_2)E[X_1]E[X_2]\\ + \operatorname{cov}(X_1, X_2)\operatorname{cov}(Y_1, Y_2),\tag{1}$$
a result that is eerily similar to the result
$$\operatorname{var}(X_iY_i) = \operatorname{var}(X_i)(E[Y_i])^2 + \operatorname{var}(Y_i)(E[X_i])^2 + \operatorname{var}(X_i)\operatorname{var}(Y_i) \tag{2}$$
quoted by the OP for independent random variables $X$ and $Y$
Thus we have
$$\operatorname{var}\left(\sum_{i=1}^n X_iY_i\right)
= \sum_{i=1}^n \operatorname{var}(X_iY_i) + 2\sum_{i=1}^{n-1}\sum_{j=i+1}^n
\operatorname{cov}(X_iY_i, X_jY_j)$$
where each term in the first sum on the right is given by $(2)$ and
each term in the double sum on the right is given by $(1)$. Plug and chug and enjoy!
A: Here's a matrix-notation version that might be more convenient to work with than the plug-and-chug required by Dilip's answer.
$\DeclareMathOperator{\E}{\mathbb E} \DeclareMathOperator{\Var}{Var} \newcommand{\tp}{^\mathsf{T}} \DeclareMathOperator{\Tr}{Tr}$
Let $\mu_X := \E[X]$, $\Sigma_X := \Var[X]$, and similarly for $\mu_Y$, $\Sigma_Y$.
Then,
using trace rotation and independence of $X$ and $Y$, we get
\begin{align}
\Var[X\tp Y]
  &= \E[(X\tp Y)^2] - \E[X\tp Y]^2
\\&= \E\left[ X\tp Y Y\tp X\right] - \left(\E[X]\tp \E[Y]\right)^2
\\&= \E\left[\Tr\left( X X\tp Y Y\tp \right)\right] - \left(\mu_X\tp \mu_Y\right)^2
\\&= \Tr\left( \E[X X\tp] \E[Y Y\tp] \right) - \left(\mu_X\tp \mu_Y\right)^2
\\&= \Tr\left(
        \left( \mu_X \mu_X\tp + \Sigma_X \right)
        \left( \mu_Y \mu_Y\tp + \Sigma_Y \right)
     \right) - \left(\mu_X\tp \mu_Y\right)^2
\\&= \Tr\left( \mu_X \mu_X\tp \mu_Y \mu_Y\tp \right)
   + \Tr\left( \mu_X \mu_X\tp \Sigma_Y \right)
   + \Tr\left( \Sigma_X \mu_Y \mu_Y\tp \right)
   + \Tr\left( \Sigma_X \Sigma_Y \right)
   - \left(\mu_X\tp \mu_Y\right)^2
\\&= \left(\mu_X\tp \mu_Y\right)^2
   + \Tr\left( \mu_X\tp \Sigma_Y \mu_X \right)
   + \Tr\left( \mu_Y\tp \Sigma_X \mu_Y \right)
   + \Tr\left( \Sigma_X \Sigma_Y \right)
   - \left(\mu_X\tp \mu_Y\right)^2
\\&= \mu_X\tp \Sigma_Y \mu_X + \mu_Y\tp \Sigma_X \mu_Y + \Tr(\Sigma_X \Sigma_Y)
.\end{align}
