记录一下工作时遇到的拟合问题，将两个数据的关系建模为最小二乘的模型：
\[y = a_0 + a_1 x + a_2 x^2 + a_3 x^3 + a_4 x^4 \]使用了python里面的numpy.linalg.lstsq函数进行拟合，以下是一个简单的示例
import numpy as np
import matplotlib.pyplot as plt
# 样本数据点
x = np.array([-10, -8, -5, -3, 0, 2, 5, 7, 9, 12])
y = np.array([1200, 800, 200, 100, 50, 80, 300, 600, 1000, 1800])
# 构建设计矩阵
X = np.vstack([x**n for n in range(5)]).T # 包含 x 的 0 次到 4 次方
# 求解最小二乘问题
coefficients, residuals, rank, singular_values = np.linalg.lstsq(X, y, rcond=None)
# 打印系数
print("拟合多项式的系数为：", coefficients)
# 计算拟合的 y 值
y_fit = X @ coefficients
# 绘制结果
plt.scatter(x, y, color='blue', label='原始数据')
plt.plot(x, y_fit, color='red', label='拟合多项式')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.title('最小二乘法多项式拟合')
plt.show()
调库拟合感觉很抽象，没办法深入理解到底做了什么事情。
于是研究了一下拟合的时候具体做了哪些步骤：
在拟合多项式 $y = a_0 + a_1 x + a_2 x^2 + a_3 x^3 + a_4 x^4 $的系数时，我们需要解决一个线性方程组，使得预测值与实际值之间的误差平方和最小。这就是典型的最小二乘问题。
numpy.linalg.lstsq 的内部操作步骤：
-
构建设计矩阵 $ A $：
-
我们将多项式的每一项视为一个特征，构建一个设计矩阵$ A $，其大小为 $ m \times n $（其中 $ m $是数据点的数量， $ n $ 是多项式系数的数量）。
-
矩阵 $ A $ 的每一列对应于 $ x $ 的不同次幂：
\[A = \begin{bmatrix} 1 & x_1 & x_1^2 & x_1^3 & x_1^4 \\ 1 & x_2 & x_2^2 & x_2^3 & x_2^4 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & x_m & x_m^2 & x_m^3 & x_m^4 \\ \end{bmatrix} \]
-
-
设定观测向量 $ b $：
-
向量 $ b $ 包含了对应的 $ y $ 值：
\[b = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_m \\ \end{bmatrix} \]
-
-
定义最小二乘目标：
-
我们的目标是找到系数向量 $ x $（即 $ a_0, a_1, a_2, a_3, a_4 $），使得残差的平方和最小：
\[\min_x \| Ax - b \|_2^2 \]
-
-
计算矩阵的奇异值分解（SVD）：
-
内部，
\[A = U \Sigma V^T \]lstsq函数使用 奇异值分解（Singular Value Decomposition, SVD） 对矩阵 $ A $ 进行分解：- $ U $ 是 $ m \times m $ 的正交矩阵。
- $ \Sigma $ 是 $ m \times n $ 的对角矩阵，对角线上的元素为 $ A $ 的奇异值。
- $ V^T $ 是 $ n \times n $ 的正交矩阵的转置。
-
-
计算矩阵 $ A $ 的伪逆 $ A^+ $：
-
基于 SVD 分解，计算 $ A $ 的 Moore-Penrose 伪逆：
\[A^+ = V \Sigma^+ U^T \]- $ \Sigma^+ $ 是将 $ \Sigma $ 中非零元素取倒数后转置而得。
-
-
求解最小二乘解 $ x $：
-
使用伪逆计算解：
\[x = A^+ b \]- 这给出了使得 $ | Ax - b |_2 $ 最小的系数向量 $ x $。
-
-
计算残差 $ r $：
-
计算残差向量：
\[r = b - Ax \] -
残差的平方和为：
\[\text{residuals} = \| r \|_2^2 \]
-
-
确定矩阵 $ A $ 的秩：
- 通过奇异值判断 $ A $ 的 有效秩（rank），这有助于了解问题的可解性和解的稳定性。
-
返回结果：
lstsq函数返回以下内容：- 系数向量 $ x $：拟合多项式的系数。
- 残差数组：残差的平方和（若方程组为超定，即方程数大于未知数）。
- 秩 $ rank $：矩阵 $ A $ 的秩。
- 奇异值数组：矩阵 $ A $ 的奇异值。
总结：
numpy.linalg.lstsq函数通过对设计矩阵 ( A ) 进行奇异值分解，计算其伪逆，然后求解最小二乘问题，得到使预测值与实际值误差平方和最小的系数。- 使用 SVD 有以下优点：
- 数值稳定性：SVD 对于病态矩阵（条件数很大）能够提供稳定的求解。
- 处理欠定和超定方程组：即使矩阵 ( A ) 不满秩，SVD 也能找到范数最小的最小二乘解。
附加说明：
-
为什么不是直接求解正常方程 $ A^T A x = A^T b $？
- 直接求解可能导致数值不稳定，尤其是当 $ A^T A $ 的条件数很大时。
- 使用 SVD 可以避免这些问题，提供更可靠的解。
-
计算复杂度：
- SVD 的计算复杂度较高，但对于中小规模的问题（如多项式拟合），现代计算机可以快速完成计算。
举例：
假设我们有以下数据点：
import numpy as np
# 样本数据点
x = np.array([-10, -8, -5, -3, 0, 2, 5, 7, 9, 12])
y = np.array([1200, 800, 200, 100, 50, 80, 300, 600, 1000, 1800])
# 构建设计矩阵 A
A = np.vstack([x**0, x**1, x**2, x**3, x**4]).T
# 使用最小二乘法求解
coefficients, residuals, rank, singular_values = np.linalg.lstsq(A, y, rcond=None)
print("拟合多项式的系数为：", coefficients)
print("残差的平方和为：", residuals)
print("矩阵 A 的秩为：", rank)
print("矩阵 A 的奇异值为：", singular_values)
详细步骤解读：
-
构建设计矩阵 $ A $：
- 每一列对应 $ x $ 的 0 到 4 次幂，共 5 列。
- $ A $ 的大小为 $ 10 \times 5 $。
-
使用
np.linalg.lstsq求解：- 函数内部对 $ A $ 进行 SVD 分解。
- 计算 $ A $ 的伪逆 $ A^+ $。
- 求解系数向量 $ x = A^+ b $。
-
解的解释：
coefficients：拟合的多项式系数 $ [a_0, a_1, a_2, a_3, a_4] $。residuals：预测值与实际 $ y $ 值之间的误差平方和。rank：矩阵 $ A $ 的秩，表示特征的线性独立数量。singular_values：用于判断矩阵是否存在病态或奇异情况。