在求解最小二乘估计的大正规方程时使用“outer”时内存不足
Out of memory when using `outer` in solving my big normal equation for least squares estimation
考虑以下 R 中的示例:
x1 <- rnorm(100000)
x2 <- rnorm(100000)
g <- cbind(x1, x2, x1^2, x2^2)
gg <- t(g) %*% g
gginv <- solve(gg)
bigmatrix <- outer(x1, x2, "<=")
Gw <- t(g) %*% bigmatrix
beta <- gginv %*% Gw
w1 <- bigmatrix - g %*% beta
如果我尝试在我的计算机中 运行 这样的东西,它会抛出内存错误(因为 bigmatrix
太大了)。
你知道我怎样才能达到同样的效果,而不 运行 解决这个问题吗?
这是一个包含 100,000 个响应的最小二乘问题。您的 bigmatrix
是响应(矩阵),beta
是系数(矩阵),而 w1
是残差(矩阵)。
bigmatrix
,以及 w1
,如果明确形成,每个成本
(100,000 * 100,000 * 8) / (1024 ^ 3) = 74.5 GB
这太大了。
由于每个响应的估计都是独立的,因此确实没有必要一次性形成bigmatrix
并尝试将其存储在RAM中。我们可以 一个接一个地拼贴,然后使用一个迭代过程:形成一个拼贴,使用一个拼贴,然后丢弃它。例如,下面考虑了一个尺寸为 100,000 * 2,000
的图块,内存大小为:
(100,000 * 2,000 * 8) / (1024 ^ 3) = 1.5 GB
通过这样的迭代过程,内存的使用得到了有效的控制。
x1 <- rnorm(100000)
x2 <- rnorm(100000)
g <- cbind(x1, x2, x1^2, x2^2)
gg <- crossprod(g) ## don't use `t(g) %*% g`
## we also don't explicitly form `gg` inverse
## initialize `beta` matrix (4 coefficients for each of 100,000 responses)
beta <- matrix(0, 4, 100000)
## we split 100,000 columns into 50 tiles, each with 2000 columns
for (i in 1:50) {
start <- 2000 * (i-1) + 1 ## chunk start
end <- 2000 * i ## chunk end
bigmatrix <- outer(x1, x2[start:end], "<=")
Gw <- crossprod(g, bigmatrix) ## don't use `t(g) %*% bigmatrix`
beta[, start:end] <- solve(gg, Gw)
}
注意,不要尝试计算残差矩阵 w1
,因为它将花费 74.5 GB。如果后期工作中需要残差矩阵,还是尽量把它分解成tiles,一个一个工作
你不需要担心这里的循环。每次迭代中的计算成本足以分摊循环开销。
考虑以下 R 中的示例:
x1 <- rnorm(100000)
x2 <- rnorm(100000)
g <- cbind(x1, x2, x1^2, x2^2)
gg <- t(g) %*% g
gginv <- solve(gg)
bigmatrix <- outer(x1, x2, "<=")
Gw <- t(g) %*% bigmatrix
beta <- gginv %*% Gw
w1 <- bigmatrix - g %*% beta
如果我尝试在我的计算机中 运行 这样的东西,它会抛出内存错误(因为 bigmatrix
太大了)。
你知道我怎样才能达到同样的效果,而不 运行 解决这个问题吗?
这是一个包含 100,000 个响应的最小二乘问题。您的 bigmatrix
是响应(矩阵),beta
是系数(矩阵),而 w1
是残差(矩阵)。
bigmatrix
,以及 w1
,如果明确形成,每个成本
(100,000 * 100,000 * 8) / (1024 ^ 3) = 74.5 GB
这太大了。
由于每个响应的估计都是独立的,因此确实没有必要一次性形成bigmatrix
并尝试将其存储在RAM中。我们可以 一个接一个地拼贴,然后使用一个迭代过程:形成一个拼贴,使用一个拼贴,然后丢弃它。例如,下面考虑了一个尺寸为 100,000 * 2,000
的图块,内存大小为:
(100,000 * 2,000 * 8) / (1024 ^ 3) = 1.5 GB
通过这样的迭代过程,内存的使用得到了有效的控制。
x1 <- rnorm(100000)
x2 <- rnorm(100000)
g <- cbind(x1, x2, x1^2, x2^2)
gg <- crossprod(g) ## don't use `t(g) %*% g`
## we also don't explicitly form `gg` inverse
## initialize `beta` matrix (4 coefficients for each of 100,000 responses)
beta <- matrix(0, 4, 100000)
## we split 100,000 columns into 50 tiles, each with 2000 columns
for (i in 1:50) {
start <- 2000 * (i-1) + 1 ## chunk start
end <- 2000 * i ## chunk end
bigmatrix <- outer(x1, x2[start:end], "<=")
Gw <- crossprod(g, bigmatrix) ## don't use `t(g) %*% bigmatrix`
beta[, start:end] <- solve(gg, Gw)
}
注意,不要尝试计算残差矩阵 w1
,因为它将花费 74.5 GB。如果后期工作中需要残差矩阵,还是尽量把它分解成tiles,一个一个工作
你不需要担心这里的循环。每次迭代中的计算成本足以分摊循环开销。