$\DeclareMathOperator{\Var}{Var}$
$\DeclareMathOperator{\Cov}{Cov}$
$\DeclareMathOperator{\tr}{tr}$
The variance of $z$ can be computed as follows (suppose $A$ is symmetric, which is standard when the quadratic form is discussed. For the non-symmetric case, rewrite $z$ as $z = x'Bx + b'x$, where ...
I get a similar expression, but off by some of the factors, and also the second-to last term I have is different.. I will do some numerical verifications to check
@becko As I wrote at the beginning of this answer, if $A$ is non-symmetric, then you need to replace $A$ in the last expression with $(A + A')/2$ -- that may be what you meant by "off some of the factors".
This is the expression I get for the variance: $\frac{1}{2}tr(\mathbb{A\Sigma A\Sigma}) + \mathbf{u}^{\top}\mathbb{A\Sigma A}\mathbf{u} + 2\mathbf{b}^{\top}\mathbb{\Sigma A}\mathbf{u} + \mathbf{b}^{\top}\Sigma\mathbf{b}$, assuming $A$ is symmetric.
I still don't agree with the variance though ... if I take my expression in the previous comment, and multiply all the $A$ by 2, you see I never get the $6\mu A\mu b\mu$ term you get.
I don't see how. I mean, it would just replace A with 2A.
I think we can at least agree on the following expression for the variance, if the x are standardized (zero mean, unit variances): $$\Tr(\mathbb{A}\mathbb{A}) + \mathbf{b}^{\top}\mathbf{b}$$
Discussion between Zhanxiong and becko
Discussion between Zhanxiong and becko