.) function wrt Xj
      dmahai = dm1i-dm2i;
      dmahaj = dm1j-dm2j;
      
      if i==j
        LdXi = L.*(dmahai + bsxfun(@plus,kdX(:,i,d),kdX(:,j,d)'));
        dSdX(i,i,d,i) = beta(:,i)'*LdXi*beta(:,j);
      else
        LdXi = L.*(dmahai + bsxfun(@plus,kdX(:,i,d),zeros(n,1)'));
        LdXj = L.*(dmahaj + bsxfun(@plus,zeros(n,1),kdX(:,j,d)'));
        dSdX(i,j,d,i) = beta(:,i)'*LdXi*beta(:,j);
        dSdX(i,j,d,j) = beta(:,i)'*LdXj*beta(:,j);
      end
      
    end % d
    
    if i==j
      dSd