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$. 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.

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.