I believe pytorch is correct. Basically the problem is this:
y=sigmoid(x @ w.T)
dy/dx=y*(1-y) @ w
where @ is matrix mul, * is element-wise mul
to calculate second order derivatives d2y/dx2, you can draw a computational graph starting from x and ending at dy/dx, then backpropagates to x. The final answer is this:
d2y/dx2 = z @ w.T * (1 - 2y) * y * (1-y) @ w
Of course, when I tried to derive the above equations, I got the wrong answer too… if anyone is good at math please let me know how it’s derived…