Learning Poisson State Space Model parameters for large number of samples

I guess at this point it would be good if you could go into more detail about your exact model equations? You started from:

\begin{align}X_{t+1} &\sim N(\mu_t, Q) \\ \mu_t &= A x_t \\\ y_t &\sim Poisson(\lambda_t) \\ \lambda_t &= C \exp(X_t) \end{align}

I made some small changes from your original equation, because I don’t make a distinction between X and U. They would just be concatenated together into what I call X, while matrices you called A and B would be block concatenated what I called A. I made some simple choices for A, C, Q, but they’re not necessarily appropriate at all. Actually, you can accomplish what I did much easier just using X_t \sim \text{GaussianRandomWalk}(\cdot).

So what kind of dynamics are you modeling? ARIMA? Seasonal cycles? Random walk with drift? Exponential smoothing?

For multivariate, the original model is multivariate. You shouldn’t be using a 3d tensor as the observed data, you just want to add more columns to X, y, or both.