To better understand SVD, I'm trying to recreate the values for U, S, and V using straight numpy, but I can't get the same results.

According to numpy's documentation for its implementation of SVD, it returns values for U, S, and V such that you can recreate your original dataset in the following way:

(U*S) @ V = original_data

My understanding is that the values for U are derived from the eigenvectors of $ XX^T $, V from the eigenvectors of $X^TX$, and they share the same eigenvalues, and their square root comprise the values for S.

So, to recreate these on my own I'm doing the following:

from sklearn.datasets import load_boston
boston = load_boston()
X = boston.data
u_vals, u_vecs = np.linalg.eig(X.dot(X.T))
v_vals, v_vecs = np.linalg.eig(X.T.dot(X))

My thinking is that u_vecs[:, :13] @ np.diag(v_vals**.5) @ v_vecs.T ought to do the trick.

However, this gives very different results.

Even more, I've checked further, and the only difference between the values I have from numpy's eigensolver and the ones returned from its svd method are the signs involved in the eigenvectors used.

However, my (perhaps naive) understanding is that this shouldn't make a difference. However, I'm clearly mistaken in something.

Is it simply a case of using a different eigensolver, or is there a deeper misunderstanding that I'm not seeing?