Polynomial least squares fitting, why is the iterative solver better?

Question

Take a look at the Wikipedia page for Polynomial_regression, the method where they use the Vandermonde matrix. I looked around the web this and method is mentioned in quite a few places but for some reason Numpy's polyfit works way better. You can go ahead see for yourself with the code below, can someone explain why and what algorithm Numpy is using?

import numpy as np
import matplotlib.pyplot as plt


x = np.linspace(0,1,10000)
y = np.exp(-80*(t-0.5)**2) #relatively sharp gaussian function, centred at 0.5

plt.figure(figsize=(9,5),dpi=150)
plt.plot(x,y,label="Gaussian")
z = np.polyfit(x, y, 40)
p = np.poly1d(z)
plt.plot(x,p(x),label="Numpy's PolyFit")

xMatrix = []
yMatrix = np.transpose(y)
n = 40

for i in range(len(x)):
    temp = []
    for j in range(n):
        temp.append(x[i]**j)
    xMatrix.append(temp)

xMatrixT = np.transpose(xMatrix)
dot1 = np.matmul(xMatrixT, xMatrix)
dot2 = np.matmul(xMatrixT, yMatrix)
dotInv = np.linalg.inv(dot1)
coefficients = np.matmul(dotInv, dot2)

p = np.poly1d(coefficients[::-1])#numpy's poly1d coefficients are descending so invert the order
plt.plot(x,p(x),label="Vandermonde Regression")
plt.legend()

Answer 1

The polyfit function in numpy ultimately ends up calling the LAPACK function GELSD(), which uses a numerically stable singular value decomposition approach to finding a minimum norm least-squares solution. Furthermore, the polyfit function scales the least-squares problem before solving it.

Your code uses the "normal equations" approach and solves

$X^{T}X \beta = X^{T}y$.

In forming the product $X^{T}X$ the condition number of the $X$ matrix is squared and the solution becomes much more sensitive to roundoff errors than solving the least squares problem using orthogonal transformation methods (QR or SVD).

You could try linalg.lstsq(X,y,-1) to solve the least-squares problem by SVD without forming the $X^{T}X$ matrix.