Pengurangan dimensi (SVD atau PCA) pada matriks besar dan jarang

/ edit: Tindak lanjut lebih lanjut sekarang Anda dapat menggunakan irlba :: prcomp_irlba

/ edit: menindaklanjuti posting saya sendiri. irlbasekarang memiliki argumen "pusat" dan "skala", yang memungkinkan Anda menggunakannya untuk menghitung komponen prinsip, misalnya:

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

Saya memiliki banyak Matrixfitur yang ingin saya gunakan dalam algoritma pembelajaran mesin:

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)

Karena matriks ini memiliki banyak kolom, saya ingin mengurangi dimensinya menjadi sesuatu yang lebih mudah dikelola. Saya dapat menggunakan paket irlba yang sangat baik untuk melakukan SVD dan mengembalikan komponen utama n pertama (5 ditampilkan di sini; Saya mungkin akan menggunakan 100 atau 500 pada dataset aktual saya):

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

Namun, saya telah membaca bahwa sebelum melakukan PCA, orang harus memusatkan matriks (kurangi rata-rata kolom dari setiap kolom). Ini sangat sulit dilakukan pada dataset saya, dan selanjutnya akan menghancurkan sparsity dari matriks.

Seberapa "buruk" untuk melakukan SVD pada data yang tidak diskalakan, dan memasukkannya langsung ke algoritma pembelajaran mesin? Adakah cara efisien yang bisa saya lakukan untuk menskala data ini, sambil menjaga sparsity dari matriks?

/ edit: A dibawa ke perhatian saya oleh B_miner, "PC" harus benar-benar:

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

Juga, saya pikir jawaban whuber seharusnya cukup mudah diimplementasikan, melalui crossprodfungsi, yang sangat cepat pada matriks jarang:

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

Sekarang saya tidak yakin apa yang harus dilakukan terhadap meansvektor sebelum dikurangi M_Mt, tetapi akan dikirim segera setelah saya mengetahuinya.

/ edit3: Berikut ini adalah versi modifikasi dari kode whuber, menggunakan operasi matriks jarang untuk setiap langkah proses. Jika Anda dapat menyimpan seluruh matriks jarang dalam memori, ia bekerja sangat cepat:

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

Jika Anda mengatur jumlah kolom menjadi 10.000 dan jumlah komponen utama menjadi 25, irlbaPCA berbasiskan membutuhkan waktu sekitar 17 menit untuk menghitung 50 perkiraan komponen utama dan mengkonsumsi sekitar 6GB RAM, yang tidak terlalu buruk.

Pertama-tama, Anda benar-benar ingin memusatkan data . Jika tidak, interpretasi geometris PCA menunjukkan bahwa komponen utama pertama akan dekat dengan vektor rata-rata dan semua PC berikutnya akan ortogonal, yang akan mencegah mereka mendekati PC mana pun yang dekat dengan vektor pertama tersebut. Kita dapat berharap bahwa sebagian besar PC yang akan datang kira-kira benar, tetapi nilai yang dipertanyakan ketika kemungkinan beberapa PC pertama - yang paling penting - akan sangat salah.

Jadi, apa yang harus dilakukan? PCA hasil dengan cara dekomposisi nilai singular dari matriks . Informasi penting akan terkandung dalam X X ' , yang dalam hal ini adalah 10000 dengan 10000 matriks: yang mungkin dikelola. Perhitungannya melibatkan sekitar 50 juta perhitungan produk titik dari satu kolom dengan yang berikutnya.$X$$X X'$$10000$$10000$

Pertimbangkan dua kolom, lalu, dan Z (masing-masingnya adalah vektor 500000 ; biarkan dimensi ini n ). Biarkan sarana mereka menjadi m Y dan m Z , masing-masing. Yang ingin Anda hitung adalah, menulis 1 untuk n- vektor 1 's,$Y$$Z$$500000$$n$$m_Y$$m_Z$$\mathbf{1}$$n$$1$

karena dan m Z = 1 ⋅ Z / n .$m_Y = \mathbf{1}\cdot Y / n$$m_Z = \mathbf{1}\cdot Z/n$

Hal ini memungkinkan Anda untuk menggunakan teknik matriks jarang untuk menghitung , yang entri memberikan nilai-nilai Y ⋅ Z , dan kemudian menyesuaikan koefisien yang didasarkan pada 10000 sarana kolom. Penyesuaian seharusnya tidak sakit, karena sepertinya tidak mungkin X X ′ akan sangat jarang.$X X'$$Y\cdot Z$$10000$$X X'$

Contoh

RKode berikut menunjukkan pendekatan ini. Ini menggunakan rintisan,, get.colyang dalam praktiknya mungkin membaca satu kolom pada suatu waktu dari sumber data eksternal, sehingga mengurangi jumlah RAM yang diperlukan (tentu saja dengan biaya dalam kecepatan komputasi, tentu saja). Ini menghitung PCA dalam dua cara: melalui SVD diterapkan pada konstruksi sebelumnya dan langsung menggunakan . Kemudian membandingkan output dari dua pendekatan. Waktu perhitungan adalah sekitar 50 detik untuk 100 kolom dan skala kira-kira secara kuadratik: bersiaplah untuk menunggu ketika melakukan SVD pada matriks 10K x 10K!$X$prcomp

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