曲线拟合和最佳逼近
Published:
曲线拟合本质上为数值逼近问题，本文将从曲线拟合的算法出发，并且分析最佳平方逼近和最佳一致逼近。
函数插值问题，是要求函数在某些点处的取值和给定曲线的函数值相等。
关于插值，可以参见，函数的插值
但插值出来的函数整体效果并不一定很好，因此可以考虑进行拟合，本质上是在给定的函数空间内寻找到逼近性质最好的一个函数。
Polynomial Fitting
多项式拟合是常见的拟合手段，对于待拟合的函数$y$,给定一组多项式基$X$，由于多项式可被这组基用系数$c$线性组合，因此等价于寻找拟合效果最好的函数$Xc$.
OLS（Ordinary Least Square）
普通最小二乘问题（OLS）等价于将噪声看作在$y$之上，而$X$是固定的。
\[\min_{y + \Delta y = Xc} \Vert \Delta y \Vert_F\]该最小二乘问题可以写成正规方程的形式，并且解得，$c = (X^T X)^{-1} X y$
通常可以使用QR分解等方法求解，$X = QR,Q^TQ=I$, 则正规方程等价于，
\[\begin{align} X^T X c&= X^T y\\ R^T R c &= R^T Q^T y \\ R c &= Q^T y \end{align}\]因此，只要解方程组$Rc = Q^T y$即可，而由于$R$为上三角阵，其求解是方便的，依次代入即可
TLS（Total Least Square）
总体最小二乘问题（TLS）将噪声看作同时作用在$X,y$之上，此时问题变为，
\[\min_{y +\Delta y = (X+ \Delta X) c} \Vert \Delta X, \Delta y \Vert_F\]该问题实际上为主成分分析问题（PCA），可以参见 最佳低秩逼近与主成分分析，由于条件等价于：
\[\begin{pmatrix} X+\Delta X & y +\Delta y \end{pmatrix} \begin{pmatrix} c \\ -1 \end{pmatrix} =0\]此时，向量$(c,-1)^T$也即最小奇异值（0奇异值）所对应右奇异向量，要对矩阵$(X,y)$进行扰动使其最小奇异值为0，
假设$(X,y) = \sum_i \sigma_i u_i v_i^T$ , 则对应的扰动为$(\Delta X,\Delta y) = -\sigma_n u_i v_i^T$
或者也可以将该问题看作最小化投影后的残差问题，也即$\min \Vert A - QQ^T A\Vert ,A = (X,y),Q^TQ=I$ ,也即PCA
Gauss-Newton Method
对于一般的函数空间，空间中的函数可能不一定能表示成基的线性组合。此时通常使用Gauss-Newton方法，在参数$c_k$的邻域内对残差函数$r(c_k) = y - f(c_k)$做Taylor展开，
\[r(c_{k+1}) = r(c_k) + J (c_{k+1} - c_k) + o(\Vert c_{k+1} - c_k \Vert ^2)\]类似于Newton方法，关于Newton方法，可以参见 非线性方程组的解法
忽略二阶无穷小量,并且希望$r(c_{k+1}) \approx 0$, 考虑最小化$\Vert r(c_k) + J(c_{k+1} -c_k) \Vert$
也即最小二乘问题，用广义逆的形式可以得到其解，也即Guass-Newton方法中的迭代公式，
\[c_{k+1} = c_k - J^{+} r(c_k)\]可见该方法和Newton方法本质上为一类方法，同样也可以有带阻尼版本、伪Newton法版本等变种。
Best Square Approximation
如果将函数空间用一组正交多项式组成的基$p_0,p_1,…p_n$来刻画，基于最小二乘与投影的关系得到
\[\min_{p = \sum_i c_i p_i} \Vert f - p \Vert = \sum_i\langle f,p_i \rangle p_i\]其中正交基满足，
\[\langle p_m ,p_n \rangle = 0,m \ne n\]此处内积可以定义为（本节主要在实数域上考虑），
\[\langle f,g \rangle = \int_{-\infty}^{\infty} \rho(x) f(x) g(x) dx\]Orthogonal Basis
正交多项式有很多重要的性质:
- 三项递推性质，$xp_n(x) = \alpha_n p_{n-1}(x) + \gamma_n p_n(x) + \beta_n p_{n+1}(x)$
- $p_n(x)$有$n$个不同的实根
- 根的交错性质，$p_{n-1}(\mu_i) = 0, p_n(\lambda_i) = 0$, 则有$\lambda_i \le \mu_i \le \lambda_{i+1}$
Three-Term Recurrence
下面逐一进行证明，证明的关键是写成矩阵形式:
首先，$p(x) \rightarrow xp(x)$ 满足线性映射，按定义验证即可，因此可以用表示矩阵$A$在多项式基$p_0,p_1,…,p_n$下表示。
此时对比系数可以得到，$H$是一个上Hessenberg矩阵，$xp(x)$仅将多项式的系数抬升1，最后加多一项新的基$p_{n+1}$
\[\begin{align} (xp_0(x),...,xp_n(x)) & = A(p_0,...,p_n) =(p_0,...,p_n) H + \beta_n p_{n+1} \end{align}\]由于$\langle xp(x),q(x) \rangle = \langle p(x),xq(x) \rangle$, 因此$A$满足对称性，也即$A = A^T$
又记$V = (p_0,…,p_n)$，则根据$V^T AV = V^TV H + V^T \beta_n p_{n+1} = H$
因此，$H = H^T$，此时的$H$退化为三对角矩阵（非零元仅在主对角线和副对角线上）
此时，也即得到了三项递推公式，
\[xp_n(x) = \alpha_n p_{n-1}(x) + \gamma_n p_n(x) + \beta_n p_{n+1}(x)\]Real Root
对于$p_n(x)$的根，因为$p_{n+1}(x) =0$，代入得到，
\[x(p_0(x),...,p_{n-1}(x)) =(p_0(x),...,p_{n-1}(x)) H\]此时，等价于求解矩阵$H$的左特征值与左特征向量，$\lambda y^T = y^T H$
求解$\lambda$则得到了$p_n(x)$的根，由于$H$为实对称矩阵，其有$n$个不同的实特征值
Interlacing Property
根据上式，$x(p_0(x),…,p_{n-1}(x)) =(p_0(x),…,p_{n-1}(x)) H_n$
对于$p_{n-1},p_n$的根，满足$p_{n-1}(\mu_i) = 0, p_n(\lambda_i) = 0$
实际上为$H_{n-1},H_n$的特征值，根据Cauthy交错定理，该定理可以参见 特征值的变分性质和谱聚类 中的对应部分。
则根也满足交错性质，$\lambda_i \le \mu_i \le \lambda_{i+1}$
Best Uniform Approximation
上一节的内容为最小二乘拟合，也即最佳均方逼近，本节讨论最佳一致逼近。最佳均方逼近可以看作在$L_2$范数下的逼近，而最佳一致逼近则是在$L_{\infty}$范数意义下的逼近，要求最大误差最小化，也即$\min \max \Vert f - p \Vert = \min \Vert f-p\Vert_{\infty}$，其中$p_n$为$n$阶多项式。
Existence
我们首先关心最佳一致逼近问题解的存在性，可以使用紧集上的连续函数比有界证明。
将逼近多项式$p$表示为多项式基的线性组合，其组合系数为$c$，首先考虑$\Vert c \Vert_2=1$的集合$\tilde p$，此时定义域为有界闭集，存在最大最小值。
\[m \le \min_{\Vert c \Vert_2=1} \Vert f - p \Vert_{\infty} \le M\]假设$p^{\star}$为最佳一致逼近，其一定可以表示为$ p = \alpha \tilde p$,下面证明对于最佳一致逼近比例系数$\alpha$也是有界的。
\[\begin{align} \Vert f - p^{\star} \Vert_{\infty} &= \Vert f - \alpha \tilde p\Vert_{\infty} \\ & \ge \Vert \alpha \tilde p \Vert_{\infty}- \Vert f \Vert_{\infty} \\ & \ge \vert \alpha m\vert - \Vert f \Vert_{\infty} \\ & \gt \Vert f \Vert_{\infty} (\text{ Let } \vert \alpha m\vert > \Vert f \Vert_{\infty}) \\ &= \Vert f - 0 \Vert_{\infty}\\ & \ge \Vert f - p^{\star} \Vert_{\infty} \end{align}\]由上式可以看出，最佳一致逼近应该满足$\vert \alpha \vert \le \frac{\Vert f \Vert_{\infty}}{m} $ ,也即定义域为$\Vert c \Vert_2 \le \frac{\Vert f \Vert_{\infty}}{m}$.
据此，最佳一致逼近定义在紧集上，存在最小值，也即存在所求的解。
De la Vallee-Poussin
De la Vallee-Poussin定理说明的是：给定多形式$q$满足函数$f-q$中存在一个长度为$n+1$正负交错的点组（或称非均匀交错点组，Nonuniform Alternating Set）$x_0,x_1,…x_{n+1}$，则可以得到最佳一致逼近的一个下界，
\[\min \Vert f(x) - p_n(x) \Vert \ge \min_i \vert f(x_i) - q(x_i) \vert\]证明的思路是反证法，假设最佳一致逼近多项式为$p^{\star}$，且满足
\[\min \Vert f(x) - p^{\star}(x) \Vert < \min_i \vert f(x_i) - q(x_i) \vert\]那么对于$n$阶多项式$p^{\star} - q$，在$x_0,x_1,…,x_{n+1}$这$n+2$个点的符号由$q(x_i)$所控制，因此多项式$p^{\star} - q$也存在长度为$n+2$的正负交错的点组，则其必有$n+1$个根，只能为0多项式，矛盾。
Chebyshev Theorem
Chebyshev定理是最佳一致逼近理论中非常深刻的结论，其给定了Chebyshev交错点组与最佳一致逼近的关系，其中Chebyshev交错点组定义为正负交错且达到函数范数的点组，也即满足$\vert f(x_i) -q(x_i) \vert = \Vert f - q \Vert_{\infty}$
此时，$p$为最佳一致逼近多项式的充要条件是，$f-p$存在长度为$n+2$的Chebyshev交错点组。
Chebyshev定理充分性的证明可以通过De la Vallee-Poussin定理得到，由于$f-p$存在长度为$n+2$的Chebyshev交错点组,因此有：
\[\Vert f(x) - p^{\star}(x) \Vert_{\infty} < \min_i \vert f(x_i) - p(x_i) \vert = \Vert f(x) - p(x) \Vert_{\infty}\]因此，此时的$p$即为最佳一致逼近$p^{\star}$.
而必要性的证明，使用反证法，假设$p$为最佳一致逼近多项式，但$f-p$的Chebyshev的交错点组$x_0,x_1,…x_{n+1}$的长度不超过$n+1$.
取每段区间$[x_i,x_{i+1}]$内最靠右的零点$t_i$，则可以得到$n$个零点，构造多项式$q = \alpha \prod_i (x -t_i)$ ,令每段区间$[x_i,x_{i+1}]$中最大值为$M_i$，次大值为$m_i (m_i < M_i)$，控制$\alpha$的大小使得$q$满足，$\Vert q \Vert_{\infty} = \min_i \frac{M_i-m_i}{2}$, 则可以发现多项式$q+p$为一个更优的逼近多项式，因为
\[\Vert f - (q+p) \Vert_{\infty} < \Vert f - p \Vert_{\infty}\]综上，Chebyshev定理给出了最佳一致逼近的充要条件。
应用在插值上，Chebyshev定理根据Chebyshev插值结点，可以给出插值的误差上界，参见 函数的插值
对于问题$\min \Vert \prod_i (x-x_i) \Vert_{\infty}$ ,也可以看作寻找$n-1$次多形式逼近首一多项式$x^n$, 也即$\min \Vert x^n - p_{n-1}(x) \Vert_{\infty}$.
而Chebyshev多项式给出了$n+1$个交错点组，因此为该问题上的最佳一致逼近多项式。
Chebyshev定理在共轭梯度法的收敛界、误差界的证明中同样起到了重要的作用。
Uniqueness
上面已经证明了最佳一致逼近的存在性，而根据Chebyshev定理还可以证明最佳一致逼近多项式的唯一性。
假设$p_1,p_2$均为最佳一致逼近多项式，令$p_0 = \frac{p_1 + p_2}{2}$，记$\Vert f - p^{\star} \Vert_{\infty} =M$
\[\begin{align} M &\le \Vert f- p_0 \Vert_{\infty} \\ &\le \Vert f- \frac{p_1+p_2}{2} \Vert_{\infty} \\ &\le \frac{1}{2} \Vert f- p_1 \Vert_{\infty} + \frac{1}{2} \Vert f- p_2 \Vert_{\infty} \\ &= \frac{1}{2} M + \frac{1}{2} M \\ &= M \end{align}\]因此，$p_0$也为最佳一致逼近多项式，根据Chebyshev定理其存在$n+2$个Chebyshev交错点组，在点组上满足，
\[p_0(x_i) = p_1(x_i) = p_2(x_i) = (-1)^i M\]而$p_1,p_2$为$n$阶多项式，在$n+2$个点上重合，也即$n$阶多项式$p_1 - p_2$有$n+2$个根，只能有$p_1 = p_2$。
因此，最佳一致逼近不仅存在，而且是唯一的。
Remez Algorithm
由于最佳一致逼近问题较为难解，Remez 算法是一种求解最佳一致逼近的近似算法，基于Chebyshev一致逼近定理。
算法流程如下：
- 初始化点组$x_0,x_1,…,x_{n+1}$，也可以用Chebyshev交错点组初始化，理由是Chebyshev插值函数的一致逼近性质可能也较好
- 解方程组$p(x_i) - f(x_i) = (-1)^i E$,得到多项式$p(x_i)$与一致逼近界$E$
- 用$p-f$的极值点替换点组$x_0,x_1,…,x_{n+1}$中的一部分点，理由是最佳一致逼近针对极值进行优化
- 迭代进行，直到$\Vert p(x_i) - f(x_i) \Vert_{\infty} \approx E$
此时$E$可以作为最佳一致逼近的一个近似解，且根据De la Vallee-Poussin定理$E$一定是一个下界。