3.4 曲线拟合的最小二乘法授课

3.4 曲线拟合的最小二乘法 深度系统讲解

最小二乘法是工程、实验、统计分析中处理离散数据的核心工具，本质是离散形式的最佳平方逼近。它解决了带测量误差的实验数据的拟合问题——区别于插值法要求曲线严格穿过所有数据点，最小二乘法追求整体误差的最小化，能有效平滑测量噪声，得到符合数据规律的拟合模型。

---

一、问题背景与核心定义

1.1 问题场景

在科学实验中，我们通常只能得到一组离散的观测数据\(\{(x_i, y_i) \mid i=0,1,\dots,m\}\)，其中\(x_i\)是自变量，\(y_i\)是带测量误差的因变量。我们的目标是：找到一个函数\(y=S^*(x)\)，尽可能贴合所有数据点，而不是严格穿过所有点。

1.2 最小二乘法的数学定义

1. 拟合函数空间：选定一组线性无关的基函数\(\varphi_0(x),\varphi_1(x),\dots,\varphi_n(x)\)，张成拟合空间\(\Phi = \text{span}\{\varphi_0(x),\varphi_1(x),\dots,\varphi_n(x)\}\)，其中\(n < m\)（若\(n=m\)则退化为插值问题）。
   拟合空间中的任意函数可表示为基函数的线性组合：

   \[S(x) = a_0 \varphi_0(x) + a_1 \varphi_1(x) + \dots + a_n \varphi_n(x) \tag{3.38} \]

2. 误差定义：对每个数据点，拟合误差为\(\delta_i = S(x_i) - y_i\)，误差向量\(\boldsymbol{\delta} = (\delta_0,\delta_1,\dots,\delta_m)^T\)。

3. 核心优化目标：找到\(S^*(x) \in \Phi\)，使得误差的平方和最小，即：

   \[\|\boldsymbol{\delta}\|_2^2 = \sum_{i=0}^m \delta_i^2 = \sum_{i=0}^m \left[S^*(x_i) - y_i\right]^2 = \min_{S(x) \in \Phi} \sum_{i=0}^m \left[S(x_i) - y_i\right]^2 \]

4. 加权最小二乘法：更通用的形式是加权平方和最小化，引入权函数\(\omega(x) \geq 0\)：

   \[\|\boldsymbol{\delta}\|_2^2 = \sum_{i=0}^m \omega(x_i) \left[S(x_i) - y_i\right]^2 \tag{3.39} \]

   ✅ 权函数的意义：\(\omega(x_i)\)代表第\(i\)个数据点的可信度/重要性，例如\(\omega(x_i)\)可以是该点的重复观测次数，权重越大，该点在拟合中的优先级越高。

---

二、核心推导：法方程组的建立

最小二乘问题的本质是多元函数的无约束极值问题，我们通过极值必要条件推导出求解系数的线性方程组。

2.1 转化为多元函数极值问题

将\(S(x) = \sum_{j=0}^n a_j \varphi_j(x)\)代入误差平方和，得到关于系数\(a_0,a_1,\dots,a_n\)的多元函数：

\[I(a_0,a_1,\dots,a_n) = \sum_{i=0}^m \omega(x_i) \left[ \sum_{j=0}^n a_j \varphi_j(x_i) - y_i \right]^2 \]

我们的目标是找到\(I\)的极小值点\((a_0^*,a_1^*,\dots,a_n^*)\)。

2.2 极值必要条件：偏导数为0

多元可微函数取得极值的必要条件是：对所有自变量的偏导数为0，即：

\[\frac{\partial I}{\partial a_k} = 0, \quad k=0,1,\dots,n \]

偏导数的详细计算

1. 用链式法则求导：
   \[\frac{\partial I}{\partial a_k} = 2 \sum_{i=0}^m \omega(x_i) \left[ \sum_{j=0}^n a_j \varphi_j(x_i) - y_i \right] \cdot \varphi_k(x_i) = 0 \]

2. 两边除以2，拆分求和项：
   \[\sum_{i=0}^m \omega(x_i) \sum_{j=0}^n a_j \varphi_j(x_i)\varphi_k(x_i) - \sum_{i=0}^m \omega(x_i) y_i \varphi_k(x_i) = 0 \]

3. 交换求和顺序，整理得：
   \[\sum_{j=0}^n \left( \sum_{i=0}^m \omega(x_i) \varphi_j(x_i)\varphi_k(x_i) \right) a_j = \sum_{i=0}^m \omega(x_i) y_i \varphi_k(x_i) \]

2.3 离散内积与法方程组

我们引入离散内积的记号，和连续函数的积分内积完全对应：

• 基函数的内积：\((\varphi_j, \varphi_k) = \sum_{i=0}^m \omega(x_i) \varphi_j(x_i)\varphi_k(x_i)\)
• 数据与基函数的内积：\((y, \varphi_k) = \sum_{i=0}^m \omega(x_i) y_i \varphi_k(x_i) \triangleq d_k\)

代入后得到核心线性方程组：

\[\sum_{j=0}^n (\varphi_k, \varphi_j) a_j = d_k, \quad k=0,1,\dots,n \tag{3.40} \]

该方程组称为法方程组（正规方程组），其矩阵形式为：

\[\boldsymbol{G}\boldsymbol{a} = \boldsymbol{d} \]

其中：

• 未知系数向量\(\boldsymbol{a} = (a_0,a_1,\dots,a_n)^T\)，右端向量\(\boldsymbol{d} = (d_0,d_1,\dots,d_n)^T\)
• 系数矩阵\(\boldsymbol{G}\)称为格拉姆（Gram）矩阵，形式为：
  \[\boldsymbol{G} = \begin{pmatrix} (\varphi_0,\varphi_0) & (\varphi_0,\varphi_1) & \dots & (\varphi_0,\varphi_n) \\ (\varphi_1,\varphi_0) & (\varphi_1,\varphi_1) & \dots & (\varphi_1,\varphi_n) \\ \vdots & \vdots & & \vdots \\ (\varphi_n,\varphi_0) & (\varphi_n,\varphi_1) & \dots & (\varphi_n,\varphi_n) \end{pmatrix} \tag{3.41} \]

✅ 核心性质：格拉姆矩阵是对称矩阵，即\((\varphi_j,\varphi_k)=(\varphi_k,\varphi_j)\)，因此\(\boldsymbol{G}^T=\boldsymbol{G}\)。

---

三、解的存在唯一性：哈尔（Haar）条件

法方程组有唯一解的前提是格拉姆矩阵\(\boldsymbol{G}\)非奇异（可逆），这里需要注意一个关键问题：基函数在连续区间上线性无关，不能推出离散点集上的格拉姆矩阵非奇异。

3.1 反例说明

例如\(\varphi_0(x)=\sin x\)，\(\varphi_1(x)=\sin2x\)，在\([0,2\pi]\)上连续线性无关；但取离散点\(x_k=k\pi\)（\(k=0,1,2\)），则\(\varphi_0(x_k)=\varphi_1(x_k)=0\)，此时格拉姆矩阵的所有元素均为0，矩阵奇异，法方程组无唯一解。

3.2 哈尔条件（解唯一的充分条件）

定义3.7 设\(\varphi_0(x),\varphi_1(x),\dots,\varphi_n(x) \in C[a,b]\)，若其任意非零线性组合在点集\(\{x_0,x_1,\dots,x_m\}\)（\(m\geq n\)）上至多只有\(n\)个不同的零点，则称该基函数在点集上满足哈尔条件。

✅ 经典特例：多项式基\(\{1,x,x^2,\dots,x^n\}\)在任意\(m\geq n\)个互异点上满足哈尔条件，因为n次非零多项式至多有n个不同的实根。

3.3 存在唯一性定理

若基函数\(\{\varphi_0,\dots,\varphi_n\}\)在点集\(\{x_0,\dots,x_m\}\)上满足哈尔条件，则格拉姆矩阵\(\boldsymbol{G}\)非奇异，法方程组存在唯一解\(\boldsymbol{a}^*=(a_0^*,a_1^*,\dots,a_n^*)\)，对应的函数：

\[S^*(x) = a_0^* \varphi_0(x) + a_1^* \varphi_1(x) + \dots + a_n^* \varphi_n(x) \]

就是最小二乘问题的唯一解，且对任意\(S(x)\in\Phi\)，都有：

\[\sum_{i=0}^m \omega(x_i)\left[S^*(x_i)-y_i\right]^2 \leq \sum_{i=0}^m \omega(x_i)\left[S(x_i)-y_i\right]^2 \]

---

四、核心应用场景与实用技巧

4.1 最常用场景：多项式拟合

当基函数取\(\varphi_k(x)=x^k\)（\(k=0,1,\dots,n\)）时，拟合函数为n次多项式\(S(x)=a_0+a_1x+\dots+a_nx^n\)，称为多项式拟合。

• 当\(n=1\)时，就是线性拟合（直线拟合），是工程中最常用的拟合形式；
• 注意：当\(n\geq3\)时，格拉姆矩阵退化为希尔伯特矩阵，高度病态，数值计算中误差会被急剧放大，因此高次多项式拟合需慎用，建议优先使用正交多项式拟合。

4.2 非线性模型的线性化技巧

很多实际问题的模型表面上是非线性的，但可以通过变量变换转化为线性模型，从而用最小二乘法求解。

经典例子：指数模型拟合

模型形式为\(S(x)=ae^{bx}\)，其中\(a,b\)为待求参数。

• 对两边取自然对数：\(\ln S(x) = \ln a + bx\)
• 令\(Y = \ln S(x)\)，\(A = \ln a\)，则模型转化为线性形式\(Y = A + bx\)
• 对变换后的数据\((x_i, \ln y_i)\)做线性拟合，得到\(A\)和\(b\)，再通过\(a=e^A\)还原原参数。

---

五、例题详解（例3.10）

问题描述

已知实验数据如下表，求加权最小二乘拟合曲线。

\(x_i\)	1	2	3	4	5
\(y_i\)	4	4.5	6	8	8.5
\(\omega_i\)	2	1	3	1	1

求解步骤

步骤1：确定拟合模型

绘制散点图，可见数据点近似分布在一条直线附近，因此选择线性拟合模型：

\[S_1(x) = a_0 + a_1 x \]

基函数\(\varphi_0(x)=1\)，\(\varphi_1(x)=x\)，\(m=4\)（5个点，索引0-4），\(n=1\)。

步骤2：计算离散内积

根据加权离散内积公式，逐个计算：

1. \((\varphi_0,\varphi_0) = \sum_{i=0}^4 \omega_i = 2+1+3+1+1 = 8\)
2. \((\varphi_0,\varphi_1) = (\varphi_1,\varphi_0) = \sum_{i=0}^4 \omega_i x_i = 2\times1 + 1\times2 + 3\times3 + 1\times4 + 1\times5 = 22\)
3. \((\varphi_1,\varphi_1) = \sum_{i=0}^4 \omega_i x_i^2 = 2\times1^2 + 1\times2^2 + 3\times3^2 + 1\times4^2 + 1\times5^2 = 74\)
4. \((\varphi_0,y) = \sum_{i=0}^4 \omega_i y_i = 2\times4 + 1\times4.5 + 3\times6 + 1\times8 + 1\times8.5 = 47\)
5. \((\varphi_1,y) = \sum_{i=0}^4 \omega_i x_i y_i = 2\times1\times4 + 1\times2\times4.5 + 3\times3\times6 + 1\times4\times8 + 1\times5\times8.5 = 145.5\)

步骤3：建立并解法方程组

代入法方程组，得到：

\[\begin{cases} 8a_0 + 22a_1 = 47 \\ 22a_0 + 74a_1 = 145.5 \end{cases} \]

用消元法求解：

1. 第一个方程两边乘11，第二个方程两边乘4，得：
   \[\begin{cases} 88a_0 + 242a_1 = 517 \\ 88a_0 + 296a_1 = 582 \end{cases} \]

2. 两式相减，得\(54a_1 = 65\)，解得\(a_1 = \frac{65}{54} \approx 1.2037\)；
3. 代入第一个方程，解得\(a_0 = \frac{47 - 22\times1.2037}{8} \approx 2.5648\)。

步骤4：得到拟合曲线

最终的加权最小二乘拟合直线为：

\[S_1^*(x) = 2.5648 + 1.2037x \]

---

六、核心总结与高频易错点

6.1 内容总结

1. 最小二乘法的本质是离散数据的最佳平方逼近，核心是通过最小化误差平方和，将拟合问题转化为线性方程组（法方程组）的求解；
2. 法方程组的系数矩阵是格拉姆矩阵，基函数满足哈尔条件时，方程组有唯一解；
3. 线性拟合是最常用的形式，非线性模型可通过变量变换转化为线性模型求解；
4. 高次多项式拟合会出现病态问题，工程中优先使用低次拟合或正交多项式拟合。

6.2 高频易错点提醒

1. 权函数的意义不能忽略：权重代表数据的可信度，不等精度的观测数据必须使用加权最小二乘法，否则会导致拟合结果偏向低精度数据；
2. 连续线性无关≠离散可逆：基函数在连续区间线性无关，不能保证离散点集上的格拉姆矩阵非奇异，必须满足哈尔条件才能保证解唯一；
3. 高次拟合的病态问题：n≥3的多项式拟合，格拉姆矩阵是希尔伯特矩阵，条件数极大，常规数值方法求解会出现严重的误差失真；
4. 非线性线性化的误差本质：指数模型等非线性模型线性化后，拟合的是变换后的变量（如\(\ln y\)），最小化的是变换后的误差，不是原变量的误差，对极端值敏感，需谨慎使用；
5. 拟合阶数的选择：不是阶数越高拟合效果越好，过高的阶数会出现“过拟合”——在已知数据点上误差很小，但对新数据的预测能力极差。

例3.11 最小二乘拟合（线性模型+可线性化的指数模型）完整详解

这个例题是最小二乘法的经典应用，完整演示了从模型选择、线性拟合、残差分析、模型修正到非线性模型线性化拟合的全流程，是工程中处理实验数据的标准流程，我们逐步骤拆解所有计算细节和逻辑。

---

一、问题与数据明确

1.1 已知数据

给定5组实验观测数据如下：

序号\(i\)	0	1	2	3	4
\(x_i\)	1.00	1.25	1.50	1.75	2.00
\(y_i\)	5.10	5.79	6.53	7.45	8.46

1.2 拟合目标

找到贴合数据的函数模型，分为两步：

1. 先尝试线性模型拟合；
2. 分析拟合缺陷后，修正为指数模型，通过变量变换转化为线性模型求解。

---

二、第一步：线性模型最小二乘拟合

2.1 模型选择

绘制\((x_i,y_i)\)的散点图（图3-10(a)），可见数据点近似分布在一条直线附近，因此先选择线性拟合模型：

\[y = A + bx \]

对应基函数\(\varphi_0(x)=1\)，\(\varphi_1(x)=x\)，权函数\(\omega(x)\equiv1\)（等精度观测），拟合空间\(\Phi=\text{span}\{1,x\}\)。

2.2 计算离散内积

根据最小二乘法的离散内积公式\((\varphi_j,\varphi_k)=\sum_{i=0}^4 \omega_i \varphi_j(x_i)\varphi_k(x_i)\)，\((y,\varphi_k)=\sum_{i=0}^4 \omega_i y_i \varphi_k(x_i)\)，逐个计算：

1. \((\varphi_0,\varphi_0) = \sum_{i=0}^4 1\times1\times1 = 5\)（共5个数据点）
2. \((\varphi_0,\varphi_1) = (\varphi_1,\varphi_0) = \sum_{i=0}^4 x_i = 1.00+1.25+1.50+1.75+2.00 = 7.50\)
3. \((\varphi_1,\varphi_1) = \sum_{i=0}^4 x_i^2 = 1^2 + 1.25^2 + 1.5^2 + 1.75^2 + 2^2 = 1+1.5625+2.25+3.0625+4 = 11.875\)
4. \((\varphi_0,y) = \sum_{i=0}^4 y_i = 5.10+5.79+6.53+7.45+8.46 = 33.33\)
5. \((\varphi_1,y) = \sum_{i=0}^4 x_i y_i\)，逐项计算：
   • \(1.00\times5.10=5.10\)，\(1.25\times5.79=7.2375\)，\(1.50\times6.53=9.795\)
   • \(1.75\times7.45=13.0375\)，\(2.00\times8.46=16.92\)
   • 求和得：\(5.10+7.2375+9.795+13.0375+16.92 = 52.09\)

2.3 建立并解法方程组

将内积代入法方程组\(\boldsymbol{G}\boldsymbol{a}=\boldsymbol{d}\)，得到：

\[\begin{cases} 5A + 7.50b = 33.33 \\ 7.50A + 11.875b = 52.09 \end{cases} \]

用消元法求解：

1. 第一个方程两边乘以1.5，得：\(7.5A + 11.25b = 49.995\)
2. 用第二个方程减去上式，消去\(A\)：
   \[(7.50A + 11.875b) - (7.5A + 11.25b) = 52.09 - 49.995 \]
   \[0.625b = 2.095 \implies b = \frac{2.095}{0.625} = 3.3520 \]

3. 将\(b=3.3520\)代入第一个方程，解得：
   \[5A = 33.33 - 7.50\times3.3520 = 8.19 \implies A = 1.6380 \]

2.4 线性拟合结果与缺陷分析

最终得到线性拟合直线：

\[y = 1.6380 + 3.3520x \]

拟合缺陷分析：
从拟合效果图（图3-10(b)）可见，中间3个数据点全部在拟合直线的同一侧，两个端点在直线的另一侧，误差呈现系统性偏差，而非随机波动。这说明线性模型不符合数据的真实规律，需要更换更合适的模型。

---

三、第二步：指数模型的最小二乘拟合

3.1 模型修正与线性化

观察数据的增长趋势，\(y\)随\(x\)的增长呈现加速上升的特征，符合指数增长规律，因此选择指数模型：

\[y = a e^{bx} \]

其中\(a,b\)为待求参数，这是一个非线性模型，无法直接用最小二乘法求解，因此通过对数变换将其线性化：

对模型两边取自然对数：

\[\ln y = \ln a + bx \]

令\(\bar{y} = \ln y\)，\(A = \ln a\)，模型转化为标准线性形式：

\[\bar{y} = A + bx \]

此时仍用基函数\(\varphi_0(x)=1\)，\(\varphi_1(x)=x\)，只需将原数据\(y_i\)替换为\(\bar{y}_i=\ln y_i\)，即可用最小二乘法求解。

3.2 变换后的数据与内积计算

先计算每个\(y_i\)对应的\(\bar{y}_i=\ln y_i\)：

\(i\)	0	1	2	3	4
\(y_i\)	5.10	5.79	6.53	7.45	8.46
\(\bar{y}_i=\ln y_i\)	1.6292	1.7561	1.8764	2.0082	2.1351

计算变换后的内积（\(x_i\)不变，因此系数矩阵和线性拟合完全一致，仅需重新计算右端项）：

1. \((\varphi_0,\bar{y}) = \sum_{i=0}^4 \bar{y}_i = 1.6292+1.7561+1.8764+2.0082+2.1351 = 9.404\)
2. \((\varphi_1,\bar{y}) = \sum_{i=0}^4 x_i \bar{y}_i\)，逐项计算：
   • \(1.00\times1.6292=1.6292\)，\(1.25\times1.7561=2.1951\)，\(1.50\times1.8764=2.8146\)
   • \(1.75\times2.0082=3.5144\)，\(2.00\times2.1351=4.2702\)
   • 求和得：\(1.6292+2.1951+2.8146+3.5144+4.2702 = 14.422\)

3.3 建立并解法方程组

法方程组为：

\[\begin{cases} 5A + 7.50b = 9.404 \\ 7.50A + 11.875b = 14.422 \end{cases} \]

同样用消元法求解：

1. 第一个方程乘以1.5得：\(7.5A + 11.25b = 14.106\)
2. 第二个方程减上式，消去\(A\)：
   \[0.625b = 14.422 - 14.106 = 0.316 \implies b = \frac{0.316}{0.625} = 0.5056 \]

3. 代入第一个方程解得：
   \[5A = 9.404 - 7.50\times0.5056 = 5.612 \implies A = 1.1224 \]

3.4 还原指数模型

由\(A = \ln a\)，还原原参数\(a\)：

\[a = e^A = e^{1.1224} \approx 3.0722 \]

最终得到指数拟合模型：

\[y = 3.0722 e^{0.5056x} \]

3.5 拟合效果验证

绘制\((x_i,\ln y_i)\)的散点图（图3-10(c)），可见数据点完美分布在拟合直线附近，无系统性偏差；还原后的指数模型曲线（图3-10(d)）完全贴合原数据，拟合效果远优于线性模型。

---

四、核心总结与注意事项

4.1 拟合流程总结

1. 模型初选：通过散点图选择初始模型；
2. 线性拟合：用最小二乘法求解，得到初始拟合结果；
3. 残差分析：通过误差分布判断模型是否合适，若存在系统性偏差，更换模型；
4. 非线性模型线性化：对指数、幂函数等可线性化的模型，通过变量变换转化为线性模型；
5. 求解与还原：对变换后的数据做线性拟合，再还原为原模型。

4.2 关键注意事项

1. 线性化的误差本质：对数变换后，最小二乘最小化的是\(\ln y\)的误差，而非原\(y\)的误差，适合原数据的误差为乘性误差的场景；若为加性误差，需用非线性最小二乘直接求解。
2. 变换的前提：对数变换要求所有\(y_i>0\)，若存在非正数据，需先做平移变换。
3. 模型选择的核心：拟合的核心是匹配数据的真实规律，而非追求高次、复杂模型，残差的随机分布是模型合适的核心标志。
4. 法方程组的复用：线性化后，若自变量\(x_i\)不变，法方程组的系数矩阵完全不变，仅需重新计算右端项，大幅减少计算量。

3.4.2 用正交多项式作最小二乘拟合 深度系统讲解

这一节的内容是多项式最小二乘拟合的终极实用解决方案，彻底解决了普通幂函数基拟合中法方程组病态、高次拟合数值失真的核心痛点，是工程中高次多项式拟合的标准、高效、稳定方法。

我们从问题根源、核心原理、递推构造、正交性证明、拟合流程、方法优势六个维度，完整拆解这部分内容。

---

一、问题根源：普通最小二乘的致命缺陷

用普通幂函数基\(\{1,x,x^2,\dots,x^n\}\)做最小二乘拟合时，法方程组的系数矩阵（格拉姆矩阵）是希尔伯特矩阵，它有一个致命缺陷：高度病态。

• 当拟合次数\(n\geq3\)时，矩阵条件数急剧增大，数值计算中微小的舍入误差会被放大成千上万倍，导致拟合系数完全失真；
• 当\(n\geq5\)时，常规计算方法得到的结果已经完全不可信，无法实现高次多项式拟合。

而解决这个问题的核心方法，就是用离散正交多项式作为拟合基函数——正交基下的法方程组是对角矩阵，条件数为1，完全消除病态问题，同时将复杂的方程组求解简化为独立的系数计算。

---

二、正交基下的最小二乘核心公式

2.1 离散正交的定义

给定离散节点\(\{x_0,x_1,\dots,x_m\}\)和权函数\(\omega(x_i)\geq0\)，若函数族\(\{\varphi_0(x),\varphi_1(x),\dots,\varphi_n(x)\}\)满足带权离散正交性：

\[(\varphi_j,\varphi_k) = \sum_{i=0}^m \omega(x_i) \varphi_j(x_i)\varphi_k(x_i) = \begin{cases} 0, & j \neq k \\ A_k > 0, & j = k \end{cases} \]

即不同基函数的离散内积为0，自身内积为正数。

2.2 正交基下的最小二乘解

当基函数满足离散正交性时，法方程组\(\boldsymbol{G}\boldsymbol{a}=\boldsymbol{d}\)完全解耦，每个系数可以独立计算，无需解联立方程组：

\[a_k^* = \frac{(y,\varphi_k)}{(\varphi_k,\varphi_k)} = \frac{\sum_{i=0}^m \omega(x_i) y_i \varphi_k(x_i)}{\sum_{i=0}^m \omega(x_i) \varphi_k^2(x_i)}, \quad k=0,1,\dots,n \]

2.3 平方误差公式

拟合的平方误差可以直接通过系数计算，无需逐个点计算误差：

\[\|\boldsymbol{\delta}\|_2^2 = \|y\|_2^2 - \sum_{k=0}^n A_k (a_k^*)^2 \]

其中\(\|y\|_2^2 = \sum_{i=0}^m \omega(x_i) y_i^2\)，\(A_k=(\varphi_k,\varphi_k)\)。

✅ 核心优势：

1. 完全避免解线性方程组，彻底消除病态问题，数值稳定性拉满；
2. 系数独立计算，增加拟合次数时，低次项的系数无需重新计算；
3. 误差计算简单，可在计算过程中动态监控拟合精度。

---

三、离散正交多项式的递推构造

正交多项式不是固定的，而是与给定的节点\(\{x_i\}\)和权\(\omega(x_i)\)严格绑定的，因此我们需要一套通用的方法，针对任意节点和权函数，构造对应的离散正交多项式。

我们采用三项递推公式构造首项系数为1的离散正交多项式\(\{p_k(x)\}\)（首一：最高次项系数为1，保证递推的一致性）。

3.1 递推公式

\[\begin{cases} p_0(x) = 1, \\ p_1(x) = (x - \alpha_1) p_0(x), \\ p_{k+1}(x) = (x - \alpha_{k+1}) p_k(x) - \beta_k p_{k-1}(x), \quad k=1,2,\dots,n-1 \end{cases} \tag{3.42} \]

3.2 递推参数计算公式

递推中的参数\(\alpha_{k+1}\)和\(\beta_k\)由离散内积直接计算：

\[\begin{cases} \alpha_{k+1} = \frac{\sum_{i=0}^m \omega(x_i) x_i p_k^2(x_i)}{\sum_{i=0}^m \omega(x_i) p_k^2(x_i)} = \frac{(x p_k, p_k)}{(p_k, p_k)}, \quad k=0,1,\dots,n-1 \\ \beta_k = \frac{\sum_{i=0}^m \omega(x_i) p_k^2(x_i)}{\sum_{i=0}^m \omega(x_i) p_{k-1}^2(x_i)} = \frac{(p_k, p_k)}{(p_{k-1}, p_{k-1})}, \quad k=1,2,\dots,n-1 \end{cases} \tag{3.43} \]

✅ 公式解读：

• \(\alpha_{k+1}\)：本质是\(x p_k(x)\)在\(p_k(x)\)上的投影系数，保证\(p_{k+1}(x)\)与\(p_k(x)\)正交；
• \(\beta_k\)：本质是\(p_k(x)\)与\(p_{k-1}(x)\)的范数比值，保证\(p_{k+1}(x)\)与\(p_{k-1}(x)\)正交；
• 三项递推是构造正交多项式的最优方式，仅需前两项即可构造下一项，计算量极小，编程实现极其简单。

---

四、正交性的严谨证明（数学归纳法）

我们用数学归纳法严格证明：由上述递推公式构造的\(\{p_k(x)\}\)，是关于节点\(\{x_i\}\)带权\(\omega(x_i)\)的正交多项式族，即对任意\(0\leq s < l \leq k\)，有\((p_s,p_l)=0\)。

步骤1：基例验证（k=1）

证明\(p_0\)与\(p_1\)正交：

\[\begin{align*} (p_0,p_1) &= (p_0, (x-\alpha_1)p_0) \\ &= (p_0,xp_0) - \alpha_1 (p_0,p_0) \end{align*} \]

代入\(\alpha_1 = \frac{(xp_0,p_0)}{(p_0,p_0)}\)，得：

\[(p_0,p_1) = (p_0,xp_0) - \frac{(xp_0,p_0)}{(p_0,p_0)} \cdot (p_0,p_0) = 0 \]

基例成立，\(p_0\)与\(p_1\)正交。

步骤2：归纳假设

假设对所有\(0\leq s < l \leq k\)，都有\((p_s,p_l)=0\)，即前\(k+1\)个多项式\(p_0,p_1,\dots,p_k\)两两正交。

步骤3：归纳递推（证明\(p_{k+1}\)与\(p_0,p_1,\dots,p_k\)全部正交）

对任意\(0\leq s \leq k\)，计算内积：

\[(p_{k+1},p_s) = \left( (x-\alpha_{k+1})p_k - \beta_k p_{k-1}, p_s \right) = (xp_k,p_s) - \alpha_{k+1}(p_k,p_s) - \beta_k(p_{k-1},p_s) \]

我们分三种情况分别证明该内积为0：

情况1：\(0\leq s \leq k-2\)

由归纳假设，\((p_k,p_s)=0\)，\((p_{k-1},p_s)=0\)。
同时，\(x p_s(x)\)是首一的\(s+1\)次多项式，而\(s+1 \leq k-1\)，因此\(x p_s(x)\)可以表示为\(p_0,p_1,\dots,p_{s+1}\)的线性组合，由归纳假设，它与\(p_k\)正交，即\((xp_k,p_s)=(p_k,xp_s)=0\)。
因此\((p_{k+1},p_s)=0\)。

情况2：\(s = k-1\)

由归纳假设，\((p_k,p_{k-1})=0\)。
同时，\(x p_{k-1}(x)\)是首一的\(k\)次多项式，可表示为\(p_k(x) + \sum_{j=0}^{k-1} c_j p_j(x)\)，因此\((xp_k,p_{k-1})=(p_k,xp_{k-1})=(p_k,p_k)\)。
代入\(\beta_k = \frac{(p_k,p_k)}{(p_{k-1},p_{k-1})}\)，得：

\[(p_{k+1},p_{k-1}) = (p_k,p_k) - 0 - \frac{(p_k,p_k)}{(p_{k-1},p_{k-1})} \cdot (p_{k-1},p_{k-1}) = 0 \]

情况3：\(s = k\)

由归纳假设，\((p_{k-1},p_k)=0\)。
代入\(\alpha_{k+1} = \frac{(xp_k,p_k)}{(p_k,p_k)}\)，得：

\[(p_{k+1},p_k) = (xp_k,p_k) - \frac{(xp_k,p_k)}{(p_k,p_k)} \cdot (p_k,p_k) - 0 = 0 \]

步骤4：结论

由数学归纳法，对所有\(0\leq s < l \leq n\)，都有\((p_s,p_l)=0\)，即递推构造的\(\{p_k(x)\}\)是离散正交多项式族。

---

五、用正交多项式做最小二乘拟合的完整流程

这个方法的核心优势是边构造正交多项式，边计算拟合系数，边累加得到最终拟合多项式，无需先构造所有多项式再计算，全程仅需循环递推，编程实现极其简单。

完整步骤

1. 初始化（k=0）

   • 构造\(p_0(x)=1\)；
   • 计算内积：\((p_0,p_0)=\sum_{i=0}^m \omega(x_i) \cdot 1^2\)，\((y,p_0)=\sum_{i=0}^m \omega(x_i) y_i \cdot 1\)；
   • 计算系数：\(a_0^* = \frac{(y,p_0)}{(p_0,p_0)}\)；
   • 初始拟合多项式：\(S(x) = a_0^* p_0(x) = a_0^*\)。

2. 构造一次多项式（k=1）

   • 计算\(\alpha_1 = \frac{(x p_0,p_0)}{(p_0,p_0)}\)；
   • 构造\(p_1(x) = (x - \alpha_1) p_0(x) = x - \alpha_1\)；
   • 计算内积：\((p_1,p_1)=\sum_{i=0}^m \omega(x_i) p_1^2(x_i)\)，\((y,p_1)=\sum_{i=0}^m \omega(x_i) y_i p_1(x_i)\)；
   • 计算系数：\(a_1^* = \frac{(y,p_1)}{(p_1,p_1)}\)；
   • 累加得到一次拟合多项式：\(S(x) = a_0^* p_0(x) + a_1^* p_1(x)\)。

3. 递推构造高次多项式（k≥2）
   对\(k=1,2,\dots,n-1\)，循环执行：

   • 计算\(\alpha_{k+1} = \frac{(x p_k,p_k)}{(p_k,p_k)}\)，\(\beta_k = \frac{(p_k,p_k)}{(p_{k-1},p_{k-1})}\)；
   • 构造\(p_{k+1}(x) = (x - \alpha_{k+1}) p_k(x) - \beta_k p_{k-1}(x)\)；
   • 计算内积：\((p_{k+1},p_{k+1})\)，\((y,p_{k+1})\)；
   • 计算系数：\(a_{k+1}^* = \frac{(y,p_{k+1})}{(p_{k+1},p_{k+1})}\)；
   • 累加得到\(k+1\)次拟合多项式：\(S(x) = S(x) + a_{k+1}^* p_{k+1}(x)\)；
   • 可同步计算拟合误差，若误差满足要求，可提前终止循环。

4. 最终结果
   循环结束后，得到的\(S(x) = \sum_{k=0}^n a_k^* p_k(x)\)就是n次最小二乘拟合多项式，可展开为标准幂函数形式\(S(x)=c_0 + c_1 x + \dots + c_n x^n\)。

---

六、方法的核心优势与工程价值

1. 彻底消除数值病态
   全程无需解线性方程组，完全避免了希尔伯特矩阵的病态问题，哪怕拟合次数\(n\)取到几十、上百，计算结果依然稳定可靠，彻底解决了高次多项式拟合的数值失真问题。

2. 计算效率极高
   仅需一次循环即可完成正交多项式构造、系数计算、拟合多项式累加，时间复杂度为\(O(nm)\)（n为拟合次数，m为数据点个数），远低于求解线性方程组的\(O(n^3)\)复杂度。

3. 灵活性极强

   • 增加拟合次数时，之前的所有计算完全复用，仅需多执行一次循环，无需重新计算；
   • 可在计算过程中动态监控误差，随时终止循环，无需事先固定拟合次数，完美适配“按需拟合”的工程场景。

4. 编程实现极简
   递推公式结构固定，仅需简单的循环即可实现，有成熟的通用程序框架，是目前工程中多项式曲线拟合的首选标准方法。

---

七、关键注意事项

1. 正交多项式的定制性
   递推构造的正交多项式是与当前节点\(\{x_i\}\)和权\(\omega(x_i)\)严格绑定的，节点或权函数改变，对应的正交多项式也会改变，不能直接套用连续区间的勒让德、切比雪夫多项式。

2. 首一性的意义
   递推公式构造的是首一多项式，保证了递推的一致性和正交性，若修改首项系数，需要同步调整递推参数，不建议修改。

3. 权函数的作用
   权函数\(\omega(x_i)\)可用于区分数据点的可信度，权重越大，该点在拟合中的优先级越高，比如\(\omega(x_i)\)可取该点的重复观测次数、测量精度的倒数等。