Let's say I need to learn 2 functions $f(x)$ and $g(z)$ (the input variables in vectors $x$ and $z$ may overlap, but in my specific case, they don't). However, I can only observe their sum, so the data is of the form $(x_k, z_k) \mapsto y_k = f(x_k) + g(z_k)$. I may use regression like $\tilde{f}(x; \theta_f)$ and $\tilde{g}(z; \theta_g)$, where $\theta_f$ and $\theta_g$ are parameters of models $\tilde{f}$ and $\tilde{g}$ of $f$ and $g$, for example they can be polynomials, NNs, etc. Given the data $\mathcal{D} = ((x_k,z_k), y_k)$ I can learn $\theta_f$ and $\theta_g$ in the usual way. With a new input $(x^\star, z^\star)$ I can predict the terms of the sum separately: $f(x^\star) \approx \tilde{f}(x^\star; \theta_f)$, similarly for $g$. First question: do you see any flaw in this method?

Now I want to use Gaussian Processes (GPs). I model $f$ as a GP $\mathcal{G}_f(0, k_f(\cdot, \cdot))$ with zero mean and kernel $k_f$, $g$ as a GP $\mathcal{G}_g(0, k_g(\cdot, \cdot))$ with zero mean and kernel $k_g$. As I understand, their sum is also a GP with zero mean and the sum of the kernels: $\mathcal{G}(0, k_f + k_g)$. I then learn the hyperparameters of $k_f$ and $k_g$ from the data as usual. With a new input $(x^\star, z^\star)$ I can predict the sum with the sum GP $\mathcal{G}$. The questions: Can I predict the individual terms $f$ and $g$? If yes, how? If I want to predict the mean of $f(x^\star)$ I would need the training observations of $f(x_k)$, something like: $\bar{f}(x^\star) = k_f(x^\star, X) K_f^{-1}(X,X) Y_f$ where everything can be calculated except $Y_f$ (it's the values of the $f$ term in $y_k$).

Thanks a lot!

Edit (to address a comment about the feasibility of this problem): while in general $f$ and $g$ are not identifiable from their sum, under certain conditions, this is possible. A simple example is if $f(x) = a x$ and $g(z) = b z$, then given a data set of $((x_k, z_k), f(x_k) + g(z_k))$, identifying $a$ and $b$, hence $f$ and $g$, is possible.

In my application, I consider discrete-time dynamic systems: $x(k+1) = f(x(k)) + g(z(k))$, where samples of $((x(k), z(k)), x(k+1))$ are given. Note that $x$ and $z$ are "independent" so $f(x)$ and $g(z)$ are "independent." If $f$ and $g$ are linear, it's well-known in system identification / control theory that, under certain conditions, $f$ and $g$ can be identified from data. For example, methods like DMDc can do this. For non-linear $f$ and $g$, methods like SINDYc can be used to learn these functions (under certain conditions).

My core question is that whether this is possible with GPs. The reason I want to separate the summands $f$ and $g$ is due to domain insights and that $f$ and $g$ will then have physical meanings.