Why does least squares linear regression perform so bad when switching from 2D to 3D line?
The matrix $X^TX$ in your example is close to a singular:

from numpy.linalg import inv, det

print(det(X.T @ X))

1.479949989718308e-06


which makes its inverse very large

inv(X.T @ X)

[[ 4.83118746e+11  2.66548273e+11 -1.83251938e+11]
 [ 2.66548273e+11  1.47061116e+11 -1.01104517e+11]
 [-1.83251938e+11 -1.01104517e+11  6.95093558e+10]]


This is the reason of the error in your question. The following expression

inv(X.T @ X) @ X.T @ X


should produce an identity matrix. But in your example it is

[[ 1.2265625   1.10289171  4.8378709 ]
 [ 0.02282715  0.61713325  0.95782902]
 [-0.01519775  0.25675456  0.33339437]] 


So, when you multiply this matrix by beta the result is far from beta.

Another way to look at this singularity is that your 3D problem is not actually 3D. It is still 2D. That's why you have a line, not a 2D surface in your plot. It can be converted to a 2D representation by a transformation of coordinates. If you make a matrix for a proper 3D problem, say

np.random.seed(3)
X = np.hstack((np.ones((100,1)),np.random.rand(100,2)*100))


then inv(X.T @ X) @ X.T @ X gives the correct result:

[[ 1.00000000e+00 -1.24233956e-13 -6.21724894e-14]
 [-1.90819582e-17  1.00000000e+00  6.38378239e-16]
 [-8.67361738e-18  9.28077060e-16  1.00000000e+00]]


and the formula $\beta=(X^TX)^{-1}X^Ty$ works well.