Let $\mathcal{G}_n = \{ N(\mu,\Sigma) ; \mu \in \mathbb{R}^n, \Sigma > 0\}$ be the collection of Gaussian distributions on $\mathbb{R}^n$ with full support.

If $f : \mathbb{R}^n \to \mathbb{R}^k$ is measurable, $k\leq n$, and $f \# \gamma \in \mathcal{G}_k$ for all $\gamma \in \mathcal{G}_n$, is it true that $f$ is affine a.e.? Here $f\#\gamma$ denotes the pushforward of $\gamma$ under $f$.

In other words, if a function sends all Gaussian inputs to a Gaussian outputs, must that function be (essentially) affine?

If it is true, I guess it should be a classically known fact. But, a quick search didn't turn up any definitive answers. Related questions can be found here and here, but these only address the case for which $f$ sends a single Gaussian input to a Gaussian output (in which case there are obvious non-affine choices of $f$ that work).

It seems one should be able to expand $f$ in the Hermite basis (since $f\in L^2(N(0,I))$ follows from the hypothesis), and show coefficients for Hermite polynomials of degree two or higher must vanish by considering, e.g., the system of equations that arise by looking at moments of $f\#\gamma$, using the assumption that $f\#\gamma$ is Gaussian (so its higher moments are expressible in terms of first two moments). However, the computations get complicated...

A reference or counterexample would be appreciated.

Let $\mathcal{G}_n = \{ N(\mu,\Sigma) ; \mu \in \mathbb{R}^n, \Sigma > 0\}$ be the collection of Gaussian distributions on $\mathbb{R}^n$ with full support.

If $f : \mathbb{R}^n \to \mathbb{R}^k$ is measurable, $k\leq n$, and $f \# \gamma \in \mathcal{G}_k$ for all $\gamma \in \mathcal{G}_n$, is it true that $f$ is affine a.e.? Here $f\#\gamma$ denotes the pushforward of $\gamma$ under $f$.

In other words, if a function sends all Gaussian inputs to a Gaussian outputs, must that function be (essentially) affine?

If it is true, I guess it should be a classically known fact. But, a quick search didn't turn up any definitive answers. Related questions can be found on CrossValidated at Gaussian-to-gaussian transformations and on MathOverflow at Non-affine smooth transformation of Gaussian is Gaussian, but these only address the case for which $f$ sends a single Gaussian input to a Gaussian output (in which case there are obvious non-affine choices of $f$ that work).

It seems one should be able to expand $f$ in the Hermite basis (since $f\in L^2(N(0,I))$ follows from the hypothesis), and show coefficients for Hermite polynomials of degree two or higher must vanish by considering, e.g., the system of equations that arise by looking at moments of $f\#\gamma$, using the assumption that $f\#\gamma$ is Gaussian (so its higher moments are expressible in terms of first two moments). However, the computations get complicated….

A reference or counterexample would be appreciated.