code: def hessian_local_log_likelihood(self, x):
          
          # Observation likelihoods
          lmbda = np.exp(np.dot(x, self.C.T) + np.dot(self.inputs, self.D.T))
          return np.einsum('tn, ni, nj ->tij', -lmbda, self.C, self.C)
docs: d/dx  (y - lmbda)^T C = d/dx -exp(Cx + d)^T C
              = -C^T exp(Cx + d)^T C