/ edit: Weitere Folgemaßnahmen können jetzt mit irlba :: prcomp_irlba durchgeführt werden
/ edit: verfolge meinen eigenen Beitrag. irlba
Verfügt nun über die Argumente "center" und "scale", mit denen Sie Hauptkomponenten berechnen können, z.
pc <- M %*% irlba(M, nv=5, nu=0, center=colMeans(M), right_only=TRUE)$v
Ich habe eine große, spärliche Anzahl Matrix
von Funktionen, die ich in einem Algorithmus für maschinelles Lernen verwenden möchte:
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)
Da diese Matrix viele Spalten hat, möchte ich ihre Dimensionalität auf etwas übersichtlicheres reduzieren. Ich kann das ausgezeichnete irlba-Paket verwenden , um SVD durchzuführen und die ersten n Hauptkomponenten zurückzugeben (5 hier gezeigt; ich werde wahrscheinlich 100 oder 500 für meinen tatsächlichen Datensatz verwenden):
library(irlba)
pc <- irlba(M, nu=5)$u
Ich habe jedoch gelesen, dass man vor der Durchführung einer PCA die Matrix zentrieren sollte (subtrahieren Sie den Spaltenmittelwert von jeder Spalte). Dies ist bei meinem Datensatz sehr schwierig und würde außerdem die Sparsamkeit der Matrix zerstören.
Wie "schlecht" ist es, SVD für die nicht skalierten Daten durchzuführen und sie direkt in einen Algorithmus für maschinelles Lernen einzuspeisen? Gibt es effiziente Möglichkeiten, diese Daten zu skalieren, während die Sparsamkeit der Matrix erhalten bleibt?
/ edit: Ein von B_miner an mich herangetragener "PC" sollte eigentlich sein:
pc <- M %*% irlba(M, nv=5, nu=0)$v
Außerdem denke ich, dass die Antwort von whuber über die crossprod
Funktion, die bei spärlichen Matrizen extrem schnell ist, ziemlich einfach zu implementieren sein sollte :
system.time(M_Mt <- crossprod(M)) # 0.463 seconds
system.time(means <- colMeans(M)) #0.003 seconds
Jetzt bin ich mir nicht ganz sicher, was ich mit dem means
Vektor tun soll, bevor M_Mt
ich ihn subtrahiere , werde aber posten, sobald ich es herausgefunden habe.
/ edit3: Hier ist eine modifizierte Version von Whubers Code, die spärliche Matrixoperationen für jeden Schritt des Prozesses verwendet. Wenn Sie die gesamte dünne Matrix im Speicher ablegen können, funktioniert dies sehr schnell:
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))
Wenn Sie die Anzahl der Spalten auf 10.000 und die Anzahl der Hauptkomponenten auf 25 festlegen, benötigt der irlba
PCA auf Basis von-etwa 17 Minuten, um 50 ungefähre Hauptkomponenten zu berechnen, und verbraucht etwa 6 GB RAM, was nicht allzu schlimm ist.
X %*% v %*% diag(d, ncol=length(d))
. Die v - Matrix in der SVD ist äquivalent zu dem „Rotation“ Elemente einesprcomp
Objekts undX %*% v
bzw.X %*% v %*% diag(d, ncol=length(d))
repräsentiert dasx
Element einesprcomp
Objekts. Werfen Sie einen Blick einstats:::prcomp.default
.Antworten:
Zunächst möchten Sie die Daten wirklich zentrieren . Wenn nicht, zeigt die geometrische Interpretation von PCA , dass sich die erste Hauptkomponente in der Nähe des Mittelwertvektors befindet und alle nachfolgenden PCs orthogonal dazu sind, was sie daran hindert, sich etwaigen PCs anzunähern, die sich diesem ersten Vektor nähern. Wir können hoffen, dass die meisten späteren PCs annähernd korrekt sind, aber der Wert ist fraglich, wenn es sich wahrscheinlich um die ersten mehreren PCs handelt - die wichtigsten -, die völlig falsch sind.
Beispiel
R
get.col
prcomp
quelle
irlba
ist, dass Sie angeben können, dassnu
der Algorithmus auf die ersten n Hauptkomponenten beschränkt werden soll, was die Effizienz erheblich steigert und (glaube ich) die Berechnung der XX'-Matrix umgeht.irlba
colMeans
Spaltenpaar in meiner Dünnschichtmatrix berechnen, das der Dünnschichtmatrix von der Skalarproduktmatrix subtrahieren und dann irlba für das Ergebnis ausführen.R