For anybody else with a similar issue, creating the distribution by explicitly specifying scale_tril is much more efficient. The GPU version is still a bit slower than the CPU version, but not by much. The documentation says that if a covariance matrix is passed to create the distribution, the corresponding triangular matrices are computed with a Cholesky decomposition. My guess is that the internal Cholesky decomposition on GPU tensors is much slower than on CPU tensors for some reason.