Solution (source code)

= Solution

Put $L=\sup_x\|\nabla f(x)\|_\infty$, the uniform bound on the <supremum norm> of the <gradient>. The <fundamental theorem of calculus> on a line segment, together with the <inner product> inequality $|\langle a,b\rangle|\leq\|a\|_\infty\|b\|_1$, gives
$$
|f(x)-f(y)|\leq L\|x-y\|_1,\qquad |f(x)|\leq|f(0)|+L\|x\|_1.
$$
Thus $f$ is <Lipschitz continuous> with respect to the <L1 norm>. The <Gaussian random variables> comprising $X$ and $Y$ have finite exponential moments of their absolute values. More explicitly, for $c\geq0$, <Hölder's inequality> gives $\mathbb E e^{c\|X\|_1}\leq\prod_{j=1}^n(\mathbb E e^{nc|X_j|})^{1/n}<\infty$. Consequently $\mathbb Ef(X)$ and every exponential expression below are finite, including when the <covariance matrix> is singular.

Use the <convex function> $\Psi(v)=e^{\lambda|v|}$. Since $Y$ is an <independent> copy of $X$, <Jensen inequality> for the <conditional expectation> gives
$$
\mathbb E\Psi(f(X)-\mathbb Ef(X))\leq\mathbb E\Psi(f(X)-f(Y)).
$$
For $0\leq\theta\leq\pi/2$, define the <Gaussian rotation of independent copies>
$$
U_\theta=X\sin\theta+Y\cos\theta,\qquad V_\theta=X\cos\theta-Y\sin\theta.
$$
Writing $\Sigma=\operatorname{Cov}(X)$, both rotated <covariance matrices> are $\Sigma$ and their cross-<covariance> is zero. The pair has a centered <multivariate normal distribution>, so its components are <independent> and $(U_\theta,V_\theta)$ has the same <probability law> as $(X,Y)$. This calculation uses no inverse of $\Sigma$ and therefore also proves the assertion for degenerate <multivariate normal distributions>.

The path $\theta\mapsto U_\theta$ runs from $Y$ to $X$. The <chain rule> and the <fundamental theorem of calculus> give the pathwise identity
$$
f(X)-f(Y)=\int_0^{\pi/2}\langle\nabla f(U_\theta),V_\theta\rangle\,d\theta.
$$
Apply <Jensen inequality> to the uniform <probability measure> on $[0,\pi/2]$. Then use <Tonelli theorem> and the <Gaussian rotation of independent copies> to obtain
$$
\begin{aligned}
\mathbb E\Psi(f(X)-f(Y))&\leq\frac2\pi\int_0^{\pi/2}\mathbb E\Psi\!\left(\frac\pi2\langle\nabla f(U_\theta),V_\theta\rangle\right)\,d\theta\\
&=\mathbb E\Psi\!\left(\frac\pi2\langle\nabla f(X),Y\rangle\right).
\end{aligned}
$$
This proves the <Gaussian rotation interpolation inequality> in the required form:
$$
\boxed{\mathbb E e^{\lambda|f(X)-\mathbb Ef(X)|}\leq\mathbb E e^{(\lambda\pi/2)|\langle\nabla f(X),Y\rangle|}.}
$$
The stated almost-sure domination now gives an <exponential moment of an absolute standard normal variable> bound. If $a=\lambda\pi/2$, completing the square in the <standard normal density> yields
$$
\mathbb E e^{a|Z|}=\frac2{\sqrt{2\pi}}\int_0^\infty e^{az-z^2/2}\,dz=2e^{a^2/2}\Phi(a)\leq2e^{a^2/2},
$$
where $\Phi$ is the <standard normal distribution function>. Hence <Markov inequality> gives, for every $\lambda>0$,
$$
\mathbb P(|f(X)-\mathbb Ef(X)|>u)\leq2\exp\!\left(-\lambda u+\frac{\pi^2\lambda^2}{8}\right).
$$
The exponent is minimized at $\lambda=4u/\pi^2>0$. Therefore the resulting <Gaussian tail bound> is
$$
\boxed{\mathbb P(|f(X)-\mathbb Ef(X)|>u)\leq2e^{-2u^2/\pi^2}\qquad(u>0).}
$$
No <independence> of $Z$ from $X$ or $Y$ is needed: only the given almost-sure domination is used.