/ edit: Plus de suivi maintenant vous pouvez utiliser irlba :: prcomp_irlba
/ edit: suivi de mon propre post. irlba
a 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 Matrix
de 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 crossprod
fonction, 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 means
vecteur 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 irlba
PCA 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.
X %*% v %*% diag(d, ncol=length(d))
. La matrice v dans le svd est équivalente à l'élément "rotation" d'unprcomp
objet etX %*% v
ouX %*% v %*% diag(d, ncol=length(d))
représente l'x
élément d'unprcomp
objet. Jetez un oeil astats:::prcomp.default
.Réponses:
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.
Exemple
R
get.col
prcomp
la source
irlba
est que vous pouvez spécifiernu
de 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 '.irlba
colMeans
matrice clairsemée de la matrice du produit scalaire, puis exécuter irlba sur le résultat?R