Inverse of the matrix with several indices

Suppose I have a square matrix $$A_{mn}$$ where $$m=(i,j)$$ and $$n=(k,l)$$. For example, if $$i,j,k,l=1..N$$ the matrix is $$N^2\times N^2$$. For simplicity, let me generate it by

n = 3;
A = RandomReal[1, {n, n, n, n}];


I want to find an inverse of this matrix: $$A^{-1}_{(i,j)(k,l)}$$. I understand how to it by introducing another square matrix $$B=A_{mn}$$ and inverting it:

B = ConstantArray[0, {n^2, n^2}];
Do[
B[[i + n ( j - 1)]][[k + n ( l - 1)]] = A[[i]][[j]][[k]][[l]];
, {i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, n}
]
Bi = Inverse[B];
Ai = Table[
Bi[[i + n ( j - 1)]][[k + n ( l - 1)]], {i, 1, n}, {j, 1, n}, {k,
1, n}, {l, 1, n}];


I believe there should be a better way of doing it(maybe using some tensor commands)?

1 Answer

B = ArrayFlatten[Transpose[A, {3, 1, 4, 2}]];
Ai = Transpose[ArrayReshape[Inverse[B], {n, n, n, n}], {2, 1, 4, 3}];

• Thanks, it works! Is there an easy way to understand why you transpose these layers? Apr 22 at 9:24
• Well, I just do it to get the ordering right. Please don't ask to figure out the right permutation! I just tried a couple of permutations until it worked! ;) Apr 22 at 9:26