The Inverse of a Fourth Order Tensor

Question

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.

Answer 1

This is a few years late, but I've been workin gwith these types of tensors for a while now, so better late than never, eh? The OP mentioned elasticity, this post answers the question in the case of three dimensions. Since we deal with tensors that have very special symmetry properties, we need not resort to such a gigantic unfolding as in Miller's response - there will be far fewer than $d^4$ components to collect.

Call the space of all rank-4 tensors with $(12)$,$(34)$, and $(13)(24)$ symmetries SYM. The first two symmetries reduce the number of independent components to $_dC_2$ for each pair of indices, which we organize in a square matrix of size $\frac{d(d-1)}{2}\times\frac{d(d-1)}{2}$. The $(13)(24)$ symmetry implies this matrix is symmetric, leaving us with $\frac{d(d+1)}{2}$ independent components. If we can find a matrix representation $\mathcal D: \text{SYM}\rightarrow GL(\frac{d(d+1)}{2})$ that preserves multiplication across double dot products i.e. $\mathcal D(\boldsymbol{A})=\mathcal D(\boldsymbol{B}:\boldsymbol{C})=\mathcal D(\boldsymbol{B})\cdot\mathcal D(\boldsymbol{C})$, then we'll be able to solve the inverse problem by just mapping it to a place where we know the answer.

Fortunately for us in $d=3$ Mandel form has been developed for just this. We map our tensors to matrices via

$$ \mathcal{D}(A)= \begin{bmatrix} A_{1111} & A_{1122} & A_{1133} & \sqrt{2}A_{1123} & \sqrt{2}A_{1131} & \sqrt{2}A_{1112}\\ & A_{2222} & A_{2233} & \sqrt{2}A_{2223} & \sqrt{2}A_{2231} & \sqrt{2}A_{2212}\\ & & A_{3333} & \sqrt{2}A_{3323} & \sqrt{2}A_{3331} & \sqrt{2}A_{3312}\\ & & & 2A_{2323} & 2A_{2331} & 2A_{2312}\\ & & & & 2A_{3131} & 2A_{3112}\\ & & & & & 2A_{1212}\\ \end{bmatrix}. $$

We'll need to introduce the symmetrization tensor $\mathrm{S}\in\text{SYM}$ which acts to symmetrize the indices it double contracts with, $\mathrm{S}:\hat{e}_i\otimes\hat{e}_j=\frac{1}{2}\left(\hat{e}_i\otimes\hat{e}_j+\hat{e}_j\otimes\hat{e}_i\right)$. It's Mandel form is the identity matrix, $\mathcal{D}(\mathrm{S})=\mathrm{1}_6$, and acts trivially on all $A\in$ SYM, $S:A=A=A:S$. Then we have

$$ A:A^{-1} = \mathbf{I} $$

Symmetrizing we get

$$ S:A:A^{-1} = A:A^{-1} = S. $$

We did this to ensure all objects are in SYM so that we can take the Mandel representation

$$ \mathcal{D}(A:A^{-1}) = \mathcal{D}(A)\cdot\mathcal{D}(A^{-1})=1_6\\ \Downarrow\\ \mathcal{D}(A^{-1}) = \mathcal{D}(A)^{-1} $$

Let's see how this works in practice for the isotropic elasticity tensor, $$ \mathcal D(\mathbf{E})= \begin{bmatrix} \kappa+\frac{4}{3}\mu & \kappa-\frac{2}{3}\mu & \kappa-\frac{2}{3}\mu & 0 & 0 & 0\\ \kappa-\frac{2}{3}\mu & \kappa+\frac{4}{3}\mu & \kappa-\frac{2}{3}\mu & 0 & 0 & 0\\ \kappa-\frac{2}{3}\mu & \kappa-\frac{2}{3}\mu & \kappa+\frac{4}{3}\mu & 0 & 0 & 0\\ 0 & 0 & 0 & 2\mu & 0 & 0\\ 0 & 0 & 0 & 0 & 2\mu & 0\\ 0 & 0 & 0 & 0 & 0 & 2\mu \end{bmatrix} $$ where $\kappa$ is the bulk modulus and $\mu$ the shear modulus. We take the inverse of the matrix $$ \mathcal{D}(\mathbf{E}^{-1})=\mathcal{D}(\mathbf{E})^{-1}= \begin{bmatrix} \frac{1}{9} \left(\frac{1}{\kappa }+\frac{3}{\mu }\right) & \frac{1}{9}\left(\frac{1}{ \kappa }-\frac{3}{2 \mu }\right) & \frac{1}{9}\left(\frac{1}{\kappa }-\frac{3}{2 \mu }\right) & 0 & 0 & 0 \\ \frac{1}{9}\left(\frac{1}{ \kappa }-\frac{3}{2 \mu }\right) & \frac{1}{9} \left(\frac{1}{\kappa }+\frac{3}{\mu }\right) & \frac{1}{9}\left(\frac{1}{ \kappa }-\frac{3}{2 \mu }\right) & 0 & 0 & 0 \\ \frac{1}{9}\left(\frac{1}{ \kappa }-\frac{3}{2 \mu }\right) & \frac{1}{9}\left(\frac{1}{ \kappa }-\frac{3}{2 \mu }\right) & \frac{1}{9} \left(\frac{1}{\kappa }+\frac{3}{\mu }\right) & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{2 \mu } & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{2 \mu } & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{2 \mu } \\ \end{bmatrix} $$ from which we can read off the rank-4 tensor components.

Answer 2

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.

Answer 3

$ \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)