def gradients_XX_diag(self, dL_dK_diag, X):
        """
        Given the derivative of the objective dL_dK, compute the second derivative of K wrt X:

        ..math:
          \frac{\partial^2 K}{\partial X\partial X}

        ..returns:
            dL2_dXdX: [NxQxQ]
        """
        dL_dK_diag = dL_dK_diag.copy().reshape(-1, 1, 1)
        assert (dL_dK_diag.size == X.shape[0]) or (dL_dK_diag.size == 1), "dL_dK_diag has to be given as row [N] or column vector [Nx1]"

        l4 =  np.ones(X.shape[1])*self.lengthscale**2
        return dL_dK_diag * (np.eye(X.shape[1]) * -self.dK2_drdr_diag()/(l4))[None, :,:]