矩阵的对称积

Question

Python：是否有任何[[function F这样的

输入：

两个矩阵A和B，两个数字

输出：

矩阵的对称积。例如：

F(A,B,1,1)=AB+BA

F(A,B,2,1)=A^2B+ABA+BA^2（2表示两个A矩阵，1表示一个B矩阵），依此类推。 
该功能对于以下问题是必需的：
计算各种m和n（0 <= m <= L，0 <= n <= L）的矩阵C
(a*A+b*B)^{m+n}=..+a^m b^n C+..

A和B是相当大的矩阵。

Answer 1

我们可以分为两个部分来重写此问题：

生成操作序列。
• 例如F(A, B, 2, 1) -> AAB+ABA+BAA（隐式矩阵乘法）

（有效地）计算这些操作。我们可以注意到，总体上将进行多次计算。例如，当我们计算AAB+ABA+BAA时，可以将所有AA乘法分组并存储结果供以后需要时使用。

要生成序列，我们可以使用more_itertoolsdistinct_permutationsfunction。通过将A编码为0并将B编码为1，它会生成我们想要的计算序列。

要执行一个序列的计算，我们应该利用以前的计算。我们可以使用记忆来记住以前的结果，并且只能执行一次。

memo = {(0,): A, 
        (1,): B
       }  # should be initialized everytimes A, or B changes.

def matmul_perm(A, B, perm):
    if perm in memo:  # If previously computed, return result
        return memo[perm]

    mid = len(perm) // 2

    memo[perm] = (matmul_perm(A, B, perm[:mid]) @ 
                  matmul_perm(A, B, perm[mid:]))  # Split computation in 2 equal part and store result
    return memo[perm]

现在我们可以定义我们的功能：

def F(A, B, na, nb):
    s = 0
    for perm in distinct_permutations((0,)*na + (1,)*nb):
        s += matmul_perm(A, B, perm)
    return s

最后测试我们的程序：

A = np.random.randn(50, 50)

B = np.random.randn(50, 50)

memo = {(1,): B, (0,): A}

np.max(np.abs(F(A, B, 2, 2) - (A@A@B@B + A@B@A@B + B@A@A@B + A@B@B@A + B@A@B@A + B@B@A@A)))
>>> 1.8189894035458565e-12

[在我的计算机上，用大小为F(A, B, 10, 10)的A和B和新的备忘录计算50×50大约需要6秒钟。 （不到一秒钟即可重新计算）