Variance of projection of vectors onto random subspace Let $x_i, y_i \in \mathbb{R}^n$ for $i=1, \dots, k < n$ satisfy
$$
\sum_{i=1}^k x_i^\top y_i = 0.
$$
Let $E$ be a random subspace of dimension $m < n$ in $\mathbb{R}^n$ distributed uniformly on the Grassmannian $G_m(\mathbb{R}^n)$ and let $P_E$ be orthogonal projection onto $E$. What is the value of
$$
\mathbb{E}\left[ \left(\sum_{i=1}^k x_i^\top P_E y_i \right)^2 \right] = \dots ?
$$
Any references or suggestions would also be very helpful.
It's quite easy to get an upper bound using Johnson-Lindenstrauss style arguments but I feel that this should be possible to compute exactly. Also (by orthogonal invariance) the first moment vanishes.
 A: It suffices to calculate $\mathbb E ( (P_E)_{ab} (P_E)_{cd})$ for all $1\leq a,b,c,d\leq n$, or, equivalently, calculate $\mathbb E( P_E \otimes P_E ) \in \mathbb R^n \otimes \mathbb R^n \otimes \mathbb R^n \otimes \mathbb R^n$.
This is an $O(n)$-invariant tensor. The space of $O(n)$-invariant tensors in $\mathbb R^n \otimes \mathbb R^n \otimes \mathbb R^n \otimes \mathbb R^n$ is three-dimensional, generated by the three "diagonal" tensors. Thus
$$\mathbb E ( (P_E)_{ab} (P_E)_{cd}) = C_1 \delta_{a=b}\delta_{c=d}+ C_2 \delta_{a=c} \delta_{b=d} + C_3 \delta_{a=c} \delta_{b=d}$$ for some $C_1, C_2, C_3 \in \mathbb R$.
We have $$ m^2 = \sum_{i=1}^n (P_E)_{ii} \sum_{j=1}^n (P_E)_{jj} = \sum_{i=1}^n \sum_{j=1}^n \mathbb E (  (P_E)_{ii}  (P_E)_{jj}) = n^2 C_1 + n C_2 + n C_3 $$
We have $$ m = \sum_{i=1}^n \sum_{j=1}^n (P_E)_{ij}^2 =  \sum_{i=1}^n \sum_{j=1}^n \mathbb E (  (P_E)_{ij}  (P_E)_{ji}) =  n C_1 + n C_2 + n ^2 C_3 $$
and similarly
$$ m = \sum_{i=1}^n \sum_{j=1}^n (P_E)_{ij}^2 =  \sum_{i=1}^n \sum_{j=1}^n \mathbb E (  (P_E)_{ij}  (P_E)_{ij}) =  n C_1 + n^2 C_2 + n C_3 $$
This gives $$C_1 = \frac{(m^2-m) (n+2) + m (n-m)}{n (n-1) (n+2) } , C_2 = \frac{m (n-m)}{ n (n-1)(n+2)}, C_3 = \frac{ m (n-m)}{ n (n-1) (n+2) } $$
We have
$$\mathbb E \left[ \left( \sum_{i=1}^k x_i^T P_E y_i\right)^2 \right] = \mathbb E\left[ \sum_{i=1}^k \sum_{j=1}^k x_i^T P_E y_i x_j^T P_E y_j \right]$$ $$ = C_1 \sum_{i=1}^k \sum_{j=1}^k (x_i \cdot y_i) (x_j \cdot y_j) + C_2 \sum_{i=1}^k \sum_{j=1}^k (x_i \cdot x_j) (y_i \cdot y_j) + C_3 \sum_{i=1}^k \sum_{j=1}^k (x_i \cdot y_j) (y_i \cdot x_j) $$
Your first term vanishes by assumption, so we get
$$ \frac{m (n-m)}{ n (n-1) (n+2) } \left( \sum_{i=1}^k \sum_{j=1}^k (x_i \cdot x_j) (y_i \cdot y_j) + C_3 \sum_{i=1}^k \sum_{j=1}^k (x_i \cdot y_j) (y_i \cdot x_j) \right).$$
Expressed in terms of the matrix $M = \sum_{i=1}^k x_i y_i^T$, this is
$$ \frac{m (n-m)}{ n (n-1) (n+2) }\operatorname{tr}( M M^T + M^2).$$
