/ edit: далее следите, теперь вы можете использовать irlba :: prcomp_irlba
/ edit: следите за своим собственным постом. irlba
теперь имеет аргументы "center" и "scale", которые позволяют использовать его для вычисления основных компонентов, например:
pc <- M %*% irlba(M, nv=5, nu=0, center=colMeans(M), right_only=TRUE)$v
У меня есть большой набор Matrix
функций, которые я хотел бы использовать в алгоритме машинного обучения:
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)
Поскольку эта матрица имеет много столбцов, я хотел бы уменьшить ее размерность до чего-то более управляемого. Я могу использовать отличный пакет irlba для выполнения SVD и вернуть первые n основных компонентов (5 показано здесь; я, вероятно, буду использовать 100 или 500 в моем фактическом наборе данных):
library(irlba)
pc <- irlba(M, nu=5)$u
Тем не менее, я прочитал, что перед выполнением PCA, необходимо отцентрировать матрицу (вычесть среднее значение столбца из каждого столбца). Это очень сложно сделать с моим набором данных, и, кроме того, это приведет к разрушению разреженности матрицы.
Насколько «плохо» выполнять SVD на немасштабированных данных и подавать их прямо в алгоритм машинного обучения? Есть ли эффективные способы, которыми я мог бы масштабировать эти данные, сохраняя разреженность матрицы?
/ edit: A, на мой взгляд B_miner, "ПК" должны быть:
pc <- M %*% irlba(M, nv=5, nu=0)$v
Кроме того, я думаю, что ответ whuber должен быть довольно простым для реализации с помощью crossprod
функции, которая очень быстро работает с разреженными матрицами:
system.time(M_Mt <- crossprod(M)) # 0.463 seconds
system.time(means <- colMeans(M)) #0.003 seconds
Теперь я не совсем уверен, что делать с means
вектором, прежде чем вычесть M_Mt
, но опубликую, как только я это выясню.
/ edit3: Вот модифицированная версия кода whuber, использующая разреженные матричные операции для каждого шага процесса. Если вы можете хранить всю разреженную матрицу в памяти, она работает очень быстро:
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))
Если вы установите число столбцов равным 10 000, а число основных компонентов - 25, irlba
PCA на основе вычислений займет около 17 минут, чтобы рассчитать 50 приблизительных основных компонентов, и потребует около 6 ГБ ОЗУ, что не так уж и плохо.
X %*% v %*% diag(d, ncol=length(d))
. Матрица v в svd эквивалентна элементу «вращения»prcomp
объекта иX %*% v
илиX %*% v %*% diag(d, ncol=length(d))
представляетx
элементprcomp
объекта. Посмотриstats:::prcomp.default
.Ответы:
Прежде всего, вы действительно хотите центрировать данные . Если нет, то геометрическая интерпретация PCA показывает, что первый главный компонент будет близок к вектору средних значений, и все последующие ПК будут ортогональны ему, что не позволит им приблизиться к любым ПК, оказавшимся близко к этому первому вектору. Мы можем надеяться, что большинство более поздних ПК будут примерно правильными, но ценность этого сомнительна, когда, вероятно, первые несколько ПК - самые важные - будут совершенно неверными.
пример
R
get.col
prcomp
источник
irlba
является то, что вы можете указатьnu
ограничение алгоритма первыми n основными компонентами, что значительно повышает его эффективность и (я думаю) обходит вычисление матрицы XX '.irlba
colMeans
разреженную матрицу из матрицы точечных произведений, а затем запустить irlba для результата?R