大矩阵的Kronecker乘积

Question

我正在寻找一种高效的方法来计算两个大矩阵的Kronecker积。我已经尝试使用方法kronecker()，如下所示：

 I = diag(700)
 data = replicate(15, rnorm(120))
 test = kronecker(I,data)

然而，执行需要很长时间，然后会出现以下错误：

 Error: cannot allocate vector of size 6.8 Gb

- Mayou

---

你正在尝试创建一个近十亿个元素的矩阵。这不是“效率”的问题 - 尽管如果你知道其中一个矩阵是对角线矩阵，Kronecker积将变得非常稀疏。也许可以一部分一部分地完成它？Kronecker积可以很容易地“细分”。 - Carl Witthoft

可能应该将上面的评论发布为答案，因为除非您拥有一台具有大量内存的计算机，否则您无法执行此操作。 - ppaulojr

这可能会给你一个思路：https://dev59.com/GnPYa4cB1Zd3GeqPfiq_ - ppaulojr

你真的需要矩阵（系数数组），还是只需要能够对其进行一些计算（乘以它，解线性系统，计算特征值等）？如果你确实需要矩阵，由于它是块对角线形式的，你可以使用Matrix包中的bdiag函数仅存储非零系数。 - Vincent Zoonekynd

实际上，我需要它作为计算的一部分。但是我喜欢您的bdiag想法。谢谢。 - Mayou

Answer 1

只要您使用Matrix::Diagonal来构建对角矩阵，您的test对象就会自动构建为稀疏矩阵。

library(Matrix)
I=Diagonal(700)
data = replicate(15,rnorm(120))
system.time(test <- kronecker(I,data))
##   user  system elapsed
##  0.600   0.044   0.671 
dim(test)
## [1] 84000 10500
format(object.size(test),"Mb")
## [1] "19.2 Mb"

- Ben Bolker

Answer 2

如果你需要计算 kron(I,A)*v，其中 v 是一个向量，你可以使用 vec(A*V) 来完成。这里的 V 将 v 重塑为一个矩阵。这个方法使用了更一般的规则：vec(ABC)=kron(C',A)*vec(B)。这避免了形成 Kronecker 乘积，并且使用更少的操作来执行计算。

请注意，根据矩阵存储方式（列与行），V 可能需要转置。

- Paul Fackler