通过预计算加速矩阵乘积

Question

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

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

其中X和是对称矩阵N x P。和是计算成本低廉的向量，并按元素输入相应的乘积。W P x Pmu0mu1N x 1

不幸的是，N可能非常大，这导致了巨大的计算X %*% W %*% t(X)需求N x 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