This question is similar to this question, but is generalized for rectangular matrices.
Suppose we have matrix $\bf{S} \in \mathbb{R}^{m \times n} $ and vector $\bf{y} \in \mathbb{R}^{n}$.
$\textbf{y}$ is sampled from a multivariate normal distribution with $\Sigma = \sigma^2 I$, where $\mu$ and $\sigma$ is known:
$$ \textbf{y} \sim \mathcal{N}(\mu, \sigma^2 I) $$
$\textbf{S}$ is a fixed, known matrix
What is the variance of the product of $\textbf{Sy}$?
What I have tried so far (attempted working off of same question): \begin{align*} \mathrm{Var}(\textbf{Sy}) &= \mathbb{E}[(\textbf{Sy})^2]-\mathbb{E}[(\textbf{Sy})]^2\\ &= \mathbb{E}[(\textbf{Sy})^T\textbf{Sy}]-\mathbb{E}[\textbf{Sy}]^T\mathbb{E}[\textbf{Sy}]\\ &= \mathbb{E}[\textbf{y}^T\textbf{S}^T\textbf{Sy}] - \mathbb{E}[\textbf{y}^T\textbf{S}^T]\mathbb{E}[\textbf{Sy}] \end{align*}
From here, it seems I would like to remove $\textbf{S}$ from the expected value, but I am not sure how to proceed.