自定义一组基函数的最佳平方逼近——matlab实现

作者：老李
时间：2019-10-26

目标：

给出一组非线性相关的函数作为基底，通过有限点集$(x_k,y_k)$，构造拟合函数

过程：

我们设这个基底为 { f i ( x ) } , i = 1 , 2 , . . . M \left \{f _{i}(x)\right \},i= 1,2,...M {fi​(x)},i=1,2,...M
我们的目标是,对于每一个$x_k$寻找一系列$c_i$,使得$E={\sum }_{k=1}^n(c_i{f}_i(x_k)-y_k{)}^2$最小。

我们使用最一般的方法：将E对每一个参数$c_i$进行求导，然后令其导函数为零，最后求出其值。

$$\frac{\partial E}{\partial c_i}=\sum _{k=1}^N((\sum _{j=1}^M{c}_j{f}_j(x_k))-y_k)f_i(x_k)=0,i=1,2,...M\text{\hspace{1em}(1)}$$

为了实现矩阵化操作，我们想要将所求的系数$c_i$分离出来，于是我们改变求和的次序，可得：$$\sum _{j=1}^M(\sum _{k=1}^N{f}_i(x_k)f_j(x_k))c_j=\sum _{k=1}^N{f}_i(x_k)y_k\text{\hspace{1em}(2)}$$

最后我们对其实现矩阵化操作：
我们令$$F=\left[\begin{array}{cccc}{f}_1(x_1)& f_2(x_1)& \cdots & f_M(x_1)\\ f_1(x_2)& f_2(x_2)& \cdots & f_M(x_2)\\ ⋮& ⋮& \ddots & ⋮\\ f_1(x_N)& f_2(x_N)& \cdots & f_M(x_N)\end{array}\right]$$可以将（2）式用矩阵表达出来$$F^{\prime }FC=F^{\prime }Y$$
求得的系数矩阵$C$为$$C=(F^{\prime }F{)}^{-1}{F}^{\prime }Y$$
我用matlab实现了这个过程，代码如下：
吐槽一点，csdn粘贴matlab代码挺麻烦的

function [ C ,err] = leastSquareApprox( X,Y,fun)
%leastSquareApprox: least square approximation by linear combination of nonlinear correlation equations
%   Input: X is the 1*n abscissa vector
%        : Y is the 1*n ordinate vector
%        : fun is the n*1 based function
%  Output: C is the coefficient list for the function
%        : err is the maximum of the least square approximation for each point
if  (length(X) ~= length(Y))
    error('wrong size');
else
    n = length(X);
    X0 = reshape(X,n,1);
    F = subs(fun, X0);
    C = (pinv(F'*F))*F'*Y;
    
    Xt = min(X):0.1:max(X);
    ft = subs(fun, Xt')*C;
    
    scatter(X,Y);
    hold on;
    plot(Xt,ft,'r');
    hold off;
    f0 = subs(fun, X)*C;
    err = max((f0-Y).^2);
end
end

这里需要特别注意的一点是输入的fun这个参数，是我们作为最小二乘拟合(最佳平方逼近)的函数。

这里我来展示一下怎么去实现这个基函数。
这里，我的基函数通过符号变量来表达，
我们可以选择多项式作为基底，当我们输入记得的数量n时，该基底可以表示为
$$\left\{1,x,x^2,x^3,\cdots ,x^{n-1}\right\}$$
(如果我们不去讨论这组基的具体在所张成的空间的定义，或者说不一定要符合在该空间多定义出的的正交这一性质的话，也就是说我们只要保证基底中的元素不能直接相互表示的话)
我们可以写一段程序来实现它：(syms和sym 用来定义符号变量)

n = input('please input the size of the data ');
syms x;
B = sym(zeros(1,n));
for i = 0:n-1   
    B(i+1) = x.^i;
end

这样，我们就把作为基底的多项式输出出来了。

效果

假设我们只输入五个点：对应的自变量X=[1，2，3，4，5]和因变量Y=[15，64，84，12，34]
这里我先选用3个基底（n=3）来进行拟合，也就是二次函数拟合：
得到的效果如下：

输出的数值为：

其中误差为d

关于误差的讨论

已知i我们输入的为5个点，n为基底中元素的个数

当我们令n=5时，由范德蒙德行列式，我们可以知道，当我们输入的点不同时，该行列式的值不为0，也就是说，插值函数是具有唯一性的。也就是说，无论我们用的是拉格朗日插值还是牛顿插值，还是这里用到的矩阵运算，我们得到的东西是一样的。而不同的插值方法的功效主要还是体现在计算机的运行效率上，比如牛顿插值用到了递归的方法从而节约了运算的时间。

当我们把n变为>5的数，我们可以保证逼近出来的曲线是过每一个点的，所得出来的误差应该为0。但由于插值函数的唯一性，我们不应该认为这是插值，但这样的一个方法应该叫什么，我现在并不清楚，也许以后会学到并理解。

我们这里把n设为8
结果如下：

写在最后

1.在计算机实现逼近的过程中，输入的变量与输出的变量一定是离散且有限的，也就是说，在计算机数值计算中，我们用的大多为$\sum$而不是$\int$。

2.假如我们输入的是N个点，假如我们使用的基底（不能相互表示）中元素的数量是>=N的，那我们逼近出来的曲线一定是过我们输入的每一个点的，（=N的时候是为插值）。

3.同样，我们也可以注意到，其实DFT算法的核心思想与我今天所写的其实是一致的。唯一不同的是他运用的是$e^{-\frac{2\pi im}{N}}$作为基底，其中$m=0,1,2,\cdots ,N-1$。