Dimensionsreduktion (SVD oder PCA) auf einer großen, dünn besetzten Matrix

31

/ edit: Weitere Folgemaßnahmen können jetzt mit irlba :: prcomp_irlba durchgeführt werden


/ edit: verfolge meinen eigenen Beitrag. irlbaVerfü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 Matrixvon 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 crossprodFunktion, 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 meansVektor tun soll, bevor M_Mtich 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 irlbaPCA auf Basis von-etwa 17 Minuten, um 50 ungefähre Hauptkomponenten zu berechnen, und verbraucht etwa 6 GB RAM, was nicht allzu schlimm ist.

Zach
quelle
Zach, neugierig, ob du das jemals gelöst hast.
B_Miner
@B_Miner: Grundsätzlich habe ich SVD durchgeführt, ohne mich zuerst um das Zentrieren oder Skalieren zu kümmern, weil ich nie einen guten Weg gefunden habe, dies zu tun, ohne meine dünne Matrix in eine dichte Matrix zu konvertieren. Die ursprüngliche Matrix% *% der V-Komponente der DVD ergibt die "Hauptkomponenten". Manchmal erhalte ich bessere Ergebnisse, wenn ich die Eigenwerte "einklappe", z. B. v% *% diag (d), wobei d der Vektor der Eigenwerte aus der SVD ist.
Zach
Behandeln Sie v% *% diag (d) für sich oder immer noch multipliziert mit der ursprünglichen Matrix X (dh X% *% v% *% diag (d))? Es scheint, als ob Sie die u-Matrix als Hauptkomponenten-Punktzahl verwenden.
B_Miner
Ich benutze X %*% v %*% diag(d, ncol=length(d)). Die v - Matrix in der SVD ist äquivalent zu dem „Rotation“ Elemente eines prcompObjekts und X %*% vbzw. X %*% v %*% diag(d, ncol=length(d))repräsentiert das xElement eines prcompObjekts. Werfen Sie einen Blick ein stats:::prcomp.default.
Zach
Ja, X% *% v ist das x-Element von prcomp. Es sieht so aus, als würden Sie bei Verwendung der u-Matrix wie in Ihrer Frage tatsächlich X% *% v% *% diag (1 / d) verwenden.
B_Miner

Antworten:

37

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.

XXX1000010000

Y.Z500000nmY.mZ1n1

(Y.-mY.1)(Z-mZ1)=Y.Z-mZ1Y.-mY.1.Z+mZmY.11=Y.Z-n(mY.mZ),

mY.=1Y./nmZ=1Z/n

XXY.Z10000XX


Beispiel

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
quelle
Vielen Dank für die ausführliche Antwort. Einer der Vorteile von irlbaist, dass Sie angeben können, dass nuder Algorithmus auf die ersten n Hauptkomponenten beschränkt werden soll, was die Effizienz erheblich steigert und (glaube ich) die Berechnung der XX'-Matrix umgeht.
Zach
1
100005000005×1091000010000108irlba
Letzteres nehme ich an. =). Ich muss also das Skalarprodukt für jedes colMeansSpaltenpaar in meiner Dünnschichtmatrix berechnen, das der Dünnschichtmatrix von der Skalarproduktmatrix subtrahieren und dann irlba für das Ergebnis ausführen.
Zach
XXRX
5
Ich habe Code hinzugefügt, um dies zu veranschaulichen.
Whuber