As luck would have it Iñigo Quilez wrote an article about this which popped up in my Facebook feed this morning.
Hardware texture intepolation is fast and convenient. It is bilinear
(plus mipmapping), and despite it can be somehow improved, it works
great for most cases. Most cases being texture mapping of surfaces
with color/albedo, normal and specular textures. Generally these types
of textures don't need lots of color fidelity (where "color" means
pixel value/content), and therefore the hardware's texture
interpolation units implement some simple interpolation based on 24.8
fixed point arithmetic. That's what all GPUs' hardware implments these
days.
A 24.8 precision texture interpolator means that there's a maximum of
256 intemediate values possible between two adjacent pixels of a
texture. 256 values are a lot for albedo textures for sure, but often
in computer graphics textures encode not only surface properties, but
they serve as LookUp Tables (LUT), heighfields (for terrain
rendering), or who knows what. In those cases, you can find yourself
easily lacking more resolution than 256 values between pixels.
Here's an example image from the article:
(Top) Hardware interpolation: visually smooth, but the barcode shows the derivative is spiky.
(Bottom) Manual interpolation: the constant colour in the derivative bar shows we're not getting step artifacts anymore.
So there's the answer. Hardware filtering is mainly intended to bridge sub-texel gaps for surface textures. So assuming our textures have enough resolution & mip maps for how close they are to the camera, we usually only have to interpolate a few steps - to deal with a texel center that's just a bit off from our rendered pixel center, or to fill a gap between two adjacent texels that are slightly more than 1 rendered pixel apart. For most cases, if we need anywhere near all 256 interpolation steps shown at once, it's because the player has crammed their camera right up against a wall, or we're trying to massively stretch a texture that's too low-res for the use it's being put to.
But if we're doing procedural graphics, it's often tempting to use a low-resolution texture stretched wide to give us low-frequency gradients, using the texture interpolator to give us free bilinear interpolation between samples. For colours this is often okay, since the human eye is pretty bad at spotting banding artifacts unless they're extreme. But if we care about the precise value at each sample - like for heightfields or sharp thresholds as you're doing here, then these staircase & plateau artifacts creep in.
So what can we do about it?
The most general-purpose solution (and what Iñigo Quilez recommends in his article) is exactly what you're doing: sample the texture four times and interpolate the results yourself in the shader, where you have full floating point precision.
If the case you've shown is representative of your real applications though, there may be more efficient methods:
Increase the resolution of your input texture, say to 16x16. Now the texture interpolation from a single fetch only needs to cover 1/15th as much ground with its 256 intermediate values, so you have an effective range of 3840 intermediates from one edge to the other. Much better! The memory taken up by this texture is still tiny, so you needn't feel bad about being wasteful.
If you only ever need to do bilinear filtering between the same 4 corners, you might as well do away with the texture and do it all in shader variables. This way you don't pay for even a single texture fetch. You can profile to double-check, but I'd wager this might be even be faster than using the "free" interpolation from the texture pipeline, since there's no contention or caching overhead.