Expansions of iterated, or nested, derivatives, or vectors--conjectured matrix computation The entry OEIS A139605 (also related OEIS A145271) has a matrix computation for the partition polynomials that represent the expansions of iterated derivatives, or vectors in differential geometry, 
$$(g(x)D_x)^n.$$
The formula section of A139605 contains the matrix formula. Multiply the $n$-th diagonal (with $n=0$  the main diagonal) of the lower triangular Pascal matrix A007318 by $g_n = D_x^n g(x)$  to obtain the matrix $VP$  with $VP_{n,k} = \binom{n}{k}g_{n-k} $. Then 
$$(g(x)D)^n = (1, 0, 0,..) [VP \dot \; S]^n (1, D, D^2, ..)^T,$$
where S is the shift matrix A129185, representing differentiation in the divided powers basis $x^n/n!$. 
Example:
$$(g(x)D_x)^3$$
$$= (1, 0, 0, 0) [VP \dot \; S]^3 (1, D, D^2, D^3)^T$$
$$= \begin{pmatrix}
1 &  0 &  0 & 0
\end{pmatrix} \begin{pmatrix}
0 & g_0 & 0 & 0 \\ 
0 & g_1 & g_0 & 0\\ 
0 & g_2 & 2g_1 & g_0 \\ 
0 & g_3 & 3g_2 & 3g_1
\end{pmatrix}^3 \begin{pmatrix}
1 \\ 
D \\ 
D^2 \\ 
D^3 \end{pmatrix} $$
$$ = [g_0g_1^2  + g_0^2 g_2] D + 3 g_0^2g_1 D^2 + g_0^3D^3 $$
And, the pdf Mathemagical Forests gives a diagrammatic method for creating forests of trees through "natural growth" that represent the partition polynomials.
I have either lost a proof of the validity of this formula or got sidetracked before I developed one.
Question:
Can someone prove this conjecture?
Some background:
The refined Eulerian numbers (RENs) of A145271 are related analytically to the compositional inversion of functions and formal generating series and to flow fields generated by tangent vectors. The $n$-th row of RENs are the numerical coefficients of the expansion of $(g(x)\frac{d}{dx})^ng(x)$ in terms of the monomials in the derivatives of $g(x)$, i.e.,
$$g_k=\frac{d^k}{dx^k}g(x).$$
For example,
$$(g(x)\frac{d}{dx})^3g(x) = 1 g_0^1 g_1^3 + 4 g_0^2 g_1^1 g_2^1 + 1 g_0^3 g_3^1.$$
With $(\omega,x) = (f(x),f^{(-1)}(\omega))$ and $g(x) = 1/f^{'}(x)$,
$$\exp[t g(x)d/dx]x = \exp[td/d\omega]f^{(-1)}(\omega) = f^{(-1)}(t+\omega)=f^{(-1)}(t+f(x)).$$
Evaluated at the origin of $x$, this gives the compositional inverse 
$$\exp[tg(x)d/dx] x |_{x=0}=f^{(-1)}(t).$$
See also
1) MO-Q Guises of the refined Eulerian numbers generated by tangent vectors
2) MO-Q Important formulas in combinatorics
3) MO-Q Why is there a connection between enumerative geometry and nonlinear waves?
 A: I have finally written up the proof in detail. It is in my note


*

*Darij Grinberg, Commutators, matrices and an identity of Copeland, also available as arXiv:1908.09179v1.


Your result is a particular case of Theorem 4.2. More precisely, you get it from Theorem 4.2 if you set $\mathbb{L}$ to be the ring of differential operators (whatever kind of differential operators you are considering), $\mathbb{K}$ to be the base ring, $a$ to be the differentiation operator (which is your $D$), $x$ to be the "multiplication by $t$" operator (you use $x$ for what I call $t$, but I keep them separate because my ring $\mathbb{L}$ doesn't have to contain the polynomial ring $\mathbb{K}\left[t\right]$), and $h$ to be the identity map.
I derive this from a more general formula (Theorem 2.7), which expresses $\left(ba\right)^n$ as a matrix product when $a$ and $b$ are two arbitrary elements of a noncommutative ring.
I work in a setting that allows the matrices to be finite or infinite as the reader prefers. This is responsible for a lot of the length of the note (as I have to set up an appropriate ring of infinite matrices on which multiplication is defined and associative, but also mess around with partial equalities in the case of finite matrices because the finite version of the shift operator doesn't behave as well as the infinite one). I hope the exposition is not completely wasted; meanwhile I trust you know what to skip and what to skim.
A: [This is a copy of Darij Grinberg's answer in the comments above.]
Nice question. Here is a sketch of a proof: Let me use infinite vectors and matrices instead of your finite ones. For each differential operator $E$, construct the infinite column vector
$$h_E =\begin{pmatrix}
D^0E \\ 
D^1E \\ 
D^2E \\ 
D^3E \\
.. \end{pmatrix}$$
and let
$$U=VP⋅S .$$
Then, your claim is that the first entry of the column vector $U^nh_1$ is $(gD)^n$. Let me claim something stronger:
$$U^nh_1=h_{(gD)^n}.$$
To prove this, it clearly suffices to show that
$$Uh_E=h_{gDE}$$
for each differential operator $E$.
How to prove this identity? Well, for each $i$, the $i$-th entry of $Uh_E$ is given by
$$(Uh_E)_i=\sum_{k=0}^i \binom{i}{k} g_{i−k}D^{k+1}$$
while the $i$-th entry of $h_{gDE}$ is given by
$$(h_{gDE})_i=D^igDE .$$
We need to prove that these two entries are equal. It is clearly enough to show that
$$\sum_{k=0}^i \binom{i}{k}g_{i−k} D^k = D^i g .$$
Now, this can be proven by straightforward induction on $i$, like the binomial formula. I want to say it also follows from the binomial formula, but right now I don't see how (probably a simple exercise in the umbral viewpoint).
A: Here is probably something close to my original thinking on the identity.
Taylor series analysis gives us, with $D_x^k g(x) = \frac{d^k}{dx^k} g(x)= g_k$,
$$g(t+x) = e^{tD_x}g(x)= \sum_{k \geq
 0} g_k\frac{t^k}{k!}.$$
Then 
$$g(x+t)g'(x+t)= g(x+t)D_xe^{tD_x}g(x)$$
$$= \sum_{k \geq
 0} g_k\frac{t^k}{k!}D_x\sum_{k \geq
 0} g_k\frac{t^k}{k!}$$
can be written as a binomial convolution (umbrally, $e^{a.t}e^{b.t}=e^{(a.+b.)t}$) in (singly infinite) matrix form in the divided powers basis  $\frac{t^k}{k!}$ as
$$ \begin{pmatrix}
g_0 & 0 & 0 & 0 \\ 
g_1 & g_0 & 0 & 0\\ 
g_2 & 2g_1 & g_0 & 0\\ 
g_3 & 3g_2 & 3g_1 & g_0
\end{pmatrix} 
\begin{pmatrix}
g_1\\ 
g_2 \\ 
g_3 \\ 
g_4 \end{pmatrix} $$
$$ = \begin{pmatrix}
g_0 & 0 & 0 & 0 \\ 
g_1 & g_0 & 0 & 0\\ 
g_2 & 2g_1 & g_0 & 0\\ 
g_3 & 3g_2 & 3g_1 & g_0
\end{pmatrix} \begin{pmatrix}
0 & 1 & 0 & 0 \\ 
0 & 0 & 1 & 0\\ 
0 & 0 & 0 & 1\\ 
0 & 0 & 0 & 0
\end{pmatrix} 
\begin{pmatrix}
1 \\ 
D \\ 
D^2 \\ 
D^3 \end{pmatrix} g(x)$$
$$= \begin{pmatrix}
0 & g_0 & 0 & 0 \\ 
0 & g_1 & g_0 & 0\\ 
0 & g_2 & 2g_1 & g_0 \\ 
0 & g_3 & 3g_2 & 3g_1
\end{pmatrix} \begin{pmatrix}
1 \\ 
D \\ 
D^2 \\ 
D^3 \end{pmatrix} g(x)$$
Repeating operation $n$  times by $g(x+t)D_x$ and ultimately setting $t=0$, or, equivalently, extracting the first term in the resulting column vector of operators acting on $g(x)$, gives us
$$[g(x)D_x]^n .$$
The full repeated (infinite) matrix operator that acts on $g(x)$ is a rep in the basis $\frac{t^k}{k!}$ of 
$$[g(x+t)D_x]^n e^{tD_x}.$$
Letting
$\omega = f(x)$, $x = f^{-1}(\omega)$, and $g(x) = 1/f'(x)$, for $h(x)$ an arbitrary analytic function,
$$  \exp[ug(x+t)D_x]e^{tD_x}h(x)= \exp[u\frac{d}{d(\omega(x+t))}]h(x+t)$$
$$= \exp[u\frac{d}{d(\omega(x+t))}]h[f^{-1}[\omega(x+t)]]= h[f^{-1}[\omega(x+t)+u]] $$
$$ = h[f^{-1}[f(x+t)+u]]  = e^{tD_x} h[f^{-1}[f(x)+u]] $$
$$ = e^{tD_x} \exp[ug(x)D_x] h(x),$$
so (as we might surmise initially), equating coefficients of $u^n$,
$$ [g(x+t)D_x]^n e^{tD_x} = e^{tD_x} [g(x)D_x]^n,$$
implying the $m-$th component of the operator matrix rep, the operator coefficient of $\frac{t^m}{m!}$ is 
$$ D_x^m [g(x)D_x]^n,$$
in agreement with Grinberg's general result.
More simply, the proof follows through a double translation--a similarity transformation;
$$  e^{-tD_x}  [g(x+t)D_x]^n e^{tD_x} h(x) =[g(x)D_x]^n h(x)$$
since $D_x = \frac{d}{dx}= \frac{d}{d(x+t)} = D_{x+t}$.
