计算方法笔记
基本概念
计算方法
数值计算方法是科学计算的核心，研究对象是那些在理论上有解，而无求解公式或计算量过大必须借助计算机实现的数学问题。
简单而言就是容忍一定的误差，通过一个改进的在计算机上迭代的计算方法获取原解的近似解。
而数值计算方法主要有两大方向：
- 数值算法的构造：针对数值问题研究可在计算机上执行且行之有效的计算公式（计算方法）
- 算法的理论分析：主要包括误差分析（数值问题的性态，数值方法的截断误差，舍入误差，稳定性，收敛性等）和复杂性分析
也就是说，我们希望找到一种合理的算法能够降低其复杂度的同时保证结果的合理性。
误差
首先明确，数值方法得到的都是近似解而非真解。
那么误差的定义就是真解与近似解之间的偏差，即 \[ e=x-x^* \] 其中 \(x^*\) 是近似解的表示符号。
误差分为四类：
- 模型误差（model error）
- 观测误差（observational error）
- 截断误差（truncation error）
- 舍入误差（round-off error）
前两者是客观存在的误差，这是原问题模型本身或测量方法所存在的误差；
后两者是计算方法所引起的误差，这是我们所需要关注的。
由简化问题所引起的误差称为截断误差或方法误差，例如我们通常会将三角函数进行泰勒展开 \[ \cos x=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+\cdots \] 随后截断一部分做近似计算 \[ \cos x\approx 1-\frac{x^2}{2!} \] 这样便会产生误差 \[ e=\frac{x^4}{4!}-\frac{x^6}{6!}+\cdots \] 通过交错级数审敛法可以得知这里的截断误差 \(e<\frac{x^4}{4!}\)。
舍入误差主要是针对计算机的存储问题，比如说标准通用 CPU 都支持的 IEEE 浮点数标准，可能会使得每次计算时都存在一定的舍入。
少量的舍入误差还好，但是一般数值方法都需要大量迭代，多次的舍入误差会导致结果出现很大的波动。
绝对误差与相对误差
假定 \(x^*\) 是准确值 \(x\) 的一个近似值。
那么其绝对误差（absolute error）为 \[ e=x-x^* \] 但是往往在实际计算中，我们无法确定准确值 \(x\) 的大小，也就是说我们无法确定绝对误差 \(e\)。
所以为了方便研究分析，实际会使用绝对误差限（absolute error bound）来断定其上限，符号为 \(\epsilon\)，定义为 \[ |e|=|x-x^*|\leq \epsilon \]
或者可以写作 \[ x^*-\epsilon\leq x\leq x^*+\epsilon \]
但是绝对误差不能满足所有的计算需要，因为近似值和真值可能十分小，那么得到的误差就会极其小，无法通过绝对误差的大小刻画近似值的准确程度，那么就需要相对误差（relative error） \[ \begin{align} e_r=\frac{e}{x}=\frac{x-x^*}{x}\\ e_r=\frac{e}{x^*}=\frac{x-x^*}{x^*} \end{align} \] 即绝对误差与真值或近似值的比值。
当 \(|e_r|\) 趋近于零时， \[ \frac{e}{x}-\frac{e}{x^*}=\frac{-e^2}{xx^*} \] 是 \(e_r\) 的高阶无穷小，可以忽略不计，所以无论哪个公式都是合理的，但一般默认是绝对误差与近似值的比值，因为这可计算。
同样的，也有相对误差限（relative error bound），符号为 \(\epsilon_r\)，定义为 \[ |e_r|=|\frac{e}{x^*}|\leq \epsilon_r \] 显然，误差限（error bound）不是唯一的，且 \(\frac{\epsilon}{|x^*|}\) 是 \(x^*\) 的一个相对误差限。
误差（error）函数写为 \[ e(x^*)=x-x^* \] 误差限（error bound）函数写为 \[ e(x^*)\leq\epsilon(x^*) \]
有效数字
从左到右，从第一个非零数字开始数，一共有几位就是几位有效数字
比如说 0.023 的有效数字是两位（23），而 0.0230 的有效数字是三位（230）。
对于有效数字是有舍入的（四舍五入），那么其误差限显然是小数点后最后一位的半个单位，一般写作 \(0.5\times 10^{-n}\)，这里的 \(n\) 不是上图的 \(n\) 个数，而是小数点后的位数。
比如说 0.023 的误差限是 0.0005，0.0230 的误差限是 0.00005。
有效数字（对于近似值而言）是有严格定义的，一定是满足四舍五入的舍入规则才能被作为有效数字，以 \(x=\sqrt{3}=1.7320508\cdots\) 为例，如果说有三个近似值 \(x_1=1.73\)，\(x_2=1.7321\)，\(x_3=1.7320\)，可以计算这三个近似值的误差限分别为 \[ \begin{align} |x-x_1|=|\sqrt{3}-1.73|=0.20508\cdots\times10^{-2}\leq0.5\times10^{-2}\\ |x-x_2|=|\sqrt{3}-1.7321|=0.508\cdots\times10^{-4}\leq0.5\times10^{-4}\\ |x-x_3|=|\sqrt{3}-1.7320|=0.508\cdots\times10^{-4}\leq0.5\times10^{-3}\\ \end{align} \] 故对于三个近似值而言，分别是
- \(x_1\) 精确到小数点后 2 位，有 3 位有效数字；
- \(x_2\) 精确到小数点后 4 位，有 5 位有效数字；
- \(x_3\) 精确到小数点后 3 位，有 4 位有效数字。
误差的传播
自变量的误差势必会传播到因变量，比如说 \(y=ax+b\)，假定自变量存在误差 \(e=x-x^*\)，那么有 \[ y-y^*=a(x-x^*) \] 即 \(e(y^*)=ae(x^*)\)，对于一个线性函数，误差会直接扩大 \(a\) 倍。
扩展到对一个连续可导函数 \(f(x)\)，有 $$ \[\begin{align} \Delta y=y-y^*=f(x)-f(x^*)\\ \end{align}\] \[ 当 $\Delta y$ 很小时，可以做一个线性近似得到 \] yf'(x)(x-x^)=f'(x)x \[ 也就是 \] e(y)f'(x^)e(x) \[ 那么其误差限也应该有 \] (y)|f'(x^)|(x) \[ 而对于相对误差，有 \] e_r(y^)=y|{y=y} \[ 当 $f$ 是多元函数时，绝对误差为 \] e(y^)=n{k=1}e(x^*) \[ 相对误差为 \] \[\begin{align} e_r(y^*)&=\frac{e(y^*)}{y^*}\\ &=\sum^n_{k=1}\frac{\partial f(x_1^*,\cdots,x_n^*)}{\partial x_k}\frac{e(x^*)}{y^*}\\ &=\sum^n_{k=1}\frac{\partial f(x_1^*,\cdots,x_n^*)}{\partial x_k}\frac{x^*}{y^*}e_r(x^*) \end{align}\] $$
四则运算的误差传播
加减法的误差传播公式 \[ \begin{cases} e(x_1\pm x_2)=e(x_1)\pm e(x_2)\\ e_r(x_1\pm x_2)=\frac{x_1}{x_1\pm x_2}e_r(x_1)\pm \frac{x_2}{x_1\pm x_2}e_r(x_2)\\ \end{cases} \] 如果扩展到 \(n\) 元，显然满足叠加定理，且相对误差按权重比例分配误差影响。
乘除法的误差传播公式 \[ \begin{cases} e(x_1x_2)\approx x_2e(x_1)\pm x_1e(x_2)\\ e_r(x_1x_2)\approx e_r(x_1)+e_r(x_2)\\ \end{cases} \]
\[ \begin{cases} e(\frac{x_1}{x_2})\approx \frac{1}{x_2}e(x_1)-\frac{x_1}{x_2^2}e(x_2)\\ e_r(x_1x_2)\approx e_r(x_1)-e_r(x_2)\\ \end{cases} \]
之所以乘除法的误差传播公式是约等于，是因为误差的传播公式是做线性近似，而这里的函数是二次的，如果对于高次函数（或者说多元乘法）误差限会很大。
误差传播例题
1
已知某物体行程 \(s\) 的近似值 \(s^*=800\mathrm{m}\)，所需时间 \(t\) 的近似值 \(t^*=35\mathrm{s}\)。如果已知 \(|s-s^*|\leq 0.5\)，\(|t-t^*|\leq 0.05\)，试求平均速度 \(v\) 的绝对误差限和相对误差限。
提取题目信息，\(\epsilon(s^*)=0.5,\epsilon(t^*)=0.05,\epsilon_r(s^*)=0.5/800,\epsilon_r(t^*)=0.05/35\)
首先计算平均速度 \(v^*\) 的绝对误差 \[ \begin{align} v^*&=\frac{s^*}{t^*}\\ e(v^*)&=e(\frac{s^*}{t^*})\\ &\approx \frac{1}{t^*}e(s^*)-\frac{s^*}{t^{*2}}e(t^*) \end{align} \] 根据绝对值不等式，有 \[ \begin{align} |e(v^*)|&\leq \frac{1}{t^*}e(s^*)+\frac{s^*}{t^{*2}}e(t^*)\\ &\leq \frac{1}{t^*}\epsilon(s^*)+\frac{s^*}{t^{*2}}\epsilon(t^*)\approx 0.0469\\ \end{align} \] 即 \(\epsilon(v^*)=0.0469\)。
那么 \[ \begin{align} \epsilon_r(v^*)=\frac{\epsilon(v^*)}{v^*}=0.00205 \end{align} \]
或者也可以套用公式 \[ \begin{align} e_r(v^*)&= e_r(\frac{s^*}{t^*})\\ &\approx e_r(s^*)-e_r(t^*) \end{align} \] 那么根据绝对值不等式，有 \[ \begin{align} |e_r(v^*)|&\leq |e_r(s^*)|+|e_r(t^*)|=0.00205 \end{align} \] 即 \(\epsilon_r(v^*)=0.00205\)。
2
已知 \(a=1.2031,b=0.978\) 是四舍五入的近似值，求 \(a+b\) 和 \(a\times b\) 各有几位有效数字？
根据有效数字的定义，得知 \[ \epsilon(a)=0.5\times 10^{-4}\\ \epsilon(b)=0.5\times 10^{-3}\\ \] 那么有 \[ \begin{align} |e(a+b)|=|e(a)+e(b)|\leq \epsilon(a)+\epsilon(b)=0.55\times 10^{-3}\leq 0.5\times 10^{-2}\\ |e(a\times b)|=|be(a)+ae(b)|\leq b\epsilon(a)+a\epsilon(b)=0.65045\times 10^{-3}\leq 0.5\times 10^{-2}\\ \end{align} \] 故全都精确到小数点后 2 位，大约计算 \(1.20+0.97>2,1.20\times 0.97>1\)，故结果均为三位有效数字。
3
设计算球体积要使相对误差限为 1%，问度量半径为 \(R\) 时允许的相对误差限是多少？
首先明确球体积公式为 \[ V(R)=\frac{4}{3}\pi R^3 \] 那么绝对误差为 \[ e[V(R)]=\frac{4}{3}\pi(3R^2)e(R)=4\pi R^2e(R) \] 那么相对误差为 \[ e_r[V(R)]=\frac{e[V(R)]}{V(R)}=3\frac{e(R)}{R}=3e_r(R) \] 那么 \[ e_r(R)=\frac{1}{3}e_r[V(R)]\leq \frac{1}{3}\epsilon[V(R)]=0.3333\% \]
误差改善
为了减少误差的影响，设计算法时需要考虑以下几个方面
- 选择数值稳定的算法
- 避免两个相近的数相减
- 避免大数吃小数现象
- 避免绝对值太小的数做除数
- 简化计算步骤，减少运算次数
选择数值稳定的算法
一个算法，如果在一定的条件下，其舍入误差在整个运算过程中能够得到控制或者舍入误差的增长不影响产生可靠的结果，则称该算法是数值稳定的，否则称为数值不稳定的。
比如说计算积分式的值 \[ I_n=\int^1_0 x^ne^{x-1}\mathrm{d}x \] 利用分部积分法可以求得递推公式 \[ I_n=1-nI_{n-1} \] 其误差为 \[ e(I_n^*)=-ne(I_{n-1}^*) \] 可以尝试数值实验：
不难发现在 \(n=8,9\) 时，误差就变得难以接受了。
显然分析可得，误差传播到 \(n\) 时， \[ |e(I_n^*)|=n!|e(I_0)^*| \] 以阶乘倍扩大误差，算法不稳定。
但是反之，从高积分值求得低积分值，即递推公式为 \[
I_{n-1}=\frac{1-I_n}{n}
\] 可知误差为 \[
e(I_{n-1}^*)=-\frac{e(I_n^*)}{n}
\] 
那么误差传播到 \(0\) 时， \[ |e(I_0^*)|=\frac{|e(I_n)^*|}{n!} \] 误差反而会越来越小。
对于同一数学问题，算法设计至关重要。只有选用数值稳定性好的算法，舍入误差才能够得到控制，求得比较精确的结果。
与数值稳定性相关的概念是良态和病态问题。
对于一个数值问题，如果输入数据的细微变化(或误差)，引起输出数据(即所求解)的很大变化，则称之为病态问题；如果输入数据的细微变化，不会引起解的很大变化，则称之为良态问题。
避免两个相近的数相减
两数差的相对误差关系式为 \[ e_r(x-y)=\frac{e(x)-e(y)}{x-y} \] 如果 \(x-y\) 趋近于 0，那么相对误差将会很大，有效数字位数将严重丢失。
一般来说会用三种方式防止这种现象产生：
因式分解 \[ \begin{align} \frac{1}{x}-\frac{1}{1+x}=\frac{1}{x(x+1)} \end{align} \]
分子分母有理化 \[ \begin{align} \sqrt{x+1}-\sqrt{x}=\frac{1}{\sqrt{x+1}+\sqrt{x}}\\ \end{align} \]
三角恒等式 \[ \begin{align} 1-\cos{x}=2\sin^2\frac{x}{2}\\ \frac{1-\cos{x}}{\sin x}=\frac{\sin x}{1-\cos{x}}\\ \end{align} \]
其他 \[ \ln{x_1}-\ln{x_2}=\ln{\frac{x_1}{x_2}} \]
例如说求二次方程 \(x^2-16x+1=0\) 的根，取 \(\sqrt{63}\approx 7.94\)，那么根据求根公式有 \[ x_1=8+\sqrt{63} \] 而第二个根为 \[ x_2=8-\sqrt{63}\approx 0.06 \] 但是如果避免减法 \[ x_2=\frac{1}{8+\sqrt{63}}\approx\frac{1}{15.94}\approx 0.0627 \] 反而精度更高。
避免大数“吃”小数现象
比如说计算机实例，使用 JavaScript
1 | 0.1+0.00000000000000001 |
结果令人吃惊
避免绝对值太小的数作除数
\[ e(\frac{x}{y})\approx \frac{1}{y}e(x)-\frac{x}{y^2}e(y)\\ \]
当 \(x\ll y\) 或 \(y\to 0\) 时，\(e(\frac{x}{y})\) 可能很大。
简化计算步骤，减少运算次数
比如说快速幂算法，本身求解 \(x^{31}\) 需要 30 次乘法运算，但是使用快速幂算法仅需要 8 次乘法运算 \[ x^{31}=x^1x^2x^4x^8x^{16} \]
插值法
插值问题
插值法是数值计算的基本课题，目标是使用简单函数为离散数组建立连续模型（连续关系）。
比如说在物理实验中测得一组被测对象的关系值，比如说压力传感器的力的大小和输出电压的大小，这是一组离散点，可能数据如下
| Force (N) | voltage (mV) |
|---|---|
| 20 | 810 |
| 21 | 892 |
| 24 | 1162 |
| 25 | 1260 |
| 27 | 1468 |
| 34 | 2322 |
| 50 | 5010 |
| 52 | 5418 |
| 55 | 6060 |
| 56 | 6282 |
这时使用简单的线性关系是无法完美描述该传感器的传感特性的，使用最小二乘法做拟合得到的结果如下
其相关系数 \(R^2=0.9926\)。
如果说这个时候应用多项式插值法，得到的拟合结果是
其相关系数 \(R^2=1\)，我们很好的拟合了该传感器的传感特性，这使得对该传感器的研究更加精确了。
上述只是为了说明插值法的意义，实际上最小二乘法是直线拟合法，不要求所有离散点必须通过该拟合直线；而插值法要求所有离散点必须通过拟合直线。
一般来说，我们将插值问题描述为，
已知 \(n+1\) 个节点 \((x_i,y_i)\)，其中 \(x_i\) 互不相同，不妨设 \(a=x_0<x_1<\cdots<x_n=b\)，那么求任一插值点 \(x^*\neq x_i\) 处的插值 \(y^*\)。
随后构造平面曲线 \(y=\varphi(x)\) 使其通过所有节点，即 \(y_i=\varphi(x_i)\)。
根据上述定义，可以得到插值法的定义为：
设函数 \(y=f(x)\) 在区间 \([a,b]\) 上有定义，且已知在点 \(a=x_0<x_1<\cdots<x_n=b\) 上的值 \(f(x_i)=y_i(i=0,1,\cdots,n)\)，若存在一简单函数 \(\varphi(x)\)，使 \[ \varphi(x_i)=y_i(i=0,1,\cdots,n) \] 成立，就称 \(\varphi(x)\) 为 \(f(x)\) 的插值函数，点 \(x_i(i=0,1,\cdots,n)\) 为插值节点，包括插值节点的区间 \([a,b]\) 称为插值区间，求插值函数 \(\varphi(x)\) 的方法称为插值法。
对于被插函数 \(f(x)\) 和插值函数 \(\varphi(x)\) 在节点 \(x_i\) 处的函数值必然相等，但在节点外 \(\varphi(x)\) 的值可能就会偏离 \(f(x)\)，因此 \(\varphi(x)\) 近似代替 \(f(x)\) 必然存在截断误差。
如果说 \(\varphi_n(x)\) 为次数不超过 \(n\) 的代数多项式 \[ \varphi_n(x)=a_0+a_1x+\cdots+a_nx^n \] 那么就称 \(\varphi_n(x)\) 为插值多项式，相应的插值法称为多项式插值法；
若 \(\varphi_n(x)\) 为分段多项式，就称为分段插值法。
一般来说会选择多项式函数作为插值法的目标插值函数，这是因为其计算简便、性质简单。
多项式插值法
唯一解
由定义可以构造方程组 \[ \begin{cases} a_0+a_1x_0+a_2x_0^2+\cdots+a_nx_0^n=y_0\\ a_0+a_1x_1+a_2x_1^2+\cdots+a_nx_1^n=y_1\\ \vdots\\ a_0+a_1x_n+a_2x_n^2+\cdots+a_nx_n^n=y_n\\ \end{cases} \] 那么该方程组的行列式为 \[ \det{(\mathcal{A})}= \begin{vmatrix} 1&x_0&\cdots&x_0^n\\ 1&x_1&\cdots&x_1^n\\ \vdots&\vdots&\ddots&\vdots\\ 1&x_n&\cdots&x_n^n\\ \end{vmatrix}=\prod(x_i-x_j) \] 只要满足任意 \(x_i\neq x_j\) 且插值多项式 \(\varphi(x)\) 的阶数 \(\leq n\)，那么该方程组就有唯一解。
误差估计
上文提到插值法会造成截断误差，而插值多项式的截断误差可以写为 \[ R_n(x)=f(x)-\varphi_n(x) \] 也称为插值多项式的余项。
截断误差的证明过程是复杂的，这里仅给出结论 \[ R_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}w_{n+1}(x),w_{n+1}(x)=\prod^n_{i=0}(x-x_i) \] 但是这个表达式只有在 \(f(x)\) 的高阶导数存在时才能应用，且 \(\xi\) 无法确定，但可以估算 \(\xi\in (a,b)\)，如果可以求出 \[ \max_{a<x<b}|f^{(n+1)}(x)|=M_{n+1} \] 那么可以得到误差限 \[ |R_n(x)|\leq \frac{M_{n+1}}{(n+1)!}|w_{n+1}(x)| \]
拉格朗日插值法
线性插值
先讨论 \(n=1\) 的情况，\(\varphi_1(x)=a_0+a_1x\)，这也被称为线性插值
那么根据两点，可以构造出方程组 \[ \begin{cases} \varphi_1(x_0)=a_0+a_1x_0=y_0=f(x_0)\\ \varphi_1(x_1)=a_0+a_1x_1=y_1=f(x_1)\\ \end{cases} \] 使用两点式公式，可以得到 \[ \varphi_1(x)=\frac{x-x_1}{x_0-x_1}y_0+\frac{x-x_0}{x_1-x_0}y_1 \] 似乎可以抽象出一个函数 \(l_i(x)\) 将插值多项式表示为 \[ \varphi_1(x)=l_0(x)y_0+l_1(x)y_1 \]
抛物插值
再讨论 \(n=2\) 的情况，\(\varphi_2(x)=a_0+a_1x+a_2x^2\)，这也被称为抛物插值
同样的，可以构造出 \[ \begin{align} \varphi_2(x)&=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}y_0\\ &+\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}y_0\\ &+\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}y_0 \end{align} \] 似乎同样可以抽象出一个函数 \(l_i(x)\) 将插值多项式表示为 \[ \varphi_2(x)=l_0(x)y_0+l_1(x)y_1+l_2(x)y_2 \]
一般形式
更一般地，可以发现插值多项式就是 \(n+1\) 个线性函数的线性组合。
将规律延伸到 \(n\) 次，那么可以得到 \[ L_n(x)=l_0(x)y_0+\cdots+l_n(x)y_n=\sum^n_{i=0}y_il_i(x) \] 其中，\(L_n\) 称为拉格朗日插值函数，\(l_i(x)\) 称为拉格朗日基函数，定义为 \[ l_i(x_j)=\begin{cases} 0,&j\neq i\\ 1,&j=i \end{cases} \] 为了计算方便，可以写出数值形式为 \[ l_i(x)=\frac{(x-x_0)\cdots(x-x_{n})}{(x_i-x_0)\cdots(x_i-x_n)}=\prod^n_{j=0,j\neq i}\frac{x-x_j}{x_i-x_j} \] 而插值多项式为 \[ L_n(x)=\sum^n_{i=0}y_il_i(x)=\sum^n_{i=0}\left(\prod^n_{j=0,j\neq i}\frac{x-x_j}{x_i-x_j}\right) \]
除此之外还有 \(w_{n+1}(x)\) 形式下的拉格朗日插值函数 \[ \begin{align} w_{n+1}(x)&=\prod^n_{j=0}(x-x_j)\\ w_{n+1}'(x_i)&=\prod^n_{j=0,j\neq i}(x_i-x_j)\\ L_n(x)&=\sum^n_{i=0}\frac{w_{n+1}(x)}{(x-x_i)w_{n+1}'(x_i)}y_i \end{align} \]
优缺点
拉格朗日插值法的公式简洁，结构紧凑，理论分析方便，但改变一个数据点时，全部的插值基函数都将发生改变，基函数没有承袭性，且数据点很多时，计算量很大，为 \(2n^2+2n\)。
牛顿插值法
规律延伸
同理首先考虑线性插值 \[ \varphi_1(x)=y_0+\frac{y_1-y_0}{x_1-x_0}(x-x_0) \] 但此时我们使用点斜式，再来观察抛物插值 \[ \varphi_2(x)=y_0+\frac{y_1-y_0}{x_1-x_0}(x-x_0)+\frac{\frac{y_2-y_0}{x_2-x_0}-\frac{y_1-y_0}{x_1-x_0}}{x_2-x_1}(x-x_0)(x-x_1) \] 我们可以从中抽象出一个函数 \(p_i\)，使得 \[ \varphi_n(x)=p_0+p_1(x-x_0)+p_2(x-x_0)(x-x_1)+\cdots+p_n\prod^{n-1}_{i=0}(x-x_i) \] 如果说接着推导 \(p_n\)，那么其形式会十分复杂，故需要引入差商的概念来简化运算。
差商
对函数 \(f(x)\) 有 \(x_0,x_1,\cdots\) 一系列互不相同的点，称 \[ f[x_i,x_j]=\frac{f(x_i)-f(x_j)}{x_i-x_j}(i\neq j) \] 为 \(f(x)\) 关于点 \(x_i, x_j\) 的一阶差商（均差）。
类比斜率。
中括号是为了有意区分函数值和差商的概念，可以认为 \(f[x_i]=f(x_i)\) 是零阶差商。
而称 \[ f[x_i,x_j,x_k]=\frac{f[x_i,x_j]-f[x_j,x_k]}{x_i-x_k}(i\neq j\neq k) \] 为 \(f(x)\) 关于点 \(x_i,x_j,x_k\) 的二阶差商（一阶差商的差商）。
那么 \(k\) 阶差商如下 \[ f[x_0,x_1,\cdots,x_k]=\frac{f[x_0,x_1,\cdots,x_{k-1}]-f[x_1,x_2,\cdots,x_{k}]}{x_0-x_k} \]
差商的性质
差商具有线性性，若 \(f(x)=a\varphi(x)+b\phi(x)\)，则对于任意正整数 \(k\) 都有 \[ f[x_0,x_1,\cdots,x_k]=a\varphi[x_0,x_1,\cdots,x_k]+b\phi[x_0,x_1,\cdots,x_k] \] \(k\) 阶差商可以由函数值的线性组合表示，即 \[ f[x_0,x_1,\cdots,x_k]=\sum^k_{i=0}\frac{f(x_i)}{(x_i-x_0)\cdots(x_i-x_{k})}=\sum^n_{i=0}\frac{f(x_i)}{w_{k+1}'(x_i)} \] 差商具有对称性，即任意调换节点的次序，差商的值不变，如 \[ f[x_0,x_1,x_2]=f[x_0,x_2,x_1]=f[x_2,x_1,x_0] \] 若 \(f(x)\) 是 \(x\) 的 \(n\) 次多项式，则其一阶差商 \(f[x,x_i]\) 是 \(n-1\) 次多项式。
若 \(f(x)\) 在 \([a,b]\) 上存在 \(n\) 阶导数，且 \(x_i\in[a,b]\)，那么 \(n\) 阶差商与导数的关系为 \[ f[x_0,x_1,\cdots,x_n]=\frac{f^{(n)}(\xi)}{n!},\xi\in(a,b) \]
牛顿多项式
用差商的形式表示为 \[ \begin{align} N_n(x)&=f(x_0)\\ &+(x-x0)f[x_0,x_1]\\ &+(x-x_0)(x-x_1)f[x_0,x_1,x_2]\\ &+\cdots\\ &+(x-x_0)\cdots(x-x_{n-1})f[x_0,\cdots,x_n] \end{align} \] 其中牛顿多项式的误差也可以用差商表示 \[ R_n(x)=(x-x_0)\cdots(x-x_n)f[x,x_0,\cdots,x_n]=w_{n+1}(x)f[x,x_0,\cdots,x_n] \] 牛顿插值的优点是，每增加一个节点，插值多项式只增加一项，因此便于递推运算和程序设计，即 \[ N_n(x)=N_{n-1}(x)+f[x_0,\cdots,x_n](x-x_0)\cdots(x-x_{n-1}) \] 且与拉格朗日插值法相比，多项式等价，误差一致。
一般来说，会根据 \(x_i,y_i\) 的值构造差分表，比如说以下数据
| \(x_i\) | \(y_i\) |
|---|---|
| 0.40 | 0.41075 |
| 0.55 | 0.57815 |
| 0.65 | 0.69675 |
| 0.80 | 0.88811 |
| 0.90 | 1.02652 |
| 1.05 | 1.25382 |
我们可以构造这样的差分表
| x_i | y_i | 一阶差商 | 二阶差商 | 三阶差商 | 四阶差商 | 五阶差商 |
|---|---|---|---|---|---|---|
| 0.40 | 0.41075 | |||||
| 0.55 | 0.57815 | 1.11600 | ||||
| 0.65 | 0.69675 | 1.18600 | 0.28000 | |||
| 0.80 | 0.88811 | 1.27573 | 0.35893 | 0.197333 | ||
| 0.90 | 1.02652 | 1.38410 | 0.43348 | 0.212952 | 0.0312381 | |
| 1.05 | 1.25382 | 1.51533 | 0.52493 | 0.228667 | 0.0314286 | 0.000380952 |
其中，对角线即为所求的差商 \(f[x_0],f[x_0,x_1],\cdots,f[x_0,\cdots,x_n]\)。
差分
定义
向前差分，Forward Difference（F-差分） \[ \Delta f_k=f_{k+1}-f_k \] 向后差分，Backward Difference（B-差分） \[ \nabla f_k=f_k-f_{k-1} \]
性质
对于向前差分而言，具有连续性 \[
\begin{align}
\Delta f_k&=f_{k+1}-f_k\\
\Delta^2 f_k&=\Delta f_{k+1}-\Delta f_k\\
&\vdots\\
\Delta^m f_k&=\Delta^{m-1} f_{k+1}-\Delta^{m-1} f_k
\end{align}
\] 
线性性 \[ f(x)=a\varphi(x)+b\phi(x)\to \Delta^m f_k=a\Delta^m \varphi_k+b\Delta^m \phi_k \] 线性表示 \[ \Delta^m f_k=\sum^m_{j=0}(-1)^j\binom{j}{m} f_{k+m-j} \] 差分与差商的关系 \[ f[x_k,\cdots,x_{k+m}]=\frac{1}{m!}\frac{1}{h^m}\Delta^m f_k \] 差分与导数的关系 \[ f^{(m)}(\xi)=\frac{\Delta^m f_k}{h^m} \]
等距节点牛顿插值法
实际使用中，通常我们得到的是等距节点下的数据值，例如说 \((0,1.002),(1,1.208),\cdots,(x_n,y_n)\)，满足 \(h=x_i-x_{i-1}\) 为常数（等差）。
不妨将差分与差商的关系代入牛顿插值公式，可以得到 \[ N_n(x)=f_0+\frac{\delta f_0}{h}w_1(x)+\frac{\Delta^2 f_0}{2!h^2}w_2(x)+\cdots+\frac{\Delta^n f_0}{n!h^n} \] 其中，\(h\) 为等距节点的等差。
为方便起见，将 \(w_{n+1}(x)\) 化为无关 \(x\) 的函数，即不妨取 \(x=x_0+th\)，其中 \(t\) 为间距倍数，使得 \[ w_{n+1}(x)=w_{n+1}(x_0+th)=(t-0)(t-1)\cdots(t-n)h^{n+1} \] 那么牛顿多项式化为 F-插值公式 \[ N_n(x)=f_0+t\Delta f_0+\frac{t(t-1)}{2!}\Delta^2 f_0+\cdots+\frac{t(t-1)\cdots(t-n+1)}{n!}\Delta^n f_0 \] 其余项（误差）化为 \[ R_n(x)=\frac{t(t-1)\cdots(t-n)}{(n+1)!}h^{n+1}f^{(n+1)}(\xi) \] 误差限为 \[ |R_n(x)|\leq \left|\frac{t(t-1)\cdots(t-n)}{(n+1)!}h^{n+1}M_{n+1}\right|,|f^{(n+1)}(\xi)|\leq M_{n+1} \] 同理，B-差分也有类似的关系，可以得到 B-插值公式 \[ N_n(x)=f_n+t\nabla f_n+\frac{t(t-1)}{2!}\nabla^2 f_n+\cdots+\frac{t(t-1)\cdots(t-n+1)}{n!}\nabla^n f_n \] 当插值点 \(x\) 接近数据表头时，一般采用向前插值法；
当插值点 \(x\) 接近数据表尾时，一般采用向后插值法。
看的是 \(x\) 接近 \(f_0\) 还是 \(f_n\)，越相近截断误差越小。
埃尔米特插值法
埃尔米特插值多项式
许多实际的插值问题不但要求节点处函数值相同，还要求它的函数有相同的一阶、二阶导数，满足这种要求的插值多项式就是埃尔米特插值多项式。
故需要 \(2n+2\) 个条件（数据点），其满足 \[ \begin{cases} H_{2n+1}(x_i)=f(x_i)\\ H'_{2n+1}(x_i)=f'(x_i) \end{cases} \] 余项为 \[ R(x)=\frac{f^{2n+2}(\xi)}{(2n+2)!}w_{n+1}^2(x) \] 为了满足条件，那么埃尔米特插值多项式的一般公式为 \[ H(x)=\sum^n_{i=0}[\alpha_i(x)y_i+\beta_i(x)y_i'] \] 其中，\(\alpha_i(x)\) 主管函数值，其导数值要为 0；\(\beta_i(x)\) 主管导数值，函数值为 0。
一个满足条件的函数为 \[ \begin{cases} \alpha_i(x)=[1-2(x-x_i)l_i'(x_i)][l_i(x)]^2\\ \beta_i(x)=(x-x_i)[l_i(x)]^2 \end{cases} \]
龙格现象
龙格（Runge）现象指的是当对函数进行高次插值时，插值多项式 \(\varphi_n(x)\) 对 \(f(x)\) 的逼近效果反而会更差。
比如说 \(f(x)=1/(1+x_i^2)\)，对该函数进行 \(n\) 次插值，可以得到不同次数下拉格朗日插值多项式的比较图
并不是插值多项式的次数越高，插值效果越好，精度也不一定是随次数的提高而升高。
故不适宜在大范围使用高次代数插值。
分段线性插值
基于龙格现象的理论，既然大范围使用代数插值效果不好，那么就说明小范围使用代数插值效果是良好的，那对于多个紧凑的数据点，不妨通过分段插值的方式，使得各个插值点用折线段连接起来，便得到了一个分段插值函数。
在一定区间内，可证明的是，计算量与 \(n\) 无关，\(n\) 越大误差越小。
分段三次埃尔米特插值
三次样条插值
函数逼近与曲线拟合
为什么需要曲线拟合
插值法的核心要求是，插值函数要经过数据点中的所有点，但实际上的实验数据并不都是完美接近理论的，总会存在一定的噪声使得数据点偏离，这时插值法实际上并不能很好地拟合原函数。
此时我们可以回归原本，我们为什么想要拟合函数经过数据点中的所有点呢？
因为希望说拟合函数逼近原函数。
这也是函数逼近问题：在实际应用中常需为解析式比较复杂的函数寻找一个简单函数来近似代替它，并要求其误差在某种度量意义下最小。
曲线拟合问题：在实际应用中，往往并不需要曲线通过给定的数据点，而只要求用曲线（函数）近似代替给定的列表函数时，其误差在某种度量意义下最小。
最佳逼近问题：要求在被插函数的定义区间上都有较好的近似。
函数逼近
我们定义，近似函数求得的近似值 \(y_i^*\) 与观测值 \(y_i\) 的差 \(\delta=y_i-y_i^*\) 称为残差。
那么为了度量近似函数的好坏，常用的准则为
- 使残差的最大绝对值 \(\max_i|\delta_i|\) 最小
- 使残差的绝对值之和 \(\sum|\delta_i|\) 最小
- 使残差的平方和 \(\sum\delta_i^2\) 最小
准则 1、2 有绝对值，计算时不方便。
准则 1 求近似函数的方法称为函数的最佳一致逼近。
准则 3 确定参数求近似函数的方法称为函数的最佳平方逼近（对连续函数），也称曲线拟合的最小二乘法（对离散数据）。
内积与内积空间
在曾经学过，两个向量的内积定义和范数定义是 \[ \begin{cases} (\mathbf{x},\mathbf{y})=\sum x_iy_i\\ \|\mathbf{x}\|_2=\sqrt{(\mathbf{x},\mathbf{x})}=\sqrt{\sum x_i^2} \end{cases} \] 在某些问题中，向量中每个分量的重要性是不相同的，比如说在机器学习中，可能前几个分量的重要性高，而后几个分量的重要性低，此时我们就会引入权重 \(w_i\) 来表征各个分量的重要性，即 \[ \begin{cases} (\mathbf{x},\mathbf{y})=\sum w_ix_iy_i\\ \|\mathbf{x}\|_2=\sqrt{(\mathbf{x},\mathbf{x})}=\sqrt{\sum w_ix_i^2} \end{cases} \] 向量是离散的，如果说扩展内积的定义，延伸到连续函数上，那么求和就变成了积分，假定两个连续函数 \(f,g\in C[a,b]\)，那么函数内积的定义为 \[ \begin{cases} (f,g)=\int^b_af(x)g(x)\mathrm{d}x\\ \|f\|_2=\sqrt{(f,f)}=\sqrt{\int^b_af^2(x)\mathrm{d}x} \end{cases} \] 带权重的函数内积同理。
在机器学习中，权函数 \(w(x)=1\) 意味着区间中所有点具有同等重要性；
权函数 \(w(x)=\frac{1}{\sqrt{1-x^2}}\) 意味着靠近区间端点处的点更加重要而区间中点处的点最不重要。
正交多项式
离散型正交多项式可作为基函数用于曲线拟合，提高算法的稳定性；
连续型正交多项式可作为基函数用于高斯型求积公式的构造。
正交的定义是内积为 0，对于函数内积也是一样的，即如果满足以下式子 \[ (f,g)=\int^b_af(x)g(x)\mathrm{d}x=0 \] 那么我们称 \(f(x)\) 与 \(g(x)\) 在 \([a,b]\) 上正交。
如果说函数族 \(\varphi_0(x),\varphi_1(x),\cdots,\varphi_n(x)\) 满足关系 \[ \begin{align} (\varphi_i,\varphi_j)=\int_a^b\varphi_i(x)\varphi_j(x)\mathrm{d}x=\begin{cases} 0,&i\neq j\\ A_j,&i=j \end{cases} \end{align} \] 那么称 \(\{\varphi_i(x)\}\) 是 \([a,b]\) 上的正交函数系。
特别地，如果 \(A_j=1\)，那么称之为标准正交函数系。
带权也是同理。
正交三角函数系
三角函数系为 \[ \{1,\cos x,\sin x,\cos 2x,\sin 2x,\cdots\},x\in[-\pi,\pi] \] 可以计算得到任意两个不同函数在 \([-\pi,\pi]\) 上的内积（积分）为 0 \[ \begin{align} \int^\pi_{-\pi}\cos {(nx)} \mathrm{d}x=0\\ \int^\pi_{-\pi}\sin {(nx)} \mathrm{d}x=0\\ \int^\pi_{-\pi}\sin {(mx)}\cos {(nx)} \mathrm{d}x=0\\ \int^\pi_{-\pi}\cos {(mx)}\sin {(nx)} \mathrm{d}x=0\\ \end{align} \] 而任意两个相同函数在 \([-\pi,\pi]\) 上的内积（积分）为正数 \[ \begin{align} (1,1)=2\pi\\ (\sin{(kx)},\sin{(kx)})=\pi\\ (\cos{(kx)},\cos{(kx)})=\pi \end{align} \]
施密特正交化
在线性代数中学习可知，对于一组线性无关向量组 \(\alpha_1,\alpha_2,\cdots,\alpha_m\)，那么可以通过以下方式计算得到 \[ \begin{cases} \beta_1=\alpha_1\\ \beta_2=\alpha_2-\frac{(\alpha_2,\beta_1)}{(\beta_1,\beta_1)}\beta_1\\ \vdots\\ \beta_m=\alpha_m-\frac{(\alpha_m,\beta_1)}{(\beta_1,\beta_1)}\beta_1-\frac{(\alpha_m,\beta_2)}{(\beta_2,\beta_2)}\beta_2-\cdots-\frac{(\alpha_m,\beta_{m-1})}{(\beta_{m-1},\beta_{m-1})}\beta_{m-1} \end{cases} \] 其中，\(\beta_1,\beta_2,\cdots,\beta_m\) 就是一组正交向量组。
如果再令 \[ e_i=\frac{\beta_i}{\|\beta_i\|} \] 那么 \(e_1,e_2,\cdots,e_m\) 就是一组标准正交向量组。
那么同样的，如果说我们对一族线性无关的幂函数 \[ \{1,x,\cdots,x^m\} \] 通过施密特正交化的方式也可以得到一族正交函数系。
可以成立如下递推关系式 \[ \varphi_{n+1}(x)=(x-\alpha_n)\varphi_n(x)-\beta_n\varphi_{n-1}(x) \] 其中， \[ \begin{cases} \varphi_0(x)=1,\varphi_{-1}(x)=0\\ \alpha_n=\displaystyle\frac{(x\varphi_n(x),\varphi_n(x))}{(\varphi_n(x),\varphi_n(x))}\\ \beta_n=\displaystyle\frac{(\varphi_n(x),\varphi_n(x))}{(\varphi_{n-1}(x),\varphi_{n-1}(x))} \end{cases} \] 在区间 \([0,1]\) 上的结果如下 \[ \begin{cases} \varphi_0(x)=1\\ \varphi_1(x)=x-\frac{1}{2}\\ \varphi_2(x)=x^2-x+\frac{1}{6}\\ \varphi_3(x)=x^3-\frac{3}{2}x^2+\frac{3}{5}x-\frac{1}{20}\\ \varphi_4(x)=x^4-2x^3+\frac{9}{7}x^2-\frac{2}{7}x+\frac{1}{70}\\ \varphi_5(x)=x^5-\frac{5}{2}x^4+\frac{20}{9}x^3-\frac{5}{6}x^2+\frac{5}{42}x-\frac{1}{252}\\ \varphi_6(x)=x^6-3x^5+\frac{75}{22}x^4-\frac{20}{11}x^3+\frac{5}{11}x^2-\frac{1}{22}x+\frac{1}{924}\\ \end{cases} \] 在区间 \([-1,1]\) 的结果会更简单些（因为奇函数积分为 0），即 \[ \begin{cases} \varphi_0(x)=1\\ \varphi_1(x)=x\\ \varphi_2(x)=x^2-\frac{1}{3}\\ \varphi_3(x)=x^3-\frac{3}{5}x\\ \varphi_4(x)=x^4-\frac{6}{7}x^2+\frac{3}{35}\\ \varphi_5(x)=x^5-\frac{10}{9}x^3+\frac{5}{21}x\\ \varphi_6(x)=x^6-\frac{15}{11}x^4+\frac{5}{11}x^2-\frac{5}{231}\\ \varphi_7(x)=x^7-\frac{231}{143}x^5+\frac{105}{143}x^3-\frac{25}{429}x \end{cases} \] 常用的正交多项式有
- 勒让德(Legendre)多项式
- 切比雪夫(Chebyshev)多项式
- 拉盖尔(Laguerre)多项式
- 第二类切比雪夫(Chebyshev)多项式
- 埃尔米特(Hermite)多项式
最佳平方逼近
在函数逼近问题中，以残差的平方 \(\delta_i^2\) 最小作为准则就是最佳平方逼近方法，而对应的拟合函数称为最佳平方逼近函数。
设 \(f(x)\) 在区间 \([a,b]\) 上连续，\(\varphi_i(x)\) 是一组线性无关的连续函数，即 \[ H=\mathrm{span}\{1,x,\cdots,x^n\} \] 若 \(\varphi(x)=a_0+a_1x+\cdots+a_nx^n\)，使得 \[ F(a_0,a_1,\cdots,a_n)=\int^b_a[f(x)-\varphi(x)]^2\mathrm{d}x \] 最小，则称 \(\varphi(x)\) 为 \(f(x)\) 在空间 \(H\) 中的 \(n\) 次最佳平方逼近多项式。
根据多元函数取极值的必要条件可得 \[ \frac{\partial F}{\partial a_i}=-2\int^b_a[f(x)-\varphi(x)]\varphi_i(x)\mathrm{d}x=0 \] 即 \[ \sum^n_{j=0}a_j\int^b_a\varphi_i(x)\varphi_j(x)\mathrm{d} x=\int^b_a f(x)\varphi_i(x)\mathrm{d}x,i=0,\cdots,n \] 将积分转化为函数内积的表现形式，为 \[ \sum^n_{j=0}(\varphi_i(x),\varphi_j(x))a_j=(f(x),\varphi_i(x)),i=0,\cdots,n \] 改写为矩阵形式为 \[ \begin{bmatrix} (\varphi_0,\varphi_0)&(\varphi_0,\varphi_1)&\cdots&(\varphi_0,\varphi_n)\\ (\varphi_1,\varphi_0)&(\varphi_1,\varphi_1)&\cdots&(\varphi_1,\varphi_n)\\ \vdots&\vdots&\ddots&\vdots\\ (\varphi_n,\varphi_0)&(\varphi_n,\varphi_1)&\cdots&(\varphi_n,\varphi_n) \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\ \vdots\\ a_n \end{bmatrix}= \begin{bmatrix} (f,\varphi_0)\\ (f,\varphi_1)\\ \vdots\\ (f,\varphi_n) \end{bmatrix} \] 由于 \(\varphi_0,\varphi_1,\cdots,\varphi_n\) 线性无关，故方程组仅有唯一解。
最佳平方逼近示例
比如说 \(f(x)=\cos{(\pi x)}\) 在 \([0,1]\) 上的一次最佳平方逼近多项式。
取 \(\varphi_0(x)=1,\varphi_1(x)=x\)，即 \[
H=\mathrm{span}\{1,x\}
\] 那么 \[
\begin{cases}
(\varphi_0,\varphi_0)=\displaystyle\int^1_0 1\mathrm{d}x=1\\
(\varphi_0,\varphi_1)=(\varphi_1,\varphi_0)=\displaystyle\int^1_0
x\mathrm{d}x=\frac{1}{2}\\
(\varphi_1,\varphi_1)=\displaystyle\int^1_0 x^2\mathrm{d}x=\frac{1}{3}\\
(f,\varphi_0)=\displaystyle\int^1_0 \cos{(\pi x)}\mathrm{d}x=0\\
(f,\varphi_1)=\displaystyle\int^1_0 x\cos{(\pi
x)}\mathrm{d}x=-\frac{2}{\pi^2}\\
\end{cases}
\] 那么列出方程 \[
\begin{cases}
a_0+\frac{1}{2}a_1=0\\
\frac{1}{2}a_0+\frac{1}{3}a_1=-\displaystyle\frac{2}{\pi^2}
\end{cases}
\] 可以解得 \[
\begin{cases}
a_0=1.2159\\
a_1=-2.4317
\end{cases}
\] 那么一次最佳平方逼近多项式为 \[
\varphi(x)=1.2159-2.4317x
\] 也可以根据极值条件有 \[
F(a_0,a_1)=\int^1_0[\cos{(\pi x)-a_0-a_1x}]^2\mathrm{d}x
\] 使得 \[
\begin{cases}
\displaystyle\frac{\partial F}{\partial a_0}=2\int^1_0-[\cos{(\pi
x)-a_0-a_1x}]\mathrm{d}x=2a_0+a_1=0\\
\displaystyle\frac{\partial F}{\partial a_1}=2\int^1_0-x[\cos{(\pi
x)-a_0-a_1x}]\mathrm{d}x=2(\frac{2}{\pi^2}+\frac{1}{2}a_0+\frac{1}{3}a_1)=0\\
\end{cases}
\] 
系数矩阵
从上文不难发现，\((\varphi_i,\varphi_j)\) 都是常数，即 \[ \begin{bmatrix} 1&1/2&\cdots&1/(n+1)\\ 1/2&1/3&\cdots&1/(n+2)\\ \vdots&\vdots&\ddots&\vdots\\ 1/(n+1)&1/(n+2)&\cdots&1/(2n+1) \end{bmatrix} \] 而此时如果 \(H\) 取正交函数系，那么可以进一步简化为 \[ \begin{bmatrix} (\varphi_0,\varphi_0)&0&\cdots&0\\ 0&(\varphi_1,\varphi_1)&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&(\varphi_n,\varphi_n) \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\ \vdots\\ a_n \end{bmatrix}= \begin{bmatrix} (f,\varphi_0)\\ (f,\varphi_1)\\ \vdots\\ (f,\varphi_n) \end{bmatrix} \] 对于不同的正交多项式，其系数矩阵均为常数。
最小二乘法
最佳平方逼近是针对于函数逼近问题的，但实际应用中我们更希望的是给予一堆数据点，然后对其作曲线拟合，往往我们并不知道原函数。
由于数据点是离散的，套用最佳平方逼近的准则是行不通的，需要做一点点改进 \[ \sum \delta_i^2=\sum [y_i-P_m(x_i)]^2 \] 使得残差平方和最小就是最小二乘法的准则，其中 \(P_m(x)\) 称为这组数据的最小二乘 \(m\) 次拟合多项式。
与最佳平方逼近同样的道理，令 \[ F(a_0,a_1,\cdots,a_n)=\sum_{i=1}^n[y_i-P_m(x_i)]^2 \] 那么由极值的必要条件，得 \[ \frac{\partial F}{\partial a_i}=-2\sum_{i=1}^n [y_i-P_m(x_i)]\varphi_j(x_i)\mathrm{d}x=0 \] 整理可得 \[ \sum_{i=1}^n P_m(x_i)=\sum_{i=1}^n y_i \] 那么转化为正规方程组 \[ \begin{cases} \displaystyle a_0n+a_1\sum_{i=1}^n x_i+\cdots+a_m\sum_{i=1}^n x_i^m=\sum_{i=1}^n y_i\\ \displaystyle a_0\sum_{i=1}^n x_i+a_1\sum_{i=1}^n x_i^2+\cdots+a_m\sum_{i=1}^n x_i^{m+1}=\sum_{i=1}^n y_ix_i\\ \vdots\\ \displaystyle a_0\sum_{i=1}^n x_i^m+a_1\sum_{i=1}^n x_i^{m+1}+\cdots+a_m\sum_{i=1}^n x_i^{2m}=\sum_{i=1}^n y_ix_i^m\\ \end{cases} \] 转化为矩阵形式为 \[ \mathbf{Q}^T\mathbf{Q}\mathbf{A}=\mathbf{Q}^T\mathbf{Y} \] 其中， \[ \begin{align} \mathbf{Q}&=\begin{bmatrix} 1&x_1&\cdots&x_1^m\\ 1&x_2&\cdots&x_2^m\\ \vdots&\vdots&\ddots&\vdots\\ 1&x_n&\cdots&x_n^m\\ \end{bmatrix}\\ \mathbf{Q}^T\mathbf{Q}&=\begin{bmatrix} n&\sum_1^n x_i&\cdots&\sum_1^n x_i^m\\ \sum_1^n x_i&\sum_1^n x_i^2&\cdots&\sum_1^n x_i^{m+1}\\ \vdots&\vdots&\ddots&\vdots\\ \sum_1^n x_i^m&\sum_1^n x_i^{m+1}&\cdots&\sum_1^n x_i^{2m}\\ \end{bmatrix}\\ \mathbf{A}&=\begin{bmatrix} a_0\\ a_1\\ \vdots\\ a_m \end{bmatrix},\mathbf{Y}=\begin{bmatrix} y_0\\ y_1\\ \vdots\\ y_n \end{bmatrix} \end{align} \] 该方程组的解对应的多项式就是满足已知数据组的最小二乘 \(m\) 次拟合多项式由于多项式拟合中的正规方程组为线性方程组，所以称其为线性最小二乘问题。
非线性最小二乘拟合
高中就学过，如何对非线性数据作最小二乘拟合。
