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.



Variance of projection of vectors onto random subspace
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).$$


