Distribución nula de similitud de subespacios, o lo que es lo mismo, la distribución de $\mathrm{tr}(AA'BB')$ ?

Question

¿Cuál es la distribución de $\mathrm{tr}(AA'BB')$ donde $A$ y $B$ son dos matrices aleatorias de $d \times k$ tamaño con columnas ortonormales?

¿Quizás el valor esperado sea más fácil de calcular? Una solución alternativa sería utilizar una simulación. ¿Cuál sería el esquema más eficaz? Valores típicos para $d$ sería alrededor de 2000, mientras que $k$ oscila entre ~10 y algunos cientos.

---

A continuación, se explica con más detalle mi problema y su contexto, cómo acabé planteando esta pregunta y qué intenté.

Contexto

Quiero comprobar si las componentes principales calculadas a partir de una muestra de proceso estocástico han convergido. Mis ideas actuales consisten en comparar los subespacios abarcados por el primer $k$ componentes principales para determinados valores de interés de $k$ ya sea para varias realizaciones del proceso estocástico o para componentes principales bootstraped. Mi criterio para la similitud del subespacio es $\mathrm{tr}(AA'BB') / k$ donde $A$ y $B$ son matrices cuyo $k$ Las columnas son las bases de los dos subespacios a comparar. Este criterio es fácil de calcular y se comporta bien, salvo por la siguiente propiedad: a medida que la dimensión de los subespacios se acerca a la dimensión total, el espacio angular restante en el que señalar las direcciones "erróneas" se reduce. Para construir un criterio más significativo, pensé en comparar esta puntuación con la obtenida al comparar dos subespacios aleatorios de dimensión $k$ .

Mi intento

Mi primer intento fue considerar que, sin pérdida de generalidad, el primer subespacio aleatorio podría tener como base el primer $k$ vectores de la base canónica.

Se puede construir una base para el otro subespacio aleatorio eligiendo primero vectores de la base canónica sin reemplazo.

La distribución resultante sería simplemente la de una ley hipergeométrica con parámetro correspondiente a $k$ se extrae en un conjunto total de $d$ vectores entre los cuales $k$ dar un resultado positivo (el $k$ primeros vectores de la base canónica), donde $d$ es la dimensión del espacio total, $\mathrm{H}(d, k, k/d)$ .

Ahora bien, no hay ninguna razón para que los vectores de las dos bases estén alineados u ortogonales. Supongo que es posible remediarlo aplicando una rotación aleatoria $R$ y mira $\mathrm{tr}(AA'RBB'R')$ . No estoy seguro de cómo una rotación en $\mathbb{R}^d$ se comporta pero tal vez utilizando las propiedades de la traza y el hecho de que $R' = R^{-1}$ ¿es posible solucionar esto?

Nota: Los proyectores ortogonales aleatorios se distribuyen según la Distribución de Wishart . Sin embargo, no sé más sobre esto.

---

Referencias relacionadas:

• Absil, Edelman & Koev, " Sobre el mayor ángulo principal entre subespacios aleatorios ", Linear Algebra and its Applications, 2006, doi:10.1016/j.laa.2005.10.004, (No he leído este)

• Björck & Golub, " Métodos numéricos para calcular los ángulos entre subespacios lineales ", Mathematics of Computation, 1973

• Ipsen & Meyer, " El ángulo entre subespacios complementarios ", The American Mathematical Monthly, 1995

• Johnstone, " Análisis multivariante y conjuntos de Jacobi: Mayor valor propio, límites de Tracy-Widom y tasas de convergencia ", Ann. Statist., 2008 (demasiado esotérico para mí)

• Liquet & Saracco, " Aplicación del enfoque Bootstrap a la elección de la dimensión y la $\alpha$ Parámetro en el {SIR} $\alpha$ Método ", Communications in Statistics Simulation and Computation, 2008

• Wang, Wang & Feng, " Análisis de la distancia subespacial con aplicación al algoritmo bayesiano adaptativo para el reconocimiento facial ", Pattern Recognition , 2006

• Zuccon, Azzopardi y van Rijsbergen, " Espacios semánticos: Medición de la distancia entre distintos subespacios ", 2009

• Hotelling, " Relaciones entre dos conjuntos de variables ", Biometrika, 1936 (En los párrafos 11 y siguientes se dan pruebas posiblemente apropiadas en casos de baja dimensión y se extienden en el párrafo 15 a dimensiones más altas).

• Bao, Hu, Pan y Zhou, " Prueba de independencia para vectores aleatorios de alta dimensión basada en la libertad de las matrices de correlación de bloques ", Revista Electrónica de Estadística, 2017

Answer 1

Alguna vez he hecho una pregunta que es esencialmente un caso especial de la suya cuando $k=1$ : Distribución de un producto escalar de dos vectores unitarios aleatorios en $\mathbb{R}^D$ . @whuber dio una excelente respuesta, donde explica que un producto punto es igual a $t=2u-1$ , donde $$u\sim \mathrm{Beta}((d-1)/2,(d-1)/2)$$ y $d$ es la dimensionalidad del espacio. De ello se desprende que $\mathrm{Var}[t]=1/d$ (también se puede mostrar esto directamente, ver la respuesta de @Student001 en el hilo enlazado).

Su pregunta es sobre una variable aleatoria $$w=\mathrm{tr}(AA^\top BB^\top) = \|A^\top B\|^2,$$ donde $A$ y $B$ son $d\times k$ matrices con columnas ortonormales. Nótese que $A^\top B$ es un $k \times k$ matriz cuadrada, donde cada elemento es un producto punto entre dos vectores unitarios aleatorios en $d$ dimensiones. Los distintos vectores no son independientes (porque tienen que ser ortogonales), pero con $d=2000 \gg k$ Espero que esto pueda ser ignorado. Entonces podemos considerar $A^\top B$ como una matriz de $k^2$ sorteos independientes de $t$ y su norma al cuadrado es una variable aleatoria $$w=\sum_{i=1}^{k^2} t^2_i = \sum_{i=1}^{k^2} (2u_i-1)^2.$$

Creo que sería difícil obtener una expresión analítica para la PDF de una suma de $k^2$ variables aleatorias con distribución beta, pero gracias al teorema del límite central, se acercará rápidamente a una distribución normal (véase un hilo relacionado en math.SE ). Para especificar esta distribución normal, necesitamos calcular su media y su varianza. La media es fácil: media de $t$ es cero, por lo que la media de $t^2$ es igual a la varianza de $t$ que es $1/d$ . Esto significa que $\mathbb E[w]=k^2/d$ .

El cálculo de la varianza es un gran lío que empecé pero no puedo terminar que vergonzosamente me llevó horas. Aquí hay algunas fórmulas auxiliares que derivé buscando la fórmula de los momentos brutos de la distribución beta $$\mathbb E[X^q] = \prod_{r=0}^{q-1}\frac{\alpha+r}{\alpha+\beta+r}, \;\;\; X \sim \mathrm{Beta}(\alpha,\beta),$$ y conectando $\alpha=\beta=(d-1)/2$ : \begin{align} \mathbb E[u]&=1/2, \\ \mathbb E[u^2] &= \frac{d+1}{4d}, \\ \mathbb E[u^3] &= \frac{d+3}{8d}, \\ \mathbb E[u^4] &= \frac{(d+3)(d+5)}{16d(d+2)}. \end{align}

A partir de esto, se puede derivar la varianza de $w$ a partir de $$\mathrm{Var}[w] = k^2 \mathrm{Var}[(2u-1)^2] = k^2 \mathbb E[(2u-1)^4]-k^2(\mathbb E[(2u-1)^2])^2.$$ Omitiré la tediosa aritmética y pasaré directamente a la respuesta: $$\mathrm{Var}[w] = k^2 \frac{2(d-1)}{d^2(d+2)}.$$ La conclusión es que asintóticamente $$w \mathrel{\dot\sim} \mathcal N\left(k^2\frac{1}{d}, \; k^2 \frac{2(d-1)}{d^2(d+2)}\right) \mathrel{\dot\sim} \mathcal N\left(\frac{k^2}{d}, \; \frac{2k^2}{d^2}\right).$$

Una rápida simulación en Matlab confirma este resultado:

Este es el código que he utilizado para producir esta figura ( $d=2000$ , $k=50$ El número de repeticiones de Monte Carlo es $n=1000$ (esto se ejecuta durante 17 segundos en mi ordenador):

d = 2000;
k = 50;
n_iter = 1000;

tic
for rep = 1:n_iter
    A = randn(d,k);
    [A,~,~] = svd(A,0);  %// orthogonalizing
    B = randn(d,k);
    [B,~,~] = svd(B,0);  %// orthogonalizing

    w(rep) = sum(sum((transpose(A)*B).^2)); %// = trace(A*A'*B*B'), but faster!
end
toc

figure
[f, xi] = ksdensity(w);
h1 = plot(xi, f, 'LineWidth', 2);
hold on
x = min(w):(max(w)-min(w))/100:max(w);
mu = k^2/d;
sigma2 = k^2 * 2*(d-1)/d^2/(d+2);
h2 = plot(x, 1/(sqrt(2*pi*sigma2)) * exp(-(x-mu).^2/2/sigma2), 'r', 'LineWidth', 2);

title(['d = ' num2str(d) ', k = ', num2str(k)])
hh = legend({['Observed density (n = ' num2str(niter) ')'], 'Predicted density'});
legend('boxoff')

Interpretación de $w$

Cosenos de ángulos principales entre los subespacios abarcados por las columnas de $A$ y $B$ vienen dados por los valores singulares de $A^\top B$ . Entonces los cuadrados de estos cosenos vienen dados por los valores propios de $A^\top BB^\top A$ o también de $AA^\top BB^\top$ . Así que, geométricamente, su traza es la suma de los cosenos al cuadrado de los ángulos principales. Si $A=B$ entonces todos los ángulos son cero y la suma de los cosenos al cuadrado es igual a $k$ . Si $A \perp B$ , entonces todos los ángulos son $90^\circ$ y la suma de los cosenos al cuadrado es cero.

Me gusta su enfoque de normalizar $w$ por $k$ , es decir, de tomar $w/k$ como principal medida de similitud. Evidentemente, no puede superar $1$ será igual a $1$ cuando los subespacios coinciden, y será cercano a cero si se eligen al azar. En efecto, $\mathbb E[w]=k^2/d$ lo que significa que $\mathbb E[w/k] = k/d$ . Cuando $k \ll d$ , este es cercano a cero.