I know this a bit of an older question, but I came across it while trying to understand the same thing, so let me try to explain how I came to think of this, which is not with tensors first, but maps and operators.
Short answer
I can't quite follow what you're suggesting in terms of tensor products, but in general I think that the calculation of the backwards pass will depend on the details of what your function is and may requires some custom derivation. You may be doing this correctly, but I'm not sure that a tensor product or matrix multiplication is always the most straightforward interpretation/implementation.
The answers you link to give a very practical point of view on how to perform reverse-mode AD, but there's a more abstract way of looking at it as well, in terms of vector-Jacobian products. I think the JAX documentation has a good explanation of this. Let me explain that on a vector-valued function before getting to matrix-valued functions
Vector-valued function
Let's say we have a function $f: \mathbb{R}^n \rightarrow \mathbb{R}^m$. The Jacobian (or tangent linear map) of this function evaluated at some "primal" point $x_0 \in \mathbb{R}^n$ is a linear map from the input space to the output space, so $\mathcal{J}: \mathbb{R}^n \rightarrow \mathbb{R}^m$. Since this is a linear operator, we could represent it numerically using a matrix $J \in \mathbb{R}^{m \times n}$.
Reverse-mode AD is really an algorithm for computing the "vector-transpose Jacobian product (VJP)" $\mathbf{v}^T J$ for some vector $\mathbf{v} \in \mathbb{R}^m$ (I think this is what you're calling the "adjoints of the parents"). If $f$ consists of a computational graph over other functions, this happens by propagating the "seed" vector $\mathbf{v}$ backwards from the final nodes. In any case, all we have to do to implement a new function is the VJP for that function, so we can ignore anything else that's happening upstream and downstream, assuming that $\mathbf{v}$ will be handed to us. For instance, if we wanted to construct the full Jacobian, we could repeat the backwards pass $m$ times using unit basis vectors, each time computing one row of $J$.
To be a little more precise, while $\mathbf{v}$ has the dimensions of an element of the output space, it is actually an element of the cotangent space, and what we evaluate is the pullback or adjoint map $J^T \mathbf{v}$. As an easy example, consider $f(x) = A \mathbf{x}$, with $A \in \mathbb{R}^{m \times n}$. Clearly $A$ is the Jacobian of this function, so the pullback $\mathcal{J}: \mathbb{R}^m \rightarrow \mathbb{R}^n$ is evaluated numerically by $A^T \mathbf{v}$.
Matrix-matrix multiplication
For scalar- and vector-valued functions this probably seems irrelevant and abstract, but I think it's helpful for more general cases. For instance, let's look at matrix-matrix multiplication $f(X) = A X$, with $A \in \mathbb{R}^{m \times n}$ and $X \in \mathbb{R}^{n \times \ell}$, so that $f: \mathbb{R}^{n \times \ell} \rightarrow \mathbb{R}^{m \times \ell}$. The tangent map is another map $\mathcal{J}: \mathbb{R}^{n \times \ell} \rightarrow \mathbb{R}^{m \times \ell}$, which can once again be represented numerically by the matrix $A$. The pullback applied to $V \in \mathbb{R}^{m \times \ell}$ is just $A^T V$. You could view this as a specific 4D tensor product, but I think it's much easier to look at it as a map.
What if we instead want to differentiate $g(A) = A X$? We still have $A \in \mathbb{R}^{m \times n}$ and $X \in \mathbb{R}^{n \times \ell}$, but now $g: \mathbb{R}^{m \times n} \rightarrow \mathbb{R}^{m \times \ell}$. This is still a linear function in $A$, but now we evaluate the tangent map by right-multiplying by the matrix $X$:
$$
\mathcal{J}(\cdot) = (\cdot) X
$$
The adjoint map is $\mathcal{J}^\dagger: \mathbb{R}^{m \times \ell} \rightarrow \mathbb{R}^{m \times n}$, which can be evaluated by
$$
\mathcal{J}^\dagger(\cdot) = (\cdot) X^T.
$$
In other words, during the backward pass, when you are given the seed $V \in \mathbb{R}^{m \times \ell}$ you compute $V X^T$.
To make a general matrix-matrix product $F(A, X) = AX$ differentiable with respect to both arguments, you evaluate the adjoint map with respect to each input (basically partial differentiation). Now $F: \mathbb{R}^{m \times n} \times \mathbb{R}^{n \times \ell} \rightarrow \mathbb{R}^{m \times \ell}$, so the tangent maps with respect to each input are $\mathcal{J}_A: \mathbb{R}^{m \times n} \rightarrow \mathbb{R}^{m \times \ell} $ and $\mathcal{J}_X = \mathbb{R}^{n \times \ell} \rightarrow \mathbb{R}^{m \times \ell}$, as before. You can derive the operators in a principled way by computing the total derivative:
$$
dF = (dA) X + A (dX)
$$
This is linear with respect to each input, and the tangent and adjoint maps can be computed as for the single-argument case.
Nonlinear matrix-valued function
So far all the example I gave are linear, so the tangent and adjoint maps are fairly obvious. Let me give a simple nonlinear example as well.
Let's take $f(X) = X^T X$, with $X \in \mathbb{R}^{m \times n}$ and $f: \mathbb{R}^{m \times n} \rightarrow \mathbb{R}^{n \times n}$. The total derivative of $f$ is
$$
df = (dX)^T X + X^T (dX),
$$
so the tangent map $\mathcal{J}: \mathbb{R}^{m \times n} \rightarrow \mathbb{R}^{n \times n}$ can be evaluated with
$$
\mathcal{J}(\cdot) = (\cdot)^T X + X^T (\cdot).
$$
To compute the adjoint we have to do a little more work and use the definition of the adjoint: for a linear operator $\mathcal{A}: \mathcal{W} \rightarrow \mathcal{V}$, the adjoint $\mathcal{A}^\dagger: \mathcal{V} \rightarrow \mathcal{W}$ is defined using the inner product $\langle \cdot, \cdot \rangle$:
$$
\langle \mathcal{A} w, v \rangle = \langle w, \mathcal{A}^\dagger v \rangle, \qquad w \in \mathcal{W}, v \in \mathcal{V}.
$$
For matrices, the standard inner product is the trace: $\langle A, B \rangle = \mathrm{tr}\{ A B^T \}$. In our case, if we choose a arbitrary elements $W$ of the tangent space $\mathbb{R}^{m \times n}$ and $V$ of the cotangent space $\mathbb{R}^{n \times n}$, then
$$
\langle \mathcal{J}(W), V \rangle = \mathrm{tr} \left\{ (W^T X + X^T W) V^T \right\}.
$$
Using some properties of the trace we can show that
$$
\mathrm{tr} \left\{ (W^T X + X^T W) V^T \right\} = \mathrm{tr} \left\{ W (V + V^T) X^T \right\} = \langle W, X (V + V^T) \rangle,
$$
from which it follows that the pullback is $\mathcal{J}^\dagger V = X(V + V^T)$.
The general recipe
- Define what spaces your function maps over, and what the tangent and linear maps will be.
- Compute the total derivative of the function and extract the tangent map.
- Use the definition of the adjoint and an appropriate inner product to derive the adjoint map.
- Implement the pullback as an evaluation of the adjoint map on the "seed".
For a great example of how to derive reverse-mode AD rules for a complicated function I'd recommend checking out the OptNet paper, where they show how to backpropagate through a quadratic program.
