Redukcja wymiarów (SVD lub PCA) na dużej, rzadkiej matrycy

31

/ edit: Dalsze działania teraz możesz użyć irlba :: prcomp_irlba


/ edit: śledzenie mojego własnego posta. irlbama teraz argumenty „środkowy” i „skalowany”, które pozwalają go używać do obliczania podstawowych składników, np .:

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


Mam dużą różnorodność Matrixfunkcji, których chciałbym użyć w algorytmie uczenia maszynowego:

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)

Ponieważ ta macierz ma wiele kolumn, chciałbym sprowadzić jej wymiar do czegoś łatwiejszego w zarządzaniu. Mogę użyć doskonałego pakietu irlba do wykonania SVD i zwrócenia pierwszych n głównych komponentów (5 tutaj pokazanych; prawdopodobnie użyję 100 lub 500 w moim rzeczywistym zestawie danych):

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

Jednak przeczytałem, że przed wykonaniem PCA należy wyśrodkować macierz (odjąć średnią z każdej kolumny). Jest to bardzo trudne do wykonania w moim zestawie danych, a ponadto zniszczyłoby rzadkość macierzy.

Jak „źle” jest wykonywać SVD na nieskalowanych danych i wprowadzać je bezpośrednio do algorytmu uczenia maszynowego? Czy istnieją wydajne sposoby skalowania tych danych przy jednoczesnym zachowaniu rzadkości macierzy?


/ edit: A, na który zwrócił moją uwagę B_miner, „komputery” powinny naprawdę być:

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

Myślę też, że odpowiedź Whubera powinna być dość łatwa do wdrożenia, dzięki crossprodfunkcji, która jest niezwykle szybka na rzadkich macierzach:

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

Teraz nie jestem pewien, co zrobić z meanswektorem przed odjęciem M_Mt, ale opublikuję, gdy tylko go wymyślę.


/ edit3: Oto zmodyfikowana wersja kodu Whubera, wykorzystująca rzadkie operacje macierzowe na każdym etapie procesu. Jeśli możesz zapisać całą rzadką macierz w pamięci, działa ona bardzo szybko:

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))

Jeśli ustawisz liczbę kolumn na 10 000, a liczbę głównych komponentów na 25, irlbaobliczenie 50 głównych komponentów na podstawie PCA zajmie około 17 minut i zużyje około 6 GB pamięci RAM, co nie jest takie złe.

Zach
źródło
Zach, ciekawy, czy kiedykolwiek to rozwiązałeś.
B_Miner
@B_Miner: Zasadniczo robiłem SVD, nie przejmując się najpierw wyśrodkowaniem lub skalowaniem, ponieważ nigdy nie znalazłem dobrego sposobu na zrobienie tego bez konwersji mojej rzadkiej macierzy do gęstej macierzy. Oryginalna macierz% *% komponentu V svd daje „podstawowe składniki”. Czasami uzyskuję lepsze wyniki, jeśli „zwijam” wartości własne, np. V% *% diag (d), gdzie d jest wektorem wartości własnych z SVD.
Zach
Czy traktujesz v% *% diag (d) jako taki, czy nadal pomnożony przez oryginalną macierz X (tj. X% *% v% *% diag (d)). Wygląda na to, że używasz macierzy u jako głównej oceny wyników?
B_Miner
Używam X %*% v %*% diag(d, ncol=length(d)). Matryca v na dysku svd jest równoważna elementowi „obracania” prcompobiektu i / X %*% vlub X %*% v %*% diag(d, ncol=length(d))reprezentuje xelement prcompobiektu. Spójrz stats:::prcomp.default.
Zach
Tak, X% *% v jest elementem x z prcomp. Wygląda na to, że kiedy używasz macierzy u jak w pytaniu, faktycznie używasz X% *% v% *% diag (1 / d).
B_Miner

Odpowiedzi:

37

Przede wszystkim naprawdę chcesz wyśrodkować dane . Jeśli nie, interpretacja geometryczna PCA pokazuje, że pierwszy główny składnik będzie zbliżony do wektora średnich, a wszystkie kolejne PC będą do niego ortogonalne, co uniemożliwi im zbliżenie się do dowolnych komputerów, które są blisko tego pierwszego wektora. Możemy mieć nadzieję, że większość późniejszych komputerów będzie w przybliżeniu poprawna, ale ich wartość jest wątpliwa, gdy prawdopodobne jest, że kilka pierwszych komputerów - najważniejszych - będzie zupełnie błędnych.

XXX1000010000

YZ500000nmYmZ1n1

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

mY=1Y/nmZ=1Z/n

XXYZ10000XX


Przykład

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
źródło
Dziękuję za szczegółową odpowiedź. Jedną z zalet irlbajest to, że można określić nuograniczenie algorytmu do pierwszych n głównych składników, co znacznie zwiększa jego skuteczność i (myślę) pomija obliczenia macierzy XX '.
Zach.
1
100005000005×1091000010000108irlba
Podejrzewam, że to drugie. =). Więc muszę obliczyć iloczyn kropkowy dla każdej pary kolumn w mojej rzadkiej macierzy, odjąć colMeansmacierz rzadką od macierzy iloczynu kropkowego, a następnie uruchomić irlba na wyniku?
Zach.
XXRX
5
Dodałem kod do zilustrowania.
whuber