Let $$ X=\begin{pmatrix} 0 & & & &\\ 1 & 0 & & &\\ & 1 & 0 & &\\ & & \ddots & \ddots &\\ & & & 1 & 0 \end{pmatrix}, D=\begin{pmatrix} 0 & 1 & & &\\ & 0 & 2 & &\\ & & \ddots & \ddots &\\ & & & 0 & n\\ & & & & 0 & \end{pmatrix}\in\mathbb{R}^{(n+1)\times(n+1)} $$ and $$ \mathbf{e}_0=(1,0,\cdots,0)^\mathrm{T},\mathbf{e}_n=(0,\cdots,0,1)^\mathrm{T}\in\mathbb{R}^{n+1} $$ Prove that $$ (X-D)^n\mathbf{e}_0=\exp\left(-\frac{D^2}{2}\right)\mathbf{e}_n $$
Background: it directly comes from two representations of (probabilists') Hermite polynomial: $H_n(x)=(x-D)^n\cdot 1$ and $H_n(x)=e^{-D^2/2}\cdot x^n$, where $D$ is the differential operator. The equivalence of the two representation is not difficult but it seems techniques like generating functions are necessary at least for me. I wonder whether the identity can be proved by linear algebra only without explicitly going into Hermite polynomials.