r - 通过预计算加速矩阵乘积

Question

我正在处理一个矩阵必须多次迭代计算的问题。必要的矩阵乘法采用以下形式

t(X) %*% ( ( X %*% W %*% t(X) * mu0 ) * mu1 )

其中 X 是 N x P 而 W 是一个对称 P x P矩阵。 mu0 和 mu1 是 N x 1 向量，计算成本低，并按元素输入各自的乘积。

不幸的是，N 可能非常大，这会导致巨大的计算需求，因为 X %*% W %*% t(X) 是 N×N。我想知道是否有任何策略或计算技巧，例如基于矩阵分解，可以用来加速这里的计算。在每次迭代中，mu0 和 mu1 都会发生变化，但是 X 和 W 是固定的，因此包括这些矩阵在内的任何预计算将工作。

基准测试

到目前为止，我能想到的最快方法是进行一些明显的预计算:

# fake data
N   = 2500
P   = 10
X   = matrix(rnorm(N*P), N, P)
W   = matrix(rnorm(P*P), P, P)
mu0 = rnorm(N)
mu1 = rnorm(N)

# precomputations
tX  = t(X)
XWX = X %*% W %*% t(X)

# functions
f_raw      = function(X, W, mu0, mu1){t(X) %*% ( ( X %*% W %*% t(X) * mu0 ) * mu1 )}
f_precomp  = function(XWX, tX, mu0, mu1){tX %*% ( ( XWX * mu0 ) * mu1 )}


# benchmark
microbenchmark::microbenchmark(f_raw(X, W, mu0, mu1), 
                               f_precomp(XWX, tX, mu0, mu1))


Unit: milliseconds
                         expr      min       lq     mean   median       uq      max neval
        f_raw(X, W, mu0, mu1) 283.5918 286.5080 299.4621 289.5151 302.9726 355.4271   100
 f_precomp(XWX, tX, mu0, mu1) 167.4169 168.7336 180.6468 171.0852 197.7475 263.8090   100

Answer 1

您的示例数据和尝试:

# fake data
N   = 2500
P   = 10
X   = matrix(rnorm(N*P), N, P)
W   = matrix(rnorm(P*P), P, P)
mu0 = rnorm(N)
mu1 = rnorm(N)

# precomputations
tX  = t(X)
XWX = X %*% W %*% t(X)

# functions
f_raw      = function(X, W, mu0, mu1){t(X) %*% ( ( X %*% W %*% t(X) * mu0 ) * mu1 )}
f_precomp  = function(XWX, tX, mu0, mu1){tX %*% ( ( XWX * mu0 ) * mu1 )}

更好的办法是预先计算

WX <- tcrossprod(W, X)

然后

f_better <- function (X, mu0, mu1, WX) crossprod(X, (mu0 * mu1) * X) %*% WX

在我的笔记本电脑上进行基准测试:

microbenchmark::microbenchmark(f_raw(X, W, mu0, mu1), 
                               f_precomp(XWX, tX, mu0, mu1),
                               f_better(X, mu0, mu1, WX))
#Unit: milliseconds
#                         expr        min         lq      mean     median
#        f_raw(X, W, mu0, mu1) 236.926979 238.984573 243.86319 242.212716
# f_precomp(XWX, tX, mu0, mu1) 151.190689 152.612059 156.25277 155.646093
#    f_better(X, mu0, mu1, WX)   1.113031   1.126434   1.17974   1.138207
#         uq        max neval
# 244.721163 270.477737   100
# 157.970378 182.935082   100
#   1.146352   5.130876   100

验证正确性:

ans_raw <- f_raw(X, W, mu0, mu1)
ans_precomp <- f_precomp(XWX, tX, mu0, mu1)
ans_better <- f_better(X, mu0, mu1, WX)

all.equal(ans_raw, ans_precomp)
# [1] TRUE
all.equal(ans_raw, ans_better)
# [1] TRUE

不知道mu0和mu1在你的实际应用中是否保证是非负的。如果是这样，则以下速度快 2 倍:

f_fastest <- function (X, mu0, mu1, WX) crossprod(sqrt(mu0 * mu1) * X) %*% WX