Higher level floating-point operations can only be expected to be accurate to numerical precision. Typical values are 1e-7ish for fp32 and 1e-15ish for fp64 per element, which seems broadly in range with your results.
Similar things would likely happen if you used numpy with different blas backends etc.
There isn’t a whole lot that PyTorch can do here.
Best regards
Thomas