我花了点时间才弄明白你在做什么，
因为z在range(3)上迭代，所以aal依次是a[:,0,0]、a[:,0,1]、a[:,0,2]。
或者一次完成：

In [178]: aaL = a[:,0,:]; aaL
Out[178]: 
array([[1, 2, 1],
       [1, 2, 1],
       [1, 2, 1]])

aam也做同样的迭代，所以用matmul/@/dot计算它们的乘积之和为：

In [179]: aaL.T@aaL
Out[179]: 
array([[ 3,  6,  3],
       [ 6, 12,  6],
       [ 3,  6,  3]])

或者在einsum中：

In [180]: np.einsum('ji,jk->ik',aaL,aaL)
Out[180]: 
array([[ 3,  6,  3],
       [ 6, 12,  6],
       [ 3,  6,  3]])

您的索引数组：

In [183]: bb
Out[183]: [(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2), (2, 0), (2, 1), (2, 2)
In [185]: np.array(bb)[:3,:]
Out[185]: 
array([[0, 0],
       [0, 1],
       [0, 2]])

如果我把它推广到bb的其余范围：

In [192]: for i in range(3):
     ...:     aaL = a[:,i]
     ...:     print(aaL.T@aaL)
     ...:     
[[ 3  6  3]
 [ 6 12  6]
 [ 3  6  3]]
[[27 36 18]
 [36 48 24]
 [18 24 12]]
[[ 75  90  45]
 [ 90 108  54]
 [ 45  54  27]]

向einsum添加尺寸：

In [195]: np.einsum('jmi,jmk->mik', a,a)
Out[195]: 
array([[[  3,   6,   3],
        [  6,  12,   6],
        [  3,   6,   3]],

       [[ 27,  36,  18],
        [ 36,  48,  24],
        [ 18,  24,  12]],

       [[ 75,  90,  45],
        [ 90, 108,  54],
        [ 45,  54,  27]]])