Estimating a function $f$ of a random vector $\mathbf{x}$ by a subset of the coordinates of $\mathbf{x}$ after a rotation of the input space Suppose I have
$$h=f(\mathbf{x})$$
with $f$ a deterministic function and $\mathbf{x}=(x_1,\ldots,x_n)$ a random vector of known distribution. I'm not using the capital letter notation for random variables, because I use it already to refer to deterministic matrices (e.g., $\mathbf{W}$), and I noticed that the two notations clash. If you have a better idea, please let me know :) 
Back to the question: let $\mathbf{W}$ be an non-random $n$ x $n$ orthogonal matrix, and let $\mathbf{W_1}$ indicate the first $m$ columns of $\mathbf{W}$ (thus it's an $n$ x $m$ matrix), while $\mathbf{W_2}$ indicate the remaining $n-m$ columns (thus it's $n$ x $(n-m)$). Then $ \mathbf{W}=[\mathbf{W_1} \quad \mathbf{W_2}]$. Since the columns of $\mathbf{W}$ are orthonormal vectors, it can be shown that
$$\mathbf{x}=\mathbf{W}\mathbf{W}^T\mathbf{x}=\mathbf{W_1}\mathbf{W_1}^T\mathbf{x}+\mathbf{W_2}\mathbf{W_2}^T\mathbf{x}=\mathbf{W_1}\mathbf{y}+\mathbf{W_2}\mathbf{z}$$
where we introduced $\mathbf{W_1}^T\mathbf{x}=\mathbf{y}$ and $\mathbf{W_2}^T\mathbf{x}=\mathbf{z}$. For some reasons, I would like to estimate $h$ with a function of $\mathbf{y}$ only. Ideally, the minimum squared loss estimation of $h$ by a function of $\mathbf{y}$ is given by the conditional expectation of $h$ with respect to $\mathbf{y}$:
$$g(\mathbf{y})=\mathbb{E}[f(\mathbf{x})|y]=\int f(\mathbf{W_1}\mathbf{y}+\mathbf{W_2}\mathbf{z})p(\mathbf{z}|\mathbf{y})d\mathbf{z}$$
To estimate this function from data, for various values of $\mathbf{y}$, a possibility would be to sample from $p(\mathbf{z}|\mathbf{y})$ and compute the integral by Monte Carlo. Sampling from $p(\mathbf{z}|\mathbf{y})$ can be quite complex. 
However, when $\mathbf{x}$ has a joint Gaussian distribution, clearly also $\mathbf{z}$ and  $\mathbf{y}$ have a joint Gaussian distribution. In a couple of practical applications, where $\mathbf{x}=(x_1,\ldots,x_n)$ were not only jointly Gaussian, but also independent, I found out that just a simple polynomial regression of $h$ on $\mathbf{y}$ gave very good results, without going through all the hassle of having to sample from $p(\mathbf{z}|\mathbf{y})$. Is this just a coincidence, or is there some good theoretical reason for this? This doesn't seem to me the classical regression setting. Consider for example linear regression: then I would have that the error term is
$$h-\boldsymbol{\beta}\cdot\mathbf{y}=f(\mathbf{W_1}\mathbf{y}+\mathbf{W_2}\mathbf{z})-\boldsymbol{\beta}\cdot\mathbf{y}= \epsilon $$
I don't think this error term is normally distributed, unless $f$ has a very simple structure.
EDIT: upon request from a user, I give more details on the method to which this question is related, active subspaces.
The idea is that we have a complex nonlinear numerical model $h=f(\mathbf{x})$, which is very expensive or complex to evaluate. For example, think of the mean Earth temperature in 50 years from now, as computed by a very complex climatology code. We would like to learn more about the model, such, for example, to understand which are the most influential inputs, which is the probability that $h$ exceeds a given threshold, etc. These kind of studies require running the black box code giving $h=f(\mathbf{x})$ multiple times for different values of $x$, and as such they're very expensive. If we can reduce the dimensionality of the problem, then these studies become much less expensive. For example, MCMC methods converge better in low dimensions, or creating a statistical model for $f(\mathbf{x})$, such as a Gaussian process, is easier in low dimensions than in high dimensions. 
In particular, we look for coordinate changes in the input space such that most of the variability of the output (not the variability of the predictors) is due to the first few coordinates in the new rotated space. In a sense, the method has some (weak) similarity to PLS regression (in the specific case where we have only one scalar response), while it's very different from PCR. Trying to keep the gory details to a minimum, the method considers the real symmetric positive semidefinite $n$ x $n$ matrix
$$\mathbf{C}=\mathbb{E}[\nabla f \nabla f^T]=\int \nabla f \nabla f^T p(\mathbf{x})d\mathbf{x}$$
Then, eigenvalues and eigenvectors of $\mathbf{C}$ are computed:
$$\mathbf{C}=\mathbf{W}\Lambda\mathbf{W}^T$$
In the above eigendecomposition, we order eigenvalues and eigenvectors so that
$$\Lambda=diag(\lambda_1,\ldots,\lambda_n)\quad with \quad\lambda_1\ge\ldots\ge\lambda_n\ge 0$$ 
We then partition $\mathbf{W}$ in $\mathbf{W}_1$ and $\mathbf{W}_2$  as above. 
$\mathbf{W}_1$ contains the $n$ leftmost columns of $\mathbf{W}$, which are related to the larger eigenvalues. It can be shown that $f(\mathbf{x})$ varies more along the $\mathbf{y}$ coordinates than along the $\mathbf{z}$ coordinates. For this reason, we try to approximate $f$ by a function of the new  $\mathbf{y}=\mathbf{W}_1^T\mathbf{x}$ coordinates only. Under some conditions, choosing $n$ according to some criteria allows to estimate $f$ accurately. However, the proofs approximate $f$ by its conditional expectation with respect to $\mathbf{y}$. In my applications I noticed that simple polynomial regression would work just fine, so I was wondering if that was just dumb luck or if there's some hidden truth here.  
 A: I will try to give you indications as of why this is not surprising theoritically. I would need more details about your results to be more precise.


*

*Polynomials can fit any function to any precision if you set a number of degree high enough. So if you are just looking at the fit and not any type of prediction metric, indeed polynomial with the correct number of degree should perform well. Considering that you are also using only top features in term of variance the fit should be decent. You are disregarding $z$ which should not contain much information and thus should not decrease much your model performances.

*Normality of errors is not required for a regression to perform well. Your model can be extremely good without having normal errors. However, when errors are normally distributed, we can derive a lot more nice results, such as optimality of the solution, Confidence interval etc... Here is a post on that topic: How does linear regression use the normal distribution? and here is the relevant wikipedia article https://en.wikipedia.org/wiki/Ordinary_least_squares
