 = dm1j-dm2j;

      <span class="keyword">if</span> i==j
        LdXi = L.*(dmahai + bsxfun(@plus,kdX(:,i,d),kdX(:,j,d)'));
        dSdX(i,i,d,i) = beta(:,i)'*LdXi*beta(:,j);
      <span class="keyword">else</span>
        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);
      <span class="keyword">end</span>

    <span class="key