矩阵的对称积

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的矩阵C (0 <=m<=L, 0<=n<=L)

(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_itertools distinct_permutations function。通过将 A 编码为 0 并将 B 编码为 1，它生成了我们想要的计算序列。

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

# /!\ IMPORTANT: This initialization is wrong in some case, see the EDIT.
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

在我的电脑上用 A 和 B 大小 50×50 计算 F(A, B, 10, 10) 大约需要 6 秒，还有一个新的备忘录。（第二次重新计算不到一秒钟）

编辑

此实现在使用 na=nb=0 调用时永远递归。最简单的修复方法是将 memo 初始化更改为

memo = {(0,): A, 
        (1,): B,
        (): np.eye(*A.shape)  # Empty product !
       }

矩阵的对称积

symmetric product of matrices

python

matrix-multiplication