¿Cómo calcular la matriz de covarianza tridiagonal aproximada para una descorrelación rápida?

Question

Dada una matriz de datos $X$ de digamos 1000000 observaciones $\times$ 100 prestaciones, ¿hay una forma rápida de construir una aproximación tridiagonal $A \approx cov(X)$ ?
Entonces uno podría factorizar $A = L L^T$ , $L$ todos 0 excepto $L_{i\ i-1}$ y $L_{i i}$ , y hacer una descorrelación rápida (blanqueamiento) resolviendo $L x = x_{white}$ . (Con "rápido" me refiero a $O( size\ X )$ .)

(Añadido, intentando aclarar): Estoy buscando un blanqueador rápido y sucio que sea más rápido que el completo $cov(X)$ pero mejor que en diagonal. Digamos que $X$ es $N$ puntos de datos $\times Nf$ características, por ejemplo 1000000 $\times$ 100, con características 0-media.

1) construir $Fullcov = X^T X$ factor Cholesky como $L L^T$ , resolver $L x = x_{white}$ para blanquear nuevos $x$ s. Esto es cuadrático en el número de características.

2) diagonal: $x_{white} = x / \sigma(x)$ ignora por completo las correlaciones cruzadas.

Un podría obtener una matriz tridiagonal a partir de $Fullcov$ simplemente poniendo a cero todas las entradas fuera de la tridiagonal, o no acumulándolas en primer lugar. Y aquí empiezo a hundirme: debe haber una aproximación mejor, ¿tal vez jerárquica, diagonal en bloque → tridiagonal?

---

(Añadido el 11 de mayo): Permítanme dividir la pregunta en dos:

1) ¿existe una aproximación rápida $cov(X)$ ?
No (whuber), hay que mirar todos ${N \choose 2}$ pares (o tener estructura, o muestra).

2) dado un $cov(X)$ ¿con qué rapidez se puede blanquear nuevo $x$ s ?
Bueno, teniendo en cuenta $cov = L L^T$ , $L$ triangular inferior, una vez, luego resolviendo $L x = x_{white}$ es bastante rápido; scipy.linalg.solve_triangular, por ejemplo, utiliza Lapack.
Estaba buscando un blanqueador aún más rápido, sigo buscando.

Answer 1

Simplemente informática La matriz de covarianza, que necesitarás para empezar, es la siguiente $O((Nf)^2)$ por lo que, asintóticamente en $N$ no se gana nada eligiendo un $O(Nf)$ algoritmo para el blanqueamiento.

Existen aproximaciones cuando las variables tienen una estructura adicional, como cuando forman una serie temporal o realizaciones de un proceso estocástico espacial en varias ubicaciones. Estas aproximaciones se basan en supuestos que nos permiten relacionar la covarianza entre un par de variables con la covarianza entre otros pares de variables, por ejemplo entre pares separados por los mismos retardos. Esta es la razón convencional para suponer que un proceso es estacionario o intrínsecamente estacionario por ejemplo. Los cálculos pueden ser $O(Nf\,\log(Nf)$ en tales casos ( Por ejemplo utilizando la transformada rápida de Fourier como en Yao y Journel 1998 ). En ausencia de tal modelo, no veo cómo se puede evitar el cálculo de todas las covarianzas entre pares.

Answer 2

En un capricho, decidí tratar de calcular (en R) la matriz de covarianza para un conjunto de datos de aproximadamente el tamaño mencionado en el OP:

z <- rnorm(1e8)
dim(z) <- c(1e6, 100)
vcv <- cov(z)

En total se tardó menos de un minuto, en un portátil bastante genérico con Windows XP de 32 bits. Probablemente se tardó más en generar z en primer lugar que calcular la matriz vcv . Y R no está especialmente optimizado para operaciones matriciales.

Teniendo en cuenta este resultado, ¿es tan importante la velocidad? Si N >> p, el tiempo necesario para calcular la aproximación probablemente no será mucho menor que para obtener la matriz de covarianza real.

Answer 3

Dos céntimos extra:

Algorítmicamente hablando, no creo que haya algoritmos más rápidos para hacer esto de forma genérica. $X$ . Si las hubiera, ya se habrían implementado en los programas hasta ahora. Sin embargo, desde una perspectiva de ingeniería de software, las velocidades pueden diferir drásticamente entre implementaciones (por ejemplo, Blas heredado, Goto Blas, MKL de Intel y OpenBlas).

Dependiendo del escenario de aplicación, puedes codificar tus implementaciones específicas para tus casos de uso, pero esto requiere hacer malabarismos de código en un lenguaje más nativo como Fortran, C y C++. Ahora, para hacer implementaciones rápidas, una cosa que definitivamente se necesita es aprovechar las nuevas características de la CPU, como AVX512, que requiere un poco de conocimiento de ASM.

Además, dependiendo de lo rudimentario que sea, una posible aproximación es utilizar la integración de Montecarlo. Por ejemplo, para una columna de 1000000 x 1 $x$ (digamos, 1000000 observaciones de algo), el cálculo literal de la media será $\sum_{i=1}^{1000000}x_i/1000000$ pero para una aproximación, se puede tomar una submuestra de sólo 1000 para obtener una estimación aproximada de la media, por ejemplo $\sum_{j=1}^{1000}x_{i(j)}/1000$ donde $x_{i(j)}$ es una submuestra de la x original. La idea se aplica al cálculo de su covarianza. Digamos que tiene dos columnas x e y; en lugar de hacer $\sum_{i=1}^{1000000}x_i*y_i$ se puede aproximar con sólo una submuestra (por ejemplo, 1000) $\sum_{j=1}^{1000}x_{i(j)}*y_{i(j)}/1000 *1000000 $ . Si la ordenación de sus filas es suficientemente aleatoria, no necesita muestrear explícitamente estos 1000 índices, puede hacerlo por un paso regular o simplemente tomar las 1000 primeras filas.

Para que tu sustitución inversa resuelva $Lx=x_w$ si L se reutiliza muchas veces, una pequeña mejora consiste en almacenar explícitamente los elementos diagonales de $L$ como su inversión (por ejemplo $1/Lii$ ); de este modo se evitará hacer la división y se hará la multiplicación en la etapa de sustitución inversa. Hoy en día, este ahorro no tiene importancia, ya que la velocidad de división y multiplicación de flotantes debería ser aproximadamente la misma. Pero este no es el caso de los antiguos CPUS porque la división solía ser mucho más lenta que la multiplicación. Una vez más, esto es sólo una cosa menor y no aceleraría el cálculo de forma espectacular, pero dependiendo de los casos de usuario, puede ayudar un poco.