Gaussian Process Regression-- reparameterizing to us class

I am trying to speed up a Gaussian Process Regression model. I notice there is a method. However, I have not been using the class in my model definition. I assume the first step to speeding up my GP model is to re parameterize using the class. I am confused on how to do this in my case. Any help or suggestions would be appreciated

with pm.Model() as GP:
   """ This model uses a precalculated distance matrix to model dependencies between individual models that are geographically close, defined in the distance matrix. """     
# ====== covariance matrix ========
        etasq = pm.HalfCauchy('etasq', 1) # sets maximum covariance ij
        rhosq = pm.HalfCauchy('rhosq', 1e-5) # determines rate of covariance decline between farms
        # will have very small posterior because distances are so large in Dmatsq
        Kij = etasq*(np.exp(-rhosq*Dmatsq)+np.diag([.01]*Num_farms))
# ========== gaussian process prior =========
        w = pm.MvNormal('w', mu=np.zeros(Num_farms), cov=Kij, shape=Num_farms) 
        a = pm.Normal('a', 0, .5, shape=Num_farms) 
 # ========== Linear Model =============== 
        u = a[zero_farm_idx]  + w[zero_farm_idx]*df.zprsum3 
#============ Likelihood ==============
        sigma = pm.Uniform('sigma', 0,2)
        y = pm.Normal('NDVI', mu=u, sd=sigma, observed=df.zNDVIsum3)

You need to first express the above as a, and then replace the MvNormal with with the kernel as an input.

Hey Thanks for the reply!

I see. I think I am running into problems because I am trying to use the GP as a part of a larger GLM varying slopes/intercepts model. I asked a better question that gets at this.

Question reframed