code: def gradient_log_joint(self, x):
          
          T, D = self.T, self.D_latent
          assert x.shape == (T, D)
  
          # Collect the Gaussian LDS prior terms
          _, h_init, _ = self.info_init_params
          _, _, _, h1, h2, _ = self.info_dynamics_params
          H_diag, H_upper_diag = self.sparse_J_prior
  
          # Compute the gradient from the prior
          g = -1 * symm_block_tridiag_matmul(H_diag, H_upper_diag, x)
          g[0] += h_init
          g[:-1] += h1
          g[1:] += h2
  
          # Compute gradient from the likelihood terms
          g += self.grad_local_log_likelihood(x)
  
          return g
docs: The gradient of the log joint probability.
  
          For the Gaussian terms, this is
  
              d/dx [-1/2 x^T J x + h^T x] = -Jx + h.
  
          For the likelihood terms, we have for each time t
  
              d/dx log p(yt | xt)