Reducción de dimensionalidad (SVD o PCA) en una matriz grande y dispersa

/ edit: Seguimiento adicional ahora puedes usar irlba :: prcomp_irlba

/ edit: siguiendo mi propio post. irlbaahora tiene argumentos de "centro" y "escala", que le permiten usarlo para calcular componentes principales, por ejemplo:

pc <- M %*% irlba(M, nv=5, nu=0, center=colMeans(M), right_only=TRUE)$v

Tengo una gran variedad Matrixde características que me gustaría utilizar en un algoritmo de aprendizaje automático:

library(Matrix)
set.seed(42)
rows <- 500000
cols <- 10000
i <- unlist(lapply(1:rows, function(i) rep(i, sample(1:5,1))))
j <- sample(1:cols, length(i), replace=TRUE)
M <- sparseMatrix(i, j)

Debido a que esta matriz tiene muchas columnas, me gustaría reducir su dimensionalidad a algo más manejable. Puedo usar el excelente paquete irlba para realizar SVD y devolver los primeros n componentes principales (se muestran 5 aquí; probablemente usaré 100 o 500 en mi conjunto de datos real):

library(irlba)
pc <- irlba(M, nu=5)$u

Sin embargo, he leído que antes de realizar PCA, uno debe centrar la matriz (restar la media de la columna de cada columna). Esto es muy difícil de hacer en mi conjunto de datos, y además destruiría la escasez de la matriz.

¿Qué tan "malo" es realizar SVD en los datos no escalados y alimentarlos directamente en un algoritmo de aprendizaje automático? ¿Hay alguna forma eficiente de escalar estos datos, preservando la escasez de la matriz?

/ edit: A me llamó la atención B_miner, las "PC" realmente deberían ser:

pc <- M %*% irlba(M, nv=5, nu=0)$v 

Además, creo que la respuesta de Whuber debería ser bastante fácil de implementar, a través de la crossprodfunción, que es extremadamente rápida en matrices dispersas:

system.time(M_Mt <- crossprod(M)) # 0.463 seconds
system.time(means <- colMeans(M)) #0.003 seconds

Ahora no estoy muy seguro de qué hacer con el meansvector antes de restarlo M_Mt, pero lo publicaré tan pronto como lo resuelva.

/ edit3: Aquí está la versión modificada del código de whuber, usando operaciones de matriz dispersa para cada paso del proceso. Si puede almacenar toda la matriz dispersa en la memoria, funciona muy rápidamente:

library('Matrix')
library('irlba')
set.seed(42)
m <- 500000
n <- 100
i <- unlist(lapply(1:m, function(i) rep(i, sample(25:50,1))))
j <- sample(1:n, length(i), replace=TRUE)
x <- sparseMatrix(i, j, x=runif(length(i)))

n_comp <- 50
system.time({
  xt.x <- crossprod(x)
  x.means <- colMeans(x)
  xt.x <- (xt.x - m * tcrossprod(x.means)) / (m-1)
  svd.0 <- irlba(xt.x, nu=0, nv=n_comp, tol=1e-10)
})
#user  system elapsed 
#0.148   0.030   2.923 

system.time(pca <- prcomp(x, center=TRUE))
#user  system elapsed 
#32.178   2.702  12.322

max(abs(pca$center - x.means))
max(abs(xt.x - cov(as.matrix(x))))
max(abs(abs(svd.0$v / pca$rotation[,1:n_comp]) - 1))

Si establece el número de columnas en 10,000 y el número de componentes principales en 25, el irlbaPCA basado en la computadora demora aproximadamente 17 minutos para calcular 50 componentes principales aproximados y consume alrededor de 6GB de RAM, lo cual no es tan malo.

En primer lugar, realmente desea centrar los datos . De lo contrario, la interpretación geométrica de PCA muestra que el primer componente principal estará cerca del vector de medios y todas las PC posteriores serán ortogonales a él, lo que evitará que se aproximen a cualquier PC que esté cerca de ese primer vector. Podemos esperar que la mayoría de las PC posteriores sean aproximadamente correctas, pero el valor de eso es cuestionable cuando es probable que las primeras PC, las más importantes, estén bastante equivocadas.

¿Entonces lo que hay que hacer? PCA procede por medio de una descomposición de valor singular de la matriz . La información esencial estará contenida en X X $X$$X X'$ , que en este caso es una matriz de por 10000 : eso puede ser manejable. Su cálculo implica unos 50 millones de cálculos de productos de punto de una columna con la siguiente.$10000$$10000$

Considere dos columnas, entonces, y Z (cada una de ellas es un vector 500000 ; deje que esta dimensión sea n ). Deje que sus medios sean m Y y m Z , respectivamente. Lo que quieres calcular es escribir $Y$$Z$$500000$$n$$m_Y$$m_Z$$\mathbf{1}$$n$$1$

$m_Y = \mathbf{1}\cdot Y / n$$m_Z = \mathbf{1}\cdot Z/n$

$X X'$$Y\cdot Z$$10000$$X X'$

Ejemplo

Rget.col$X$ a la vez desde una fuente de datos externa, reduciendo así la cantidad de RAM requerida (a cierto costo en la velocidad de cálculo, por supuesto). Calcula PCA de dos maneras: a través de SVD aplicado a la construcción anterior y directamente usando prcomp. Luego compara la salida de los dos enfoques. El tiempo de cálculo es de aproximadamente 50 segundos para 100 columnas y escalas aproximadamente cuadráticamente: ¡prepárese para esperar cuando realice SVD en una matriz de 10K por 10K!

m <- 500000 # Will be 500,000
n <- 100    # will be 10,000
library("Matrix")
x <- as(matrix(pmax(0,rnorm(m*n, mean=-2)), nrow=m), "sparseMatrix")
#
# Compute centered version of x'x by having at most two columns
# of x in memory at any time.
#
get.col <- function(i) x[,i] # Emulates reading a column
system.time({
  xt.x <- matrix(numeric(), n, n)
  x.means <- rep(numeric(), n)
  for (i in 1:n) {
    i.col <- get.col(i)
    x.means[i] <- mean(i.col)
    xt.x[i,i] <- sum(i.col * i.col)
    if (i < n) {
      for (j in (i+1):n) {
        j.col <- get.col(j)
        xt.x[i,j] <- xt.x[j,i] <- sum(j.col * i.col)
      }    
    }
  }
  xt.x <- (xt.x - m * outer(x.means, x.means, `*`)) / (m-1)
  svd.0 <- svd(xt.x / m)
}
)
system.time(pca <- prcomp(x, center=TRUE))
#
# Checks: all should be essentially zero.
#
max(abs(pca$center - x.means))
max(abs(xt.x - cov(x)))
max(abs(abs(svd.0$v / pca$rotation) - 1)) # (This is an unstable calculation.)