Thanks! It’s interesting that shape=(1,1) works - and that at least gives a valid workaround for the problem.
I just ran across another workaround - replacing sd*np.eye(D) with (sd*np.eye(D))[:].
I’m pretty sure sd*np.eye(1) should be a 1x1 two-dimensional matrix (scalar*matrix=matrix), though I’m not sure how to debug whether it really is - certainly sd*np.eye(2) seems to be two-dimensional, and I don’t see why it should be any different.