Why does numpy SVD return a non-normal vector when given a set of points exactly in a plane?

Question

I am attempting to fit a plane using numpy's SVD method np.linalg.svd(). As a test, I will use two sets of points in R³. In both cases, all points have value 100 in the 3rd dimension. Since the set of points I created is perfectly within the Z=100 plane, I expect that:

1. The third singular value will be 0 (within machine-precision).
2. The third column of vh will be [0, 0, 1] (within machine-precision).

Set 1

For some points in this set, the first and second values have magnitudes much larger than the third value.

pts = [[2345,-124, 100], [981, -123, 100], [4987,12345, 100], [-1324, 0, 100]]
svd = np.linalg.svd(pts)

The result here is roughly as-expected:

svd[1] produces array([13349.56221861, 2705.21722461, 158.26983058]). I would expect the third singular value to be closer to 0, since my points fit perfectly into a plane, but it's at least clear enough to indicate that the third column of svd[2] will be my plane normal vector.

svd[2] produces the following:

array([[-0.38833201, -0.92148669, -0.00778029],
       [-0.92117922,  0.38840419, -0.02389612],
       [ 0.02504185, -0.00211259, -0.99968417]])

Again, it's close. I would expect the first two dimensions of the 3rd column to be closer to zero (more like within machine-precision) but this is workable for my fitting application.

Set 2

For all points in this set, the first and second values have magnitudes smaller than the third value.

pts = [[57, 37, 100], [34, 37, 100], [11, -37, 100], [-11, 38, 100]]
svd = np.linalg.svd(pts)

This is where things started to look pretty weird.

svd[1] produces array([209.35774076, 64.78329726, 46.58820429]). This is surprising. The third singular value should be closer to 0.

svd[2] produces the following:

array([[-0.23396901, -0.2018738 , -0.95105493],
       [ 0.31100035,  0.91126958, -0.26993803],
       [-0.92116084,  0.35893555,  0.15042602]])

This is extremely unexpected. The third column of vh is quite far from [0, 0, 1]. Certainly well outside machine precision. It's actually closer to [1, 0, 0].

What is going on here? Is there something about how the SVD is implemented in numpy that does not give higher precision results? Am I just not using it right or misinterpreting the results?

Answer 1

EDIT: Looking at your edited post.

I think what you want from your SVD is the result after having shifted your points to their mean. You should get the exact results you want after doing so.

pts = np.array([[57, 37, 100], [34, 37, 100], [11, -37, 100], [-11, 38, 100]])
pts = pts - pts.mean(axis=0)
svd = np.linalg.svd(pts)

(array([[-0.45221016,  0.51373144,  0.15971725,  0.71139045],
        [-0.31606121,  0.06525962,  0.8397871 , -0.43658232],
        [ 0.83162358,  0.2416381 ,  0.42832508,  0.25797456],
        [-0.06335221, -0.82062916,  0.29289188,  0.48658876]]),
 array([67.13174076, 47.06197384,  0.        ]),
 array([[-0.39738768, -0.91765082, -0.        ],
        [ 0.91765082, -0.39738768, -0.        ],
        [ 0.        ,  0.        ,  1.        ]]))

The reason why this is different is that SVD finds the best linear (not affine) transformation to fit your points, so it cannot handle a shifted plane. In your case when you set the distance in z further away, you're incentivizing the plane which has to go through the origin to lay itself upright and cut through the points, rather than lie flat.

Thus when the values of your z coordinates are much higher than your x and y coordinates you'll get planes that look unexpected.