Python无内置方法实现最小二乘多项式拟合的问题求助

我正在完成一项练习：从含噪数据集中用最小二乘法找到指定阶数的最优拟合多项式。目前卡在推导到求解线性方程组的步骤，希望了解具体推导过程，或者获取能生成可传入我现有LU分解/高斯分解程序的矩阵的Python代码。注：我已经实现了三次样条、LU分解、高斯分解的Python程序。

我曾直接对数据集应用高斯消元/LU分解，但发现遗漏了关键步骤；同时我也不清楚三次样条在此问题中的作用。

编辑1：尝试高斯消元与LU分解代码 #

我先写了高斯消元的代码，原本以为输入数据集和阶数就能输出拟合多项式的系数，但教授的评测工具显示误差很大：

import numpy as np
import math

def swapRows(v,i,j):
    if len(v.shape) == 1:
        v[i],v[j] = v[j],v[i]
    else:
        v[[i,j],:] = v[[j,i],:]
def swapCols(v,i,j):
    v[:,[i,j]] = v[:,[j,i]]
def gaussPivot(a,b,tol=1.0e-12):
    n = len(b)
    # Set up scale factors
    s = np.zeros(n)
    for i in range(n):
        s[i] = max(np.abs(a[i,:]))
    for k in range(0,n-1):
        # Row interchange, if needed
        p = np.argmax(np.abs(a[k:n,k])/s[k:n]) + k
        if abs(a[p,k]) < tol: error.err('Matrix is singular')
        if p != k:
            swapRows(b,k,p)
            swapRows(s,k,p)
            swapRows(a,k,p)
        # Elimination
        for i in range(k+1,n):
            if a[i,k] != 0.0:
                lam = a[i,k]/a[k,k]
                a[i,k+1:n] = a[i,k+1:n] - lam*a[k,k+1:n]
                b[i] = b[i] - lam*b[k]
    if abs(a[n-1,n-1]) < tol: error.err('Matrix is singular')
    # Back substitution
    b[n-1] = b[n-1]/a[n-1,n-1]
    for k in range(n-2,-1,-1):
        b[k] = (b[k] - np.dot(a[k,k+1:n],b[k+1:n]))/a[k,k]
    return b
def polyFit(xData,yData,m):
    a = np.zeros((m+1,m+1))
    b = np.zeros(m+1)
    s = np.zeros(2*m+1)
    for i in range(len(xData)):
        temp = yData[i]
        for j in range(m+1):
            b[j] = b[j] + temp
            temp = temp*xData[i]
        temp = 1.0
        for j in range(2*m+1):
            s[j] = s[j] + temp
            temp = temp*xData[i]
    for i in range(m+1):
        for j in range(m+1):
            a[i,j] = s[i+j]
    return gaussPivot(a,b)
degree = 10 # can be any degree
polyFit(xData,yData,degree)    

之后我改用LU分解代码，结果有所改善，但仍未达到要求：

import numpy as np
def swapRows(v,i,j):
    if len(v.shape) == 1:
        v[i],v[j] = v[j],v[i]
    else:
        v[[i,j],:] = v[[j,i],:]
def swapCols(v,i,j):
    v[:,[i,j]] = v[:,[j,i]]
def LUdecomp(a,tol=1.0e-9):
    n = len(a)
    seq = np.array(range(n))
    # Set up scale factors
    s = np.zeros((n))
    for i in range(n):
        s[i] = max(abs(a[i,:]))
    for k in range(0,n-1):
        # Row interchange, if needed
        p = np.argmax(np.abs(a[k:n,k])/s[k:n]) + k
        if abs(a[p,k]) < tol: error.err('Matrix is singular')
        if p != k:
            swapRows(s,k,p)
            swapRows(a,k,p)
            swapRows(seq,k,p)
    # Elimination
        for i in range(k+1,n):
            if a[i,k] != 0.0:
                lam = a[i,k]/a[k,k]
                a[i,k+1:n] = a[i,k+1:n] - lam*a[k,k+1:n]
                a[i,k] = lam
        return a,seq
def LUsolve(a,b,seq):
    n = len(a)
    # Rearrange constant vector; store it in [x]
    x = b.copy()
    for i in range(n):
        x[i] = b[seq[i]]
    # Solution
    for k in range(1,n):
        x[k] = x[k] - np.dot(a[k,0:k],x[0:k])
    x[n-1] = x[n-1]/a[n-1,n-1]
    for k in range(n-2,-1,-1):
        x[k] = (x[k] - np.dot(a[k,k+1:n],x[k+1:n]))/a[k,k]
    return x

编辑2：尝试切比雪夫方法后的问题 #

我尝试了切比雪夫方法，写出如下代码，但遇到两个问题：

import numpy as np

def chebyshev_transform(x, n):
    """
    Transforms x-coordinates to Chebyshev coordinates
    """
    return np.cos(n * np.arccos(x))

def chebyshev_design_matrix(x, n):
    """
    Constructs the Chebyshev design matrix
    """
    x_cheb = chebyshev_transform(x, n)
    T = np.zeros((len(x), n+1))
    T[:,0] = 1
    T[:,1] = x_cheb
    for i in range(2, n+1):
        T[:,i] = 2 * x_cheb * T[:,i-1] - T[:,i-2]
    return T
degree =10
f = lambda x: np.cos(X)
xdata = np.linspace(-1,1,num=100)
ydata = np.array([f(i) for i in xdata])
M = chebyshev_design_matrix(xdata,degree)
D_x ,D_y = np.linalg.qr(M)
D_x, seq = LUdecomp(D_x)
A = LUsolve(D_x,D_y,seq)

1. 我的程序不能使用np.linalg.qr，上面只是测试用，同时我不清楚评论中提到的“慢速方法”公式；
2. 该程序仅支持[-1,1]范围内的x值，是否有归一化方法解决此限制？

感谢帮助。

内容的提问来源于stack exchange，提问作者Alex Itenberg