      = Cdm*Mds(i,:) + Cds*Sds(ii,:);

  S(i,k) = S(i,i)*C; S(k,i) = S(i,k)';                      <span class="comment">% off-diagonal</span>
  SS = kron(eye(length(k)),S(i,i)); CC = kron(C',eye(length(i)));
  Sdm(ik,:) = SS*dCdm + CC*Sdm(ii,:); Sdm(ki,:) = Sdm(ik,:);
  Sds(ik,:) = SS*dCds + CC*Sds(ii,:); Sds(ki,:) = Sds(ik,:);
<span class="keyword">end</span>

<span class="comment">% 4. Calculate loss!</span>
L = 0; dLdm = zeros(1,D0); dLds = zeros(1,D0*D0); S2 = 0;
<span class="keyword">for</span> i = 1: