$\frac{\partial\left( F^{T}F\right)}{\partial F}$ in tensor notation

Question

In index notation, it my calculations are correct the result should be $$\left(\frac{\partial \left(F^{T}F\right)}{\partial F}\right)_{ijkl} = \frac{\partial\left( F_{mi}F_{mj}\right)}{\partial F_{kl}} = F_{mi}\frac{\partial F_{mj}}{\partial F_{kl}} + F_{mj}\frac{\partial F_{mi}}{\partial F_{kl}} = F_{mi}\delta_{mk}\delta_{jl}+F_{mj}\delta_{mk}\delta_{il} = F_{ki}\delta_{jl}+F_{kj}\delta_{il}.$$

I need to add this to some FEM code that uses tensor notation all the way through, so if possible I'd like to have this in tensor notation too, but I have no idea if, or how, this could be written out. It doesn't look like any combination of index notation that I could find in the definitions.

(Just to be clear, by "tensor notation" I mean things like $F\otimes I^T+F^T\otimes I$, for example.)

Answer 1

First, note the two conventions:

1. $A⊗B = (A_{ik}B_{jl})_{ij,kl} ↭ AXB^𝖳 = (A⊗B)⋅X ↭ \frac{𝐝 AXB^𝖳}{𝐝 X}=A⊗B$
2. $A⊗B = (A_{jl}B_{ik})_{ij,kl} ↭ AXB^𝖳 = (B⊗A)⋅X ↭ \frac{𝐝 AXB^𝖳}{𝐝 X}=B⊗A$

For example, matrixcalulus.org uses (2). I'll be using (1). So let $F$ be $m×n$, then:

$$\begin{aligned} \frac{𝐝 F'F}{𝐝F} &= \frac{\partial F'F}{\partial (F', F)} ⋅ \frac{\partial(F', F)}{\partial F} \\&= \begin{bmatrix} 𝕀_n⊗F' & F'⊗𝕀_n \end{bmatrix}⋅\begin{bmatrix} 𝕋_{m, n} \\ 𝕀_{m, n} \end{bmatrix} \\&= (𝕀_n⊗F')⋅𝕋_{m, n} + (F'⊗𝕀_n)⋅𝕀_{m, n} \\&= 𝕋_{n, n}⋅(F'⊗𝕀_n) + 𝕀_{n, n}⋅(F'⊗𝕀_n) \\&= (𝕋_{n, n} + 𝕀_{n, n})⋅(F'⊗𝕀_n) \end{aligned}$$

In particular the directional (Gâteaux-) derivative is given by:

$$ 𝐃f(F)⋅H = (𝕋_{n, n} + 𝕀_{n, n})⋅(F'⊗𝕀_n)⋅H = (𝕋_{n, n} + 𝕀_{n, n})⋅F'H = H'F + F'H $$

Which agrees with the direct way of computing it via $\frac{𝖽 f(F+εH)}{𝖽ε}\big|_{ε=0}$.

But let me explain the details step by step:

1. In this context "$⋅$" generally does not mean matrix multiplication, but appropriate tensor contraction. For the 4D tensors involved here, typically $A⋅B = (∑_{kl} A_{ij, kl}B_{kl, mn})_{ij, mn}$ which is really just regular old matrix multiplication but with multi-indices.
   • In particular, $(A⊗B)⋅(C⊗D)=(AC⊗BD)$ when dimensions match.
2. $𝕀_{m, n} = (δ_{ik}δ_{jl})_{ij, kl} = 𝕀_m ⊗ 𝕀_n$ is the identity tensor of shape $(m×n, m×n)$.
   • If $A$ is $m×n$ and $B$ is $m'×n'$ then $𝕀_{m, m'}·(A ⊗ B) = A ⊗ B = (A ⊗ B)·𝕀_{n, n'}$
3. $𝕋_{m, n} = (δ_{il}δ_{kj})_{ij, kl}$ is the transpose tensor of shape $(n×m, m×n)$.
   • It cannot be written as a pure tensor of the form $A⊗B$.
   • This follows from $A^𝖳 = ∑_{mn} (e_me_n^𝖳) A (e_me_n^𝖳) = ∑_{mn} (E_{mn} ⊗ E_{nm})⋅A = 𝕋_{m, n}⋅A$
   • I.e. $∑_{mn} (E_{mn} ⊗ E_{nm}) = \sum_{mn}(δ_{im}δ_{kn}δ_{nj}δ_{ml})_{ij, kl} = (δ_{il}δ_{kj})_{ij, kl} = 𝕋_{m, n}$
   • It satisfies $𝕋_{m, n}^𝖳 = 𝕋_{n, m}$, where $(A⊗B)^𝖳 = (A^𝖳⊗B^𝖳)$ is the transpose (for tensors)
   • If $A$ is $m×n$ and $B$ is $m'×n'$ then $𝕋_{m, m'}·(A ⊗ B) = (B⊗A)·𝕋_{n, n'}$

---

Source Code Demo

In python, using https://github.com/google/jax for automatic jacobian computation.

import jax
import numpy as np

def otimes(A, B):
    """Tensor-product A ⊗ B = (Aᵢₖ·Bⱼₗ)ᵢⱼ,ₖₗ"""
    assert A.ndim==2 and B.ndim==2
    return np.einsum('ij,kl -> ikjl', A, B)



def II(m,n):
    """Identity tensor (𝕀ₘ,ₙ)ᵢⱼ,ₖₗ = (δᵢₖ·δⱼₗ)ᵢⱼ,ₖₗ = (𝕀ₘ ⊗ 𝕀ₙ)ᵢⱼ,ₖₗ"""
    I = np.zeros( (m, n, m, n) )
    for i, j, k, l in np.ndindex(I.shape):
        if i==k and j==l:
            I[i,j,k,l] = 1
    return I

def TT(m,n):
    """Transpose tensor (𝕋ₘ,ₙ)ᵢⱼₖₗ = (δᵢₗ·δₖⱼ)ᵢⱼₖₗ"""
    T = np.zeros( (n, m, m, n) )
    for i, j, k, l in np.ndindex(T.shape):
        if i==l and j==k:
            T[i,j,k,l] = 1
    return T

def f(X):
    return X.T @ X

g = jax.jacfwd(f)

def g_manual_intermediate(X):
    l = np.tensordot(otimes(np.eye(n), X.T), TT(m,n))
    r = np.tensordot(otimes(X.T, np.eye(n)), II(m,n))
    return l+r

def g_manual(X):
    m, n = X.shape    
    return np.tensordot(TT(n,n) + II(n,n), otimes(X.T, np.eye(n)))

m,n = 5,4
X = np.random.randn(5,4)

assert (g(X) == g_manual_intermediate(X)).all()
assert (g(X) == g_manual(X)).all()