机器算法验证 - 大型稀疏矩阵的降维（SVD 或 PCA） - 吾爱随笔录

大型稀疏矩阵的降维（SVD 或 PCA）

机器算法验证 r 主成分分析降维 svd 矩阵分解

2022-01-17 23:16:15

/edit: 现在进一步跟进你可以使用irlba::prcomp_irlba

/编辑：跟进我自己的帖子。 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：B_miner 引起了我的注意，“PC”应该是：

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 个近似主成分，并且消耗大约 6GB 的 RAM，这还不错。

1个回答

首先，您确实想将数据居中。如果不是，PCA 的几何解释表明第一个主成分将接近均值向量，并且所有后续 PC 将与其正交，这将阻止它们逼近任何恰好接近该第一个向量的 PC。我们可以希望后来的大多数 PC 大致正确，但是当最初的几台 PC（最重要的 PC）可能完全错误时，其价值值得怀疑。

那么该怎么办？PCA 通过矩阵的奇异值分解进行 $X$ . 基本信息将包含在 $X X'$ , 在这种情况下是 $10000$ 经过 $10000$ 矩阵：这可能是可以管理的。它的计算涉及一列与下一列的点积的大约 5000 万次计算。

考虑任意两列，然后， $Y$ 和 $Z$ （每一个都是一个 $500000$ -向量; 让这个维度成为 $n$ ）。让他们的手段 $m_Y$ 和 $m_Z$ ，分别。你要计算的是，写 $\mathbf{1}$ 为了 $n$ -向量 $1$ 的，

(Y - m_{Y} 1) \cdot (Z - m_{Z} 1) = Y \cdot Z - m_{Z} 1 \cdot Y - m_{Y} 1 . Z + m_{Z} m_{Y} 1 \cdot 1 = Y \cdot Z - n (m_{Y} m_{Z}),

$(Y - m_Y\mathbf{1}) \cdot (Z - m_Z\mathbf{1}) = Y\cdot Z - m_Z\mathbf{1}\cdot Y - m_Y\mathbf{1}.Z + m_Z m_Y \mathbf{1}\cdot \mathbf{1}\\ = Y\cdot Z -n (m_Ym_Z),$

因为 $m_Y = \mathbf{1}\cdot Y / n$ 和 $m_Z = \mathbf{1}\cdot Z/n$ .

这允许您使用稀疏矩阵技术来计算 $X X'$ ，其条目提供的值 $Y\cdot Z$ , 然后根据 $10000$ 列的意思。调整不应该受到伤害，因为它似乎不太可能 $X X'$ 会很稀疏。

例子

下面的R代码演示了这种方法。它使用一个存根，get.col在实践中可能会读取一列 $X$ 一次从外部数据源，从而减少所需的 RAM 量（当然，以计算速度为代价）。它以两种方式计算 PCA：通过应用于前面构造的 SVD 和直接使用prcomp. 然后比较两种方法的输出。100 列的计算时间大约为 50 秒，并且大约按二次缩放：在 10K x 10K 矩阵上执行 SVD 时，请准备好等待！

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

其它你可能感兴趣的问题

上一篇为什么需要变量选择？下一篇为什么整流线性单元被认为是非线性的？