wolffd@0: function dPsi = cpGradientDiag(X, S, batchsize)
wolffd@0: 
wolffd@0:     dPsi    = diag(X * S * X') / batchsize;
wolffd@0: 
wolffd@0: end