The Inverse of a Fourth Order Tensor Suppose that we have a fourth order tensor ${\bf{A}}$
$${\bf{A}}=A_{ijkl} {\bf{e}}_i \otimes {\bf{e}}_j \otimes {\bf{e}}_k \otimes {\bf{e}}_l$$
in the orthonormal basis $\{{\bf{e}}_1,{\bf{e}}_2,{\bf{e}}_3\}$ for $\mathbb{R}^3$. Here, we have used Einstein summation convention which assume sum over any repeated index. Then we define the inverse of ${\bf{A}}$ denoted by ${\bf{B}}$ as follows
$${\bf{A}} : {\bf{B}} = {\bf{B}} : {\bf{A}} = {\bf{I}}$$
where ${\bf{I}}$ is the fourth order identity tensor
$$\begin{align}
{\bf{I}} &=I_{ijkl} {\bf{e}}_i \otimes {\bf{e}}_j \otimes {\bf{e}}_k \otimes {\bf{e}}_l\\
&=\delta_{ik} \delta_{jl} {\bf{e}}_i \otimes {\bf{e}}_j \otimes {\bf{e}}_k \otimes {\bf{e}}_l
\end{align}$$
where $\delta_{ij}$ is the Kronecker's Delta
$$\delta_{ij}=
\begin{cases}
1 & i=j \\
0 & i \ne j
\end{cases}$$
and $:$ is the double contraction defined by
$${\bf{A}} : {\bf{B}}=A_{ijmn}B_{mnkl} {\bf{e}}_i \otimes {\bf{e}}_j \otimes {\bf{e}}_k \otimes {\bf{e}}_l$$
I want to write a code to compute ${\bf{B}}$. However, I could not find any good resource on the net that gives the elements of ${\bf{B}}$ in terms of the elements of ${\bf{A}}$. How should I compute ${\bf{B}} = {\bf{A}}^{-1}$?
An application of this can be found in the theory of elasticity, where the fourth order tensor have the following symmetries
$$A_{ijkl}=A_{jikl}=A_{ijlk}=A_{klij}.$$
That would be great if you can also touch upon this.
 A: In general
$$
\sum_{m=1}^N\sum_{n=1}^N A_{ijmn}B_{mnkl} = \delta_{ik}\delta_{jl} 
$$
for $i,j,k,l = 1,\ldots,N$. Now suppose the tensor $\mathbf{A}$ is unfolded as the $N^2\times N^2$ matrix given by
$$
A = \left[
\begin{array}{ccccccccc}
A_{1111} & A_{1112} & \cdots & A_{111N} & \cdots & A_{11N1} & A_{11N2} & \cdots & A_{11NN}\\
A_{1211} & A_{1212} & \cdots & A_{121N} & \cdots & A_{12N1} & A_{12N2} & \cdots & A_{12NN}\\
\vdots & \vdots & &\vdots & & \vdots & \vdots & & \vdots\\
A_{NN11} & A_{NN12} & \cdots & A_{NN1N} & \cdots & A_{NNN1} & A_{NNN2} & \cdots & A_{NNNN}\\
\end{array}
\right]
$$
If the same is done for $\mathbf{B}$, then the matrix product $AB$ corresponds to the unfolded tensor product $\mathbf{A} : \mathbf{B}$. Let $E$ be the same unfolding of the tensor identity matrix $\mathbf{I}$. Then 
$$
\mathbf{A} : \mathbf{B} = \mathbf{I} \iff AB = E
$$
Similarly,
$$
\mathbf{B} : \mathbf{A} = \mathbf{I} \iff BA = E
$$
In order for such a matrix $B$ to exist it follows that $AB = BA$ must hold, i.e., the product of $A$ and $B$ must commute. Now suppose $A$ is invertible. Then $B$ must satisfy the equations
$$
B = A^{-1}E \quad\text{and}\quad B = EA^{-1} \iff A^{-1}E = EA^{-1} \iff EA = AE
$$
Thus, in the case that $A$ is invertible, the product of $A$ and $E$ must commute in order for such a matrix $B$ to exist. If these conditions are met, then you should be able to compute a unique $B$ from the equation $AB = E$ via the LU factorization or by some other matrix factorization.
A: $
\def\bbR#1{{\mathbb R}^{#1}}
\def\A{{\cal A}}\def\B{{\cal B}}\def\I{{\cal I}}
\def\vc{{\rm vec}}
\def\qiq{\quad\implies\quad}
$The tensor $\A\in\bbR{3\times 3\times 3\times 3}$ can be turned into a vector $a\in\bbR{81}$ by column stacking, and subsequently unstacked into a matrix $A\in\bbR{9\times 9}$ which can be inverted
$$\eqalign{
a &= [\A_{1111}, \;\A_{1112}, \;\A_{1113}, \;\A_{1121}, \;\A_{1122}, \;\A_{1123}, \;\A_{1131}, \,\cdots,\,\A_{3332},\;\A_{3333}]^T \\
A &= {\rm reshape}(a, 9,9)
}$$
And don't forget that a double dot product becomes a single dot product when switching between the tensor and the matrix representation, i.e.
$$\I=\A:\B \quad\iff\quad I=A\cdot B$$
As for coding, assuming that you have the $81$ values loaded into a tensor variable Ta (or a matrix variable A or a vector a) then you can code the whole problem in one line in Julia or Matlab
Tb = reshape(inv(reshape(Ta,9,9)), 3,3,3,3)

but you should probably break it down into the more readable
A = reshape(Ta, 9,9)
B = inv(A)
Tb = reshape(B, 3,3,3,3)

