I see, that does make things harder! I’d say there’s no perfect answer here, but some suggestions below.
The first thing I would probably check is whether it might be reasonable to make a Gaussian approximation to your Poisson likelihood. If your counts are large, this might be reasonable (see e.g. Poisson normal approximation error). Then you could marginalise out the GP, at which point MAP estimation for the hyperparameters usually works fine. By the way, I just noticed your toy data is apparently continuous – it should be discrete for a Poisson! That might also be causing some numerical issues.
Another option could be to see if you could estimate the hyperparameters with NUTS on some smaller subset of data, and then turn them into a point estimate, e.g. by just using their mean, when fitting / predicting on your full set of time series. It might also be worth checking out the JAX-based samplers (see here) and trying NUTS with those on a GPU – there can definitely be big speedups.
You could also try a sparse approximation to the GPs: Sparse Approximations — PyMC3 3.10.0 documentation . The idea there would roughly be that you approximate the GPs with 300+ points with ones with say 100, which can speed things up.
So no great answer from me here, sorry. Maybe someone else has a clear answer, but scaling GPs is often a bit tricky. You could do things with variational inference too, there are special techniques developed for GPs (e.g. here) that can work well, but you would likely have to do this outside pymc3.