$$ Y_i = \sum_j B_{ij} X_j$$
$$ covar(X_i, X_j) = V_X = \delta_{ij} var(X_i)^2 $$
$$ covar(B_{ij},B_{kl}) \neq 0 $$
$$ V_Y = ? $$
I know, that if $B$ was fixed, it is straight forward, but I would like to generalize it. I managed only by inventing some notation myself, which a physicist should not do I guess. I will check my method later with Monte Carlo, but maybe there is a formal way to write it down?
Maybe the answer should be generalized from variance-of-product-of-2-independent-random-vector?
Edit
We can also assume that $covar(X_i,B_{kl})=0$.
Edit2
Replaced Einstein Summation rule with a sum. I don't just skip the indexes as I want to to stress the dimensions and make the latter more readable.