SVD Inversion, Moore Penrose and and LSQ give different answers using Numpy-CodePudding

I am solving a matrix using different methods. According to my interpretation of the numpy descriptions, all three of my tested methods (SVD inversion, moore-penrose inversion, and Least Squares) should result in the same answer. However, the SVD inversion results in a very different answer. I cannot find a mathematical reason for this in Numerical Recipes. Is there a Numpy implementation nuance that is causing this?

I am using the following code on Python 3.8.10, Numpy 1.21.4, in a jupyter notebook

y = np.array([176, 166, 194])
x = np.array([324, 322, 376])
x = np.stack([x, np.ones_like(x)], axis=1)

# Solve the matrix using singular value decomposition
u, s, vh = np.linalg.svd(x, full_matrices=False)
s = np.where(s < np.finfo(s.dtype).eps, 0, s)
manual_scale, manual_offset = vh @ np.linalg.inv(np.diag(s)) @ u.T @ y
display(manual_scale, manual_offset, manual_scale * x   manual_offset)

# Solve the matrix using Moore-Penrose Inversion
# Manually
manual_scale, manual_offset = np.linalg.inv(x.T @ x) @ x.T @ y
display(manual_scale, manual_offset, manual_scale * x   manual_offset)
# Using supplied numpy methods
manual_scale, manual_offset = np.linalg.pinv(x) @ y
display(manual_scale, manual_offset, manual_scale * x   manual_offset)

# Solve using lstsq
((manual_scale, manual_offset), residuals, rank, s) = np.linalg.lstsq(x, y)

display(manual_scale, manual_offset, manual_scale * x   manual_offset)

The output (edited for clarity) is then

'SVD'
0.6091639943577222
29.167637174498772
array([[226.53677135,  29.77680117],
       [225.31844336,  29.77680117],
       [258.21329905,  29.77680117]])

'Manual Moore-Penrose'
0.4388335704125341
29.170697012800005
array([[171.35277383,  29.60953058],
       [170.47510669,  29.60953058],
       [194.17211949,  29.60953058]])

'Moore-Penrose'
0.43883357041251736
29.170697012802187
array([[171.35277383,  29.60953058],
       [170.47510669,  29.60953058],
       [194.17211949,  29.60953058]])

'LSTSQ'
/tmp/ipykernel_261995/387148285.py:24: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
  ((manual_scale, manual_offset), residuals, rank, s) = np.linalg.lstsq(x, y)
0.43883357041251814
29.17069701280214
array([[171.35277383,  29.60953058],
       [170.47510669,  29.60953058],
       [194.17211949,  29.60953058]])

As you can see three later methods get the same result, yet the manual svd calculation is different. What is going on?

CodePudding user response：

You are missing a transpose of vh. The SVD solution should be

manual_scale, manual_offset = vh.T @ np.linalg.inv(np.diag(s)) @ u.T @ y

By the way, you can simplify the inverse of the diagonal factor:

manual_scale, manual_offset = vh.T @ np.diag(1/s) @ u.T @ y

(That assumes there are no zeros in s.)

CodePudding user response：

For the next person who needs this, the fixed code is below. Thanks Warren!

y = np.array([176, 166, 194])
x = np.array([324, 322, 376])
x = np.stack([x, np.ones_like(x)], axis=1)

# Solve the matrix using singular value decomposition
u, s, vh = np.linalg.svd(x, full_matrices=False)
s = np.where(s < np.finfo(s.dtype).eps, 0, s)
manual_scale, manual_offset = vh.T @ np.diag(1/s) @ u.T @ y
display('SVD')
display(manual_scale, manual_offset, manual_scale * x   manual_offset)

# Solve the matrix using Moore-Penrose Inversion
# Manually
manual_scale, manual_offset = np.linalg.inv(x.T @ x) @ x.T @ y
display('Manual Moore-Penrose')
display(manual_scale, manual_offset, manual_scale * x   manual_offset)
# Using supplied numpy methods
manual_scale, manual_offset = np.linalg.pinv(x) @ y
display('Moore-Penrose')
display(manual_scale, manual_offset, manual_scale * x   manual_offset)

# Solve using lstsq
display('LSTSQ')
((manual_scale, manual_offset), residuals, rank, s) = np.linalg.lstsq(x, y)

display(manual_scale, manual_offset, manual_scale * x   manual_offset)