code: def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S, target):
          
          self._psi_computations(Z, mu, S)
          d_var = 2.*self._psi2 / self.variance
          # d_length = 2.*self._psi2[:, :, :, None] * (self._psi2_Zdist_sq * self._psi2_denom + self._psi2_mudist_sq + S[:, None, None, :] / self.lengthscale2) / (self.lengthscale * self._psi2_denom)
          d_length = -2.*self._psi2[:, :, :, None] * (self._psi2_Zdist_sq * self._psi2_denom + self._psi2_mudist_sq + S[:, None, None, :] * self.inv_lengthscale2) / (self.inv_lengthscale * self._psi2_denom)
          target[0] += np.sum(dL_dpsi2 * d_var)
          dpsi2_dlength = d_length * dL_dpsi2[:, :, :, None]
          if not self.ARD:
              target[1] += dpsi2_dlength.sum() # *(-self.lengthscale2)
          else:
              target[1:] += dpsi2_dlength.sum(0).sum(0).sum(0)
docs: Shape N,num_inducing,num_inducing,Ntheta