How well can the function $f(x_1,x_2)$ be approximated by $f_1(x_1)+f_2(x_2)$? Is it possible to find the upper bound of the following quantity?
$\min_{f_1,f_2}\int_0^1\int_0^1|f^*(x_1,x_2)-f_1(x_1)-f_2(x_2)|^2dx_1dx_2$.
where $f_1,f_2$ can be any continous functions.
For example,
$\min_{f_1,f_2}\int_0^1\int_0^1|x_1x_2-f_1(x_1)-f_2(x_2)|^2dx_1dx_2$.
Are there any related works in the literature?
Any comments are welcome.
 A: You can get a general solution by first neglecting the continuity constraints.
Fact: By properties of mean-square-error, we know that if $Z$ is a random variable with finite mean and variance, then the constant $c$ that minimizes
$$ E[(Z-c)^2]$$
is $c^*= E[Z]$.

Your problem (without continuity constraints) reduces to this:
Problem: Given an integrable function $f:[0,1]^2\rightarrow\mathbb{R}$, find a constant $c \in \mathbb{R}$ and functions $h_1:[0,1]\rightarrow \mathbb{R}, h_2:[0,1]\rightarrow\mathbb{R}$ to minimize:
$$ \int_0^1\int_0^1 (f(x,y) - c - h_1(x) - h_2(y))^2dxdy$$
subject to the constraints:
$$ \int_0^1 h_1(x)dx= \int_0^1 h_2(y)dy = 0$$
Solution:
The answer is
\begin{align}
c^* &= \int_0^1\int_0^1 f(x,y)dxdy\\
h_1^*(x) &=-c^*+\int_0^1 f(x,y)dy \quad \forall x \in [0,1]\\
h_2^*(y) &= -c^*+\int_0^1 f(x,y)dx \quad \forall y \in [0,1]
\end{align}
Proof: Define random variables $X, Y$ that are independent and uniformly distributed over $[0,1]$.
Suppose we are given $f,h_1, h_2$, where $h_1$ and $h_2$ meet the constraints.  Then
$$E[h_1(X)]=E[h_2(Y)]=0$$
We want to choose the best constant $c$  to minimize
$$ E[(f(X,Y)-h_1(X) - h_2(Y)-c)^2]$$
By the above fact we have
$$ c^* = E[f(X,Y)- h_1(X)-h_2(Y)] = E[f(X,Y)] = \int_0^1\int_0^1 f(x,y)dxdy$$
Now suppose we are given $c^*, f, h_2$, where $h_2$ satisfies the constraint $E[h_2(Y)]=0$. We want to choose $h_1$ to minimize:
$$ E[(f(X,Y)-h_2(Y)-c^* - h_1(X))^2] = \int_0^1E[(f(x,Y)-h_2(Y)-c^* - h_1(x))^2]dx $$
For each $x \in [0,1]$ we simply choose $h_1(x)$ to minimize the expectation inside the integral.  By the above fact we have for all $x \in [0,1]$:
\begin{align} 
h_1^*(x) &= E[f(x,Y)-h_2(Y)-c^*] \\
&= -c^*+\int_0^1f(x,y)dy 
\end{align}
and this also satisfies the desired constraint $\int_0^1 h_1^*(x)dx=0$.
Similarly, if we are given $c^*$, $f$, $h_1$ (where $h_1$ satisfies the constraint $E[h_1(X)]=0$) then for each $y \in [0,1]$ the best $h_2(y)$ is:
$$ h_2^*(y) = -c^*+\int_0^1 f(x,y)dx  $$
and this satisfies the desired constraint $\int_0^1 h_2^*(y)dy=0$.
Adding continuity:
If $h_1^*$ and $h_2^*$ are already continuous then we are done.  If they are piecewise continuous with a finite number of segments, we can approximate them arbitrarily closely by continuous functions, so the resulting optimal mean-square-error over continuous functions can get arbitrarily close to the case without the continuity constraint.
Notes:

*

*Note that this solution is consistent with the previous answer in the case when Fourier analysis is used.


*If we want, we can of course define $f_1^*(x) = h_1^*(x) +c^*$ and $f_2^*(y) = h_2^*(y)$.  Or we can add $c^*$ to the $h_2^*$ function, or add $c^*/2$ to both, etc.
A: Let's assume we restrict our attention to "nice" functions, where "nice" means "has a Fourier expansion on the interval $[0,1]$." In particular, write out the following expansions:
$$ f^{\ast}(x_1,x_2)=\sum_{k,\ell}c_{k,\ell}e^{2\pi i (kx_1+\ell x_2)} $$
$$ f_1(x_1)=\sum_k a_k e^{2\pi i kx_1}, \quad f_2(x_2)=\sum_{\ell} b_{\ell}e^{2\pi i\ell x_2}  $$
Then the squared $L^2$ distance between $f^{\ast}(x_1,x_2)$ and $f_1(x_1)+f_2(x_2)$ is
$$ \int_0^1\int_0^1 |f^{\ast}(x_1,x_2)-f_1(x_1)-f_2(x_2)|^2\,\mathrm{d}x_1\mathrm{d}x_2 $$
$$ =\sum_{k,\ell\ne0} |c_{k,\ell}|^2+\sum_{k\ne0} |c_{k,0}-a_k|^2+\sum_{\ell\ne0}|c_{0,\ell}-b_{\ell}|^2+|c_{0,0}-a_0-b_0|^2. $$
because the exponentials are an orthonormal basis for $L^2[0,1]$.
Thus, $a_k=c_{k,0}$ and $b_{\ell}=c_{0,\ell}$ for $k,\ell\ne0$ and $a_0+b_0=c_{0,0}$. In other words, $f_1(x_1)+f_2(x_2)$ ought to be all the terms of $f^{\ast}(x_1,x_2)$'s Fourier expansion that depend on only one of $x_1$ or $x_2$.
In particular, $\int_0^1 x\,\mathrm{d}x=\frac{1}{2}$ so $x$ has a Fourier expansion $x=\frac{1}{2}+(\cdots)$ so $xy$ is
$$ xy=(\frac{1}{2}+\cdots)(\frac{1}{2}+\cdots)=\frac{1}{2}\cdot\frac{1}{2}+\frac{1}{2}(x-\frac{1}{2})+(y-\frac{1}{2})\frac{1}{2}+(\cdots) $$
and therefore $f(x)+g(y)=-\frac{1}{4}+\frac{1}{2}x+\frac{1}{2}y$ minimizes the $L^2$ distance to $xy$, which is $\frac{1}{12}$.
Similarly, the terms of $f^{\ast}(x_1,x_2)$'s Fourier expansion that depend on only one of $x_1,x_2$ become:
$$ f_1(x_1)+f_2(x_2)=\int_0^1 f(x_1,t_2)\,\mathrm{d}t_2+\int_0^1 f(t_1,x_2)\,\mathrm{d}t_1-\int_0^1\int_0^1 f(t_1,t_2)\,\mathrm{d}t_1\mathrm{d}t_2. $$
