可能不是最有效的，但这里有一个手动解决方案：

M = np.random.randint(0, 10, (2, 2, 2, 2))

def matmatmul(a, b):
    
    output = np.zeros((a.shape[0], b.shape[1]), dtype = object)
    
    for i in range(output.shape[0]):
        for j in range(output.shape[1]):
            row = a[i]
            col = b[:, j]
            
            output[i,j] = sum([r @ c for r,c in zip(row, col)])
            
    return output

def matmatpow(a, n):
    if n == 1:
        return a
    else:
        output = matmatmul(a, a)
        
        for i in range(2, n):
            output = matmatmul(output, a)
        
        return output
    
M2 = matmatpow(M, 2)

print(M2)

[[A, B], [C, D]] = M

assert np.all(M2[0,0] == A@A + B@C)
assert np.all(M2[0,1] == A@B + B@D)
assert np.all(M2[1,0] == C@A + D@C)
assert np.all(M2[1,1] == C@B + D@D)

您只是在计算4x4分块矩阵的正规矩阵积，该分块矩阵是通过连接较小的矩阵A到D创建的。
在MATLAB中，你的期望结果使用一些任意矩阵：

A = [1, 2
     3, 4];
B = [5, 6
     7, 8];
C = [ 9, 10
     11, 12];
D = [13, 14
     15, 16];

res = [A*A + B*C, A*B + B*D
       C*A + D*C, C*B + D*D]

res =
   118   132   174   188
   166   188   254   276
   310   356   494   540
   358   412   574   628

4 × 4分块矩阵及其平方：
一个二个一个一个

定义一组（2，2）数组及其组合：

In [45]: A,B,C,D = [np.arange(i,i+4).reshape(2,2) for i in range(4)]    
In [46]: M=np.array([[A,B],[C,D]])

您需要的M^2阵列：

In [47]: np.array([[A@A + B@C, A@B + B@D],
    ...:       [C@A + D@C, C@B + D@D]])
Out[47]: 
array([[[[12, 16],
         [28, 40]],

        [[16, 20],
         [40, 52]]],

       [[[28, 40],
         [44, 64]],

        [[40, 52],
         [64, 84]]]])

使用einsum也是一样，其中j和l是积和维：

In [48]: np.einsum('ijkl,jmln->imkn',M,M)
Out[48]: 
array([[[[12, 16],
         [28, 40]],

        [[16, 20],
         [40, 52]]],

       [[[28, 40],
         [44, 64]],

        [[40, 52],
         [64, 84]]]])

matmul等价于'ijkl，ijlm-〉ijkm'，其中ij是batch的维数，l是乘积和。通常，einsum可以通过一些重新整形和广义转置来重现。但我将把这个留给其他人来探索。
通过使用einsum索引、转置和整形数组，我可以得到等效的结果：

In [56]: np.matmul(M.transpose(0,2,1,3).reshape(4,4),M.transpose(0,2,1,3).reshape(4,4))
Out[56]: 
array([[12, 16, 16, 20],
       [28, 40, 40, 52],
       [28, 40, 40, 52],
       [44, 64, 64, 84]])

再多一点的按摩就变成了所期望的（4，4，4，4）

In [57]: np.matmul(M.transpose(0,2,1,3).reshape(4,4),M.transpose(0,2,1,3).reshape(4,4)).reshape(2,2,2,2).transpose(0,2,1,3)
Out[57]: 
array([[[[12, 16],
         [28, 40]],

        [[16, 20],
         [40, 52]]],

       [[[28, 40],
         [44, 64]],

        [[40, 52],
         [64, 84]]]])