What is $\frac{\partial\left( F^{T}F\right)}{\partial F}$ in tensor notation? 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.)
 A: First, note the two conventions:

*

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

*$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:

*

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



*$_{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'}$



*$_{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()

