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

Question

/edición: Para seguir con el seguimiento, ahora se puede utilizar irlba::prcomp_irlba

---

/edición: seguimiento de mi propio post. irlba tiene ahora argumentos "centro" y "escala", que permiten utilizarlo para calcular los componentes principales, por ejemplo

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

---

Tengo un gran y escaso Matrix de 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)

Como esta matriz tiene muchas columnas, me gustaría reducir su dimensionalidad a algo más manejable. Puedo utilizar el excelente paquete irlba para realizar la SVD y devolver los primeros n componentes principales (5 mostrados 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 el PCA hay que centrar la matriz (restar la media de cada columna). Esto es muy difícil de hacer en mi conjunto de datos, y además destruiría la dispersión de la matriz.

¿Qué tan "malo" es realizar la SVD en los datos sin escalar, y alimentar directamente a un algoritmo de aprendizaje automático? ¿Hay alguna forma eficiente de escalar estos datos, preservando la dispersión de la matriz?

---

/edit: A traído a mi atención por B_miner, el "PC" debe ser realmente:

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

Además, creo que la respuesta de whuber debe ser bastante fácil de implementar, a través de la crossprod 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 sé muy bien qué hacer con el means vector antes de restar de M_Mt pero lo publicaré tan pronto como lo resuelva.

---

/edit3: Aquí está la versión modificada del código de whuber, utilizando 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 se fija el número de columnas en 10.000 y el número de componentes principales en 25, el irlba -El PCA basado en la tecnología de la información tarda unos 17 minutos en calcular 50 componentes principales aproximados y consume unos 6 GB de RAM, lo que no está nada mal.

Answer 1

En primer lugar, realmente quieres centrar los datos . Si no, el interpretación geométrica del PCA muestra que el primer componente principal estará cerca del vector de medias y todos los 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, sean bastante incorrectas.

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

Considera entonces dos columnas cualesquiera, $Y$ y $Z$ (cada uno de ellos es un $500000$ -vector; que esta dimensión sea $n$ ). Que sus medios sean $m_Y$ y $m_Z$ respectivamente. Lo que quiere para calcular es, escribiendo $\mathbf{1}$ para el $n$ -vector de $1$ 's,

$$(Y - m_Y\mathbf{1}) \cdot (Z - m_Z\mathbf{1}) = Y\cdot Z - m_Z\mathbf{1}\cdot Y - m_Y\mathbf{1}.Z + m_Z m_Y \mathbf{1}\cdot \mathbf{1}\\ = Y\cdot Z -n (m_Ym_Z),$$

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

Esto permite utilizar técnicas de matrices dispersas para calcular $X X'$ cuyas entradas proporcionan los valores de $Y\cdot Z$ y luego ajustar sus coeficientes en función de la $10000$ significa la columna. El ajuste no debería perjudicar, porque parece poco probable $X X'$ será muy escaso.

---

Ejemplo

Lo siguiente R demuestra este enfoque. Utiliza un stub, get.col que en la práctica puede ser una columna de $X$ a la vez desde una fuente de datos externa, reduciendo así la cantidad de RAM necesaria (con cierto coste en la velocidad de cálculo, por supuesto). Calcula el PCA de dos maneras: a través de la SVD aplicada a la construcción anterior y directamente utilizando prcomp . A continuación, compara el resultado de los dos enfoques. El tiempo de cálculo es de unos 50 segundos para 100 columnas y se escala aproximadamente de forma cuadrática: ¡prepárese para esperar cuando realice la 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.)