<span class="comment">% j</span>

  S(i,i) = S(i,i) + 1e-06;    <span class="comment">% add small diagonal jitter for numerical reasons</span>
<span class="keyword">end</span> <span class="comment">% i</span>

dSdX(:,:,D+1,:) = -dSdX(:,:,D+2,:);
dSdX(:,:,D+1,:) = -dSdX(:,:,D+2,:);

<span class="comment">% 4) centralize moments</span>
S = S - M*M';
<span class="comment">%S(diag(S)&lt;0,diag(S)&lt;0) = 1e-6;</span>

<span class="comment">% 5) Vectorize derivatives</span>
dMds=reshape(dMds,[E D*D]);
dSdm=resha