Approximating interactions in OLS Suppose I have a simple linear model with two variables and their interactions:
$y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_1 x_2 + \beta_4 x_1^2 + \beta_5 x_2^2 + \epsilon$
where the $\beta_i$ can, of course, be estimated by OLS. Assume that the standard Gauss-Markov assumptions of OLS hold, and that data points are sampled i.i.d. My question is the following: Is there a way to approximate this model such that the approximation is separable in $x_1, x_2$? 
Clearly, I could just omit the term $\beta_3 x_1 x_2$ and obtain a separable model at the cost of accuracy, but I'm wondering if there is an approach which would better approximate the original model? I'm only interested in predicting $y$, so the metric by which an approximation would be assessed would be the additional error in introduces in predicting $y$. In other words, I'm interested in an approximation which minimizes the difference in sum-squared-residuals relative to OLS.
The approximation need not be linear in its parameters (or even parametric). This seems related to the idea behind partitioned regression. I also know that orthogonal polynomials are commonly applied to function approximation, but again, am not sure if it would be useful here.
If anyone has any suggestions about where to look for some ideas/theory I'd be happy to hear them.
 A: If you really need separability for some reason, you can often get a reasonable approximation by transformation of both response and predictors; however the resulting relationships - while still smooth - will generally require additive modeling. 
You might, for example look into approaches like ACE (alternating conditional expectations -- see Breiman, L. and Friedman, J. H. (1985). "Estimating optimal transformations for multiple regression and correlation.",
JASA
80
, 580-78. This is implemented as ace in package acepack for R), for example. 
You can do simpler things; if you fit reasonably flexible models in the two predictors (such as via splines for example), you can use Box-Cox transformations (and even estimate the $\lambda$ parameter) to sometimes get reasonable transformations for the response. While less flexible this may be easier to interpret and explain.
There are a wide variety of other approaches to getting approximate transformations for additivity -- indeed it's essentially the same notion as Lafay anamorphosis in nomography which dates back well over a century.
