Réduction dimensionnelle (SVD ou PCA) sur une grande matrice clairsemée

31

/ edit: Plus de suivi maintenant vous pouvez utiliser irlba :: prcomp_irlba


/ edit: suivi de mon propre post. irlbaa maintenant des arguments "center" et "scale", qui vous permettent de calculer les principaux composants, par exemple:

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


J'ai un grand nombre Matrixde fonctionnalités que j'aimerais utiliser dans un algorithme d'apprentissage automatique:

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)

Parce que cette matrice a de nombreuses colonnes, je voudrais réduire sa dimensionnalité à quelque chose de plus gérable. Je peux utiliser l'excellent package irlba pour effectuer SVD et retourner les n premiers composants principaux (5 affichés ici; j'utiliserai probablement 100 ou 500 sur mon ensemble de données réel):

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

Cependant, j'ai lu qu'avant d'effectuer l'ACP, il faut centrer la matrice (soustraire la moyenne des colonnes de chaque colonne). C'est très difficile à faire sur mon jeu de données, et cela détruirait en outre la rareté de la matrice.

À quel point est-il «mauvais» d'exécuter SVD sur les données non mises à l'échelle et de les alimenter directement dans un algorithme d'apprentissage automatique? Existe-t-il des moyens efficaces de mettre à l'échelle ces données tout en préservant la rareté de la matrice?


/ edit: A porté à mon attention par B_miner, les "PC" devraient vraiment être:

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

De plus, je pense que la réponse de whuber devrait être assez facile à implémenter, via la crossprodfonction, qui est extrêmement rapide sur les matrices clairsemées:

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

Maintenant, je ne sais pas trop quoi faire du meansvecteur avant de le soustraire M_Mt, mais je posterai dès que je le comprendrai.


/ edit3: Voici une version modifiée du code de whuber, utilisant des opérations matricielles clairsemées pour chaque étape du processus. Si vous pouvez stocker toute la matrice clairsemée en mémoire, cela fonctionne très rapidement:

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 vous définissez le nombre de colonnes sur 10 000 et le nombre de composants principaux sur 25, le irlbaPCA basé prend environ 17 minutes pour calculer 50 composants principaux approximatifs et consomme environ 6 Go de RAM, ce qui n'est pas trop mal.

Zach
la source
Zach, curieux de savoir si vous avez déjà résolu cela.
B_Miner
@B_Miner: Fondamentalement, je fais de la SVD sans prendre la peine de centrer ou de mettre à l'échelle d'abord, car je n'ai jamais trouvé un bon moyen de le faire sans convertir ma matrice clairsemée en une matrice dense. La matrice d'origine% *% le composant V du svd donne les "composants principaux". Parfois, j'obtiens de meilleurs résultats si je "replie" les valeurs propres, par exemple v% *% diag (d), où d est le vecteur des valeurs propres du SVD.
Zach
Traitez-vous v% *% diag (d) par lui-même ou encore multiplié par la matrice d'origine X (c'est-à-dire X% *% v% *% diag (d)). Il semble que ci-dessus, vous utilisez la matrice u comme score de composante principale?
B_Miner
J'utilise X %*% v %*% diag(d, ncol=length(d)). La matrice v dans le svd est équivalente à l'élément "rotation" d'un prcompobjet et X %*% vou X %*% v %*% diag(d, ncol=length(d))représente l' xélément d'un prcompobjet. Jetez un oeil a stats:::prcomp.default.
Zach
Oui, X% *% v est l'élément x de prcomp. Il semble que lorsque vous utilisez la matrice u comme dans votre question, vous utilisez en fait X% *% v% *% diag (1 / d).
B_Miner

Réponses:

37

Tout d'abord, vous voulez vraiment centrer les données . Sinon, l' interprétation géométrique de l'ACP montre que la première composante principale sera proche du vecteur de moyennes et que tous les PC suivants lui seront orthogonaux, ce qui les empêchera d'approximer les PC qui se trouvent être proches de ce premier vecteur. Nous pouvons espérer que la plupart des PC ultérieurs seront approximativement corrects, mais la valeur de cela est discutable alors qu'il est probable que les premiers PC - les plus importants - seront tout à fait faux.

XXX1000010000

YZ500000nmYmZ1n1

(YmY1)(ZmZ1)=YZmZ1YmY1.Z+mZmY11=YZn(mYmZ),

mY=1Y/nmZ=1Z/n

XXYZ10000XX


Exemple

Rget.colXprcomp

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.)
whuber
la source
Merci pour la réponse détaillée. L'un des avantages de irlbaest que vous pouvez spécifier nude limiter l'algorithme aux n premiers composants principaux, ce qui augmente considérablement son efficacité et (je pense) contourne le calcul de la matrice XX '.
Zach
1
100005000005×1091000010000108irlba
Je suppose que ce dernier. =). Je dois donc calculer le produit scalaire pour chaque paire de colonnes dans ma matrice clairsemée, soustraire la colMeansmatrice clairsemée de la matrice du produit scalaire, puis exécuter irlba sur le résultat?
Zach
XXRX
5
J'ai ajouté du code pour illustrer.
whuber