/ edit: Dalsze działania teraz możesz użyć irlba :: prcomp_irlba
/ edit: śledzenie mojego własnego posta. irlba
ma 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ść Matrix
funkcji, 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 crossprod
funkcji, 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 means
wektorem 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, irlba
obliczenie 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.
X %*% v %*% diag(d, ncol=length(d))
. Matryca v na dysku svd jest równoważna elementowi „obracania”prcomp
obiektu i /X %*% v
lubX %*% v %*% diag(d, ncol=length(d))
reprezentujex
elementprcomp
obiektu. Spójrzstats:::prcomp.default
.Odpowiedzi:
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.
Przykład
R
get.col
prcomp
źródło
irlba
jest to, że można określićnu
ograniczenie algorytmu do pierwszych n głównych składników, co znacznie zwiększa jego skuteczność i (myślę) pomija obliczenia macierzy XX '.irlba
colMeans
macierz rzadką od macierzy iloczynu kropkowego, a następnie uruchomić irlba na wyniku?R