why my Linear Least-Squares does not fit right the data-points

Question

I'm getting not very fitted 2D-plane with such a code:

import numpy as np
from scipy import linalg as la
from scipy.linalg import solve

# data
f1 = np.array([1., 1.5, 3.5, 4.])
f2 =  np.array([3., 4., 7., 7.25])
# z = np.array([6., 6.5, 8., 9.])
A= X= np.array([f1, f2]).T

b= y= np.array([0.5, 1., 1.5, 2.]).T

##################### la.lstsq

res= la.lstsq(A,b)[0]
print(res)    
##################### custom lu

#custom OLS 
def ord_ls(X, y):
    A = X.T @ X
    b = X.T @ y
    beta = solve(A, b, overwrite_a=True, overwrite_b=True,
                 check_finite=True)
    return beta

res = ord_ls(X, y)
print(res)

##################### plot

# use the optimized parameters to plot the fitted curve in 3D space.
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Create 3D plot of the data points and the fitted curve
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(f1, f2, y, color='blue')
x_range = np.linspace(0, 7, 100)
y_range = np.linspace(0, 7,100)

X, Y = np.meshgrid(x_range, y_range)
Z = res[0]*X + res[1]

ax.plot_surface(X, Y, Z, color='red', alpha=0.5)
ax.set_xlabel('feat.1')
ax.set_ylabel('feat.2')
ax.set_zlabel('target')
plt.show()

# [0.2961165  0.09475728]
# [0.2961165  0.09475728]

though coeffs seems to be the same , still plot is distorting. Is there any explanation or correction? or some regularization, like in least-squares is needed ? or two features are collinear & that is the reason ? (I'm not very familiar with Linalg yet)

p.s. scipy-0.18.0 docs, - thought of LU-factorization at p.184 p, l, u = la.lu(A)

Answer 1

renamed topic, LS-fit (Mult.Lin.Regr.) for x, y,

working with intercept_column (reposted from my comment) - ok, I think that is because of 0 sum of squared errors (or residuals) for Lin Regr when we take into consideration Intersept - thus, minimization of SSE (of equation including intercept) becomes obviously fitted to OLS -- as from Normal Equation intercept, *coef = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y) we use equation (X.T.dot(X)).dot(β) = (X.T).dot(y) that minimizes L2-norm(xβ-y) (or in linalg terms (ax-b)):

import numpy as np
from scipy import linalg as la
from scipy.linalg import solve

# data
f1 = np.array([1., 1.5, 3.5, 4.])
f2 =  np.array([3., 4., 7., 7.25])
# z = np.array([6., 6.5, 8., 9.])
X= np.array([f1, f2]).T
y= np.array([0.5, 1., 1.5, 2.]).T
# adds x0 = 1 to each instance
X = np.column_stack((np.ones(len(y)),X))    # CORRECTION

##################### la.lstsq

x= la.lstsq(X,y)[0]
print("lstsq:", x)

##################### LU-decomposition
from scipy.linalg import lu
p, l, u = lu(X)
# reconstruct
print(p @ l @ u)
print(np.allclose(X - p @ l @ u, np.zeros((4, 3))))

#####################
from scipy.linalg import lu_factor, lu_solve

lu, piv = lu_factor(X.T@X)  # if not cov.matrix (as is square) taken then ValueError: expected square matrix
x = lu_solve((lu, piv), X.T @ y)
print("lu_solve: ", x)

##################### custom ls

# LR: y= ax-b, MLR: beta= (X.T @ X)^(-1)(X.T @ y)
A = X.T @ X # take covariance matrix to make it square for solve !
b = X.T @ y
beta = solve(A, b, overwrite_a=True, overwrite_b=True,
                 check_finite=True)
print("solve: ", beta)

############
# calculate normal equation (NOT STABLE)
theta_best = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
# best values for theta
intercept, *coef = theta_best
print(f"Intercept: {intercept}\n\
Coefficients: {coef}")

##################### plot

# use the optimized parameters to plot the fitted curve in 3D space.
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Create 3D plot of the data points and the fitted curve
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(f1, f2, y, color='blue')
x_range = np.linspace(0, 7, 100)
y_range = np.linspace(0, 7,100)

X, Y = np.meshgrid(x_range, y_range)
Z = beta[0] + beta[1]*X + beta[2]*Y

ax.plot_surface(X, Y, Z, color='red', alpha=0.5)
ax.set_xlabel('feat.1')
ax.set_ylabel('feat.2')
ax.set_zlabel('target')
plt.show()

N.B.

Don't calculate matrix inverse to solve linear systems The professional algorithms don't solve for the matrix inverse. It's slow and introduces unnecessary error. It's not a disaster for small systems, but why do something suboptimal? Basically anytime you see the math written as: x = A^-1 * b you instead want: x = np.linalg.solve(A, b)