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

• I would try to transform the problem to an equivalent problem stated in terms of matrices. Then you can use existing codes for inverting matrices and transform back. Since a 4th order tensor corresponds to a block matrix, it seems like this approach may work. Can you give more detail regarding the tensor-tensor product you are using to define the inverse. Commented Jan 24, 2016 at 15:18
• @K.Miller: Thanks for the attention. I added the definition to the question. :) Commented Jan 24, 2016 at 15:35
• I assume by double scalar product you mean the sum over indexes $m,n$? Commented Jan 24, 2016 at 16:17
• @K.Miller: Yes. :) Whatever terminology you like. :) Commented Jan 24, 2016 at 16:32
• There are three $4^{th}$-order isotropic tensors, but only one of them $(E_{ijkl}=\delta_{ik}\delta_{jl})$ acts like an identity with respect to the double contraction product \eqalign{E:A &= A:E = A} The other two $(F_{ijkl}=\delta_{il}\delta_{jk},\,G_{ijkl}=\delta_{ij}\delta_{kl})$, produce a trace or a transpose under the product\eqalign{F:A &= A:F = A^T \cr G:A &= A:G = {\rm tr}(A)\,I \cr}
– lynn
Commented Jan 24, 2016 at 20:24

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.

• I just noticed that if the tensor ${\bf{A}}$ has the symmetries $A_{ijkl}=A_{jikl}=A_{ijlk}=A_{klij}$ then the unfolded matrix form $A$ will be singular! :) Commented Mar 5, 2016 at 14:43

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.

• (+1) Thanks for your answer Damian. :) Commented Aug 3, 2022 at 12:46

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

• (+1): Thanks for your contribution Iynn. I think the crucial point here is that how properties of the tensor transfer to the reshaped form. Two of the important ones are the tensor product and tensor inverse as you mentioned. Commented Dec 2, 2022 at 9:16
• @HoseinRahnama Agreed. The critical part is how the reshape() (or as K Miller calls it unfolding) operation works. I prefer using $a={\rm vec}(A)={\rm vec}({\cal A})$ because the indices on the vector $a$ are always ordered lexicographically, reflecting the fact that in most computer languages matrices (and multidimensional arrays) are stored as a contiguous linear array with special indexing rules. This means that the vec() and reshape() operations have essentially zero cost.
– lynn
Commented Dec 3, 2022 at 17:38