@bwengals and @junpenglao thank you! This actually works well. @bwengals the theta is same for each x_i so what @junpenglao wrote works perfect! In general \theta \in \mathbb{R}^d, but I have already figured out how to modify the code so it works for a vector as well.
The mean function of my GP f is also a little bit peculiar. It looks like this:
m(x, \theta) = c(x) + \theta^Tv(x),
where c(x) is a constant and v(x) is a vector in \mathbb{R}^d.
I am implementing this by adding column to input_2, since I can calculate it explicitly for [x2, theta_star] and then having a linear mean function for gp2 with coeff = [1, 0, ...,0] and intercept = 0. I will also let the covariance function act only on the appropriate dimensions.
In the case of gp1, I am planning to include v(x)^T as rows in the design matrix x1 and then have a linear mean function with coeff = tt.concatenate([theta, [0,...,0]]) again with covariance matrix acting on particular dimensions.