上节课主要介绍了特征值与特征向量的概念，低阶矩阵的特征值可以通过列出特征方程求解，高阶矩阵则可以通过幂法与反幂法迭代求解出最大特征值与最小特征值（模），要求出矩阵的全部特征值则需要借助矩阵的 QR分解来将矩阵相似化为一个上三角矩阵，相似化过程不改变矩阵的特征值，因此转化后的上三角矩阵的对角线元素即为原矩阵的特征值。本节课主要介绍曲线拟合（机器学习中的线性回归算法）与插值算法。

1. 拟合与插值（Fitting, Interpolation）

• 曲线拟合是给定一个数据集，构建一个数学方程能够最好的匹配这些数据点
• 插值则是构建一个多项式方程来完美地契合给定的数据集，当数据点比较少的时候，一个简单的多项式方程就可以实现插值，而数据点比较多的时候，更多地采用的是样条插值，在各个区间构建不同的多项式方程

Fitting and Interpolation.png

2. 线性曲线拟合

线性曲线拟合是使用如下线性方程（一元一阶多项式，对应机器学习线性回归算法中只有一个特征的情况）来拟合给定的数据集： 若数据集中只有两个点，则该线性方程可以准确求解出来，通过这两个点；当数据集超过两个点时，该线性方程无法通过所有的数据点，此时需要找出拟合数据点最好的那个方程。

image.png

image.png

image.png

image.png

function [a1,a0] = LinearRegression(x, y)
% LinearRegration calculates the coefficients a1 and a0 of the linear
% equation y = a1*x + a0 that best fit n data points.
% Input variables:
% x    A vector with the coordinates x of the data points.
% y    A vector with the coordinates y of the data points.
% Output variable:
% a1   The coefficient a1.
% a0   The coefficient a0.

nx = length(x);
ny = length(y);
if nx ~= ny
    disp('ERROR: The number of elements in x must be the same as in y.')
    a1 = 'Error';
    a0 = 'Error';
else
Sx = sum(x);
Sy = sum(y);
Sxy = sum(x.*y);
Sxx = sum(x.^2);
a1 = (nx*Sxy - Sx*Sy)/(nx*Sxx - Sx^2);
a0 = (Sxx*Sy - Sxy*Sx)/(nx*Sxx - Sx^2);
end

3. 非线性曲线拟合

image.png

3.1 常见的几种非线性方程

（这里同线性拟合类似，只考虑单变量的情况，对应线性回归算法中的一个特征）： 对非线性拟合，很容易将其改写为线性形式： 将数据集进行相应的特征变换后，采用线性最小二乘法求解相应的参数。然而，对于给定的数据集，如何判断使用哪一种非线性方程进行拟合呢？一个可行的方案是首先通过绘制数据集中的点（当特征很多时，很难可视化），确定是选择线性还是非线性进行拟合，如果选择非线性拟合，具体选择哪一种非线性方程则要取决于其背后的物理现象和数学原理。同样的，可以通过以特定的方式绘制数据点并检查这些点是否符合直线来预测所提出的非线性函数是否具有很好的拟合能力。
MATLAB实现：

texp = 2:2:30;
vexp = [9.7 8.1 6.6 5.1 4.4 3.7 2.8 2.4 2.0 1.6 1.4 1.1 0.85 0.69 0.6];
vexpLOG = log(vexp);
R = 5E6;
[a1,a0] = LinearRegression(texp, vexpLOG)
b = exp(a0)
C = -1/(R*a1)
t = 0:0.5:30;
v = b*exp(a1*t);
plot(t,v,texp,vexp,'or')

image.png

3.2 二次及高阶多项式拟合

多项式是指具有如下形式的函数（单变量）： 称作多项式的阶数。如果给定数据集中有个点，则可以确定个从1阶到阶的拟合多项式，其中阶多项式可以实现通过每一个数据点，但并不是多项式的阶数越高越好，这取决于实际问题，当数据集点本身就具有噪声污染时，提高拟合的阶数并没有意义。（同机器学习中的过拟合，overfitting类似）

image.png

clear, clc
x = 0:0.4:6;
%x=1:4
y = [0 3 4.5 5.8 5.9 5.8 6.2 7.4 9.6 15.6 20.7 26.7 31.1 35.6 39.3 41.5];
n = length(x);
m = 4;
for i = 1:2*m
    xsum(i) = sum(x.^(i));
end
% Beginning of Step 3
a(1,1) = n;
b(1,1)= sum(y);
for j= 2:m+1
    a(1,j)=xsum(j-1);
end
for i = 2:m+1
    for j=1:m+1
        a(i,j)= xsum(j+i-2);
    end
    b(i,1) = sum(x.^(i-1).*y);
end
% Step 4
p=(a\b)'
for i = 1:m+1
    Pcoef(i) = p(m+2-i);
end
epsilon=0:0.1:6;
stressfit=polyval(Pcoef,epsilon);
plot(x,y,'ro',epsilon,stressfit,'k','linewidth',2)
xlabel('Strain','fontsize',20)
ylabel('Stress (MPa)','fontsize',20)
%%%%%%%%%%%%%%%%%%%%
>> help polyval
polyval - 多项式计算

  此 MATLAB 函数 返回在 x 处计算的 n 次多项式的值。输入参数 p 是长为 n+1 的向量，
其元素是按要计算的多项式降幂排序的系数。

    y = polyval(p,x)
    [y,delta] = polyval(p,x,S)
    y = polyval(p,x,[],mu)
    [y,delta] = polyval(p,x,S,mu)

    See also polyder, polyfit, polyint, polyvalm

    Reference page for polyval
    Other functions named polyval

image.png

上述多项式拟合可以扩展到更一般的情况： 定义总误差为 令偏导为 改为线性系统形式，有 矩阵形式为
MATLAB实现：

clear all
x = [0.6 0.8 0.85 0.95 1.0 1.1 1.2 1.3 1.45 1.6 1.8];
y = [0.08 0.06 0.07 0.07 0.07 0.06 0.06 0.06 0.05 0.05 0.04];
a(1,1) = sum(1./x.^2);
a(1,2) = sum(exp(-2*x.^2)./x.^2);
a(2,1) = a(1,2);
a(2,2) = sum(exp(-4*x.^2)./x.^2);
b(1,1) = sum(y./x);
b(2,1) = sum((y.*exp(-2*x.^2))./x);
AB = a\b
xfit = 0.6:0.02:1.8;
yfit = AB(1)./xfit + AB(2)*exp(-2*xfit.^2)./xfit;
plot(x,y,'o',xfit,yfit)

image.png

4. 多项式插值

如前所述，插值是让拟合的曲线通过所有给定的数据点（总误差为0，机器学习中的为），对个点，总是能找到一个阶的多项式函数通过所有的点： 矩阵形式

4.1 多项式插值的拉格朗日形式：

image.png

考虑两个点, ，其一阶拉格朗日多项式具有如下形式： 令该函数通过这两个点，得 因此多项式为： 类似地，对三个点，可得其二阶拉格朗日具有形式 因此，拉格朗日多项式的通解形式为： MATLAB实现：

function Yint = LagrangeINT(x,y,Xint)
% LagrangeINT fits a Lagrange polynomial to a set of given points and
% uses the polynomial to determines the interpolated value of a point.
% Input variables:
% x  A vector with the x coordinates of the given points.
% y  A vector with the y coordinates of the given points.
% Xint  The x coordinate of the point to be interpolated.
% Output variable:
% Yint  The interpolated value of Xint.

n = length(x);
for i = 1:n
    L(i) = 1;
    for j = 1:n
        if j ~= i
            L(i)= L(i)*(Xint-x(j))/(x(i)-x(j));
        end
    end
end
y.*L
Yint = sum(y.*L);
%%%%%%%%
>> x = [1 2 4 5 7];
>> y = [52 5 -5 -40 10];
>> Yinterpolated = LagrangeINT(x,y,3)

ans =

   -5.7778    2.6667   -4.4444   13.3333    0.2222

Yinterpolated =

    6.0000

>> i=1;
>> for x1=0:0.1:10
y1(i)=LagrangeINT(x,y,x1);
i = i+1;
end
>> x1 = 0:0.1:10;
>> plot(x,y,'*',x1,y1)

image.png

4.2 牛顿多项式插值

上述拉格朗日多项式插值有一个缺陷，当数据集增加一个插值点时，每个系数都需要重新计算，而牛顿多项式插值则不需要重新计算：

image.png

image.png

function Yint = NewtonsINT(x,y,Xint)
% NewtonsINT fits a Newtons polynomial to a set of given points and
% uses the polynomial to determines the interpolated value of a point.
% Input variables:
% x  A vector with the x coordinates of the given points.
% y  A vector with the y coordinates of the given points.
% Xint  The x coordinate of the point to be interpolated.
% Output variable:
% Yint  The interpolated value of Xint.

n = length(x);
a(1) = y(1);
for i = 1:n-1
    divDIF(i,1) = (y(i+1)-y(i))/(x(i+1)-x(i));
end
for j = 2:n-1
    for i = 1:n-j
        divDIF(i,j) = (divDIF(i+1,j-1) - divDIF(i,j-1))/(x(j+i) - x(i));
    end
end
for j=2:n
    a(j)=divDIF(1,j-1);
end
Yint=a(1);
xn=1;
for k=2:n
    xn=xn*(Xint-x(k-1));
    Yint=Yint+a(k)*xn;
end

5. 样条插值

5.1 线性样条插值

当数据点比较少时，使用低阶多项式就可以实现插值，然而当数据点比较多时，插值的多项式需要很高的阶数，这会造成插值点不准确（同回归中预测不准确）。因此，这个时候可以采用区间插值的方法，即每两个或三个相邻点之间，构建具有不同参数的插值多项式（阶数一致）。

线性插值就是每两个相邻的点之间构建各自的线性函数进行插值。此时，第一个点与第二个点之间的插值函数(拉格朗日插值公式)为：，同样地，很容易得到第个点与第个点之间的插值函数为：

image.png

function Yint = LinearSpline(x, y, Xint)
% LinearSpline calculates interpolation using linear splines.
% Input variables:
% x    A vector with the coordinates x of the data points.
% y    A vector with the coordinates y of the data points.
% Xint The x coordinate of the interpolated point. 
% Output variable:
% Yint The y value of the interpolated point.

n = length(x);
for i = 2:n
    if Xint < x(i)
        break
    end
end
Yint = (Xint - x(i))*y(i-1)/(x(i-1)-x(i)) + (Xint - x(i-1))*y(i)/(x(i)-x(i-1));

5.2 二次样条插值

给定个数据点，二次样条插值是在个区间内各自构建一个二次曲线插值数据： 可以看到每个多项式有个待定系数，共有个待定系数，因此求解规则包括：

• 每个二次曲线在区间两端点与的值为： 共构建出个方程
• 每两个相邻区间的节点处，对应多项式的斜率相同（保证导数连续）： 共构建个方程
• 上述条件共构建了个方程，想要求解个参数，还少一个条件，因此，通常假设第一个点或者最后一个点的二阶导数为零： 即第一个或最后一个多项式为线性多项式
  
  image.png

5.3 三次样条插值

对个数据点，三次样条插值是在个区间内各自构建一个三次多项式进行插值。三次样条插值包含标准形式与拉格朗日形式。考虑如下标准形式： 共有个待定参数。求解规则：

• 每个三次曲线在区间两端点与的值为： 共构建出个方程
• 每两个相邻区间的节点处，对应多项式的斜率相同（保证一阶导数连续）： 共构建个方程
• 每两个相邻区间的节点处，对应多项式的二阶导数相同（当从一条曲线转换为下一条曲线时，斜率的变化是连续的）： 共构建个方程
• 上述条件共构建了个方程，还差两个条件，通常要求第一个点和最后一个点的二阶导数为零(自然三次样条插值，natural cubic splines)：

image.png

• 从三次多项式的二阶导数出发，三次多项式的二阶导数是一个线性函数，采用拉格朗日多项式插值形式（点）： 积分可得： 再次积分，有 令 解得
  
  因此 定义 ，则 进一步利用 ，可得： 注意是待定数，然而这里只有个方程，还需要两个条件，利用标准形式中的最后一个条件，即第一个点和最后一个点的二阶导数为零，可得。

MATLAB实现：

function Yint = CubicSplines(x,y,xint)
n = length(x);
a1=0;
an=0;
A = zeros(n-2,n-2);
b = zeros(n-2,1);
for i=1:n-2
    h(i) = x(i+1)-x(i);
    h(i+1) = x(i+2)-x(i+1);
    if i == 1
       A(i,i) = 2*(h(i)+h(i+1));
       A(i,i+1) =  h(i+1);
       b(i,1) = 6*((y(i+2)-y(i+1))/h(i+1)-(y(i+1)-y(i))/h(i));
    elseif i == n-2
        A(i,i-1) = h(i);
        A(i,i) = 2*(h(i)+h(i+1));
        b(i,1) = 6*((y(i+2)-y(i+1))/h(i+1)-(y(i+1)-y(i))/h(i));
    else
        A(i,i-1) = h(i); 
        A(i,i) = 2*(h(i)+h(i+1));
        A(i,i+1) = h(i+1);
        b(i,1) = 6*((y(i+2)-y(i+1))/h(i+1)-(y(i+1)-y(i))/h(i));
    end
end
 p = A\b;
 a(1) = 0;
 a(n) = 0;
 for i=1:n-2
    a(1+i) = p(i);
 end
 
 for i = 2:n
    if xint <= x(i)
        break
    end
 end
  h(i) = x(i)-x(i-1);
  Yint = a(i-1)/(6*h(i))*(x(i)-xint)^3+a(i)/(6*h(i))*(xint-x(i-1))^3+(y(i-1)/h(i)-a(i-1)*h(i)/6)*(x(i)-xint)+(y(i)/h(i)-a(i)*h(i)/6)*(xint-x(i-1));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x = [5 7.5 10 12.5 15 17.5 20 22.5 25 27.5 30];
h = [3.3 7.5 41.8 51.8 61 101.1 132.9 145.5 171.4 225.8 260.9];


xint = 5:0.5:30;
l = length(xint);
yint = zeros(l,1);
for i = 1:l
    yint(i) = CubicSplines(x,h,xint(i));
end
figure
subplot(1,2,1)
plot(x,h,'k*',xint,yint,'r-')
title('cubic')

p = interp1(x,h,xint,'spline');
subplot(1,2,2)
plot(x,h,'k*',xint,p,'r-')
title('spline')

image.png

6. 总结

本节课主要介绍了曲线的拟合与插值，与机器学习中的回归算法相似。给定一组数据集，拟合是找到最好的曲线来匹配数据集中的点（线性回归，非线性回归），插值则是找到通过所有数据点的曲线。线性回归是最简单的曲线拟合算法，非线性回归可以通过特征变换转化为线性回归，而其中的高阶多项式拟合还可以进一步扩展到更一般的情形，求解方法均是基于最优点处的导数为零得到。对于插值，任意个点都可以找到一个阶的多项式函数通过所有数据点，求解方法有拉格朗日法和牛顿法， 拉格朗日法计算简单，但增加数据点时，所有参数均需重新计算，而牛顿法则不需要重新计算，但计算较为复杂。另一方面，当数据点比较少时，多项式插值较为简单，但对于数据点较多的情况，多项式的阶数需要很高，计算复杂，并且插值时并不可靠（过拟合）。样条插值是在每个区间内构建各自的多项式插值函数，一般有线性样条插值、二次样条插值以及三次样条插值。其中为求解简便，三次样条插值除了标准形式外，还有一个便于求解的拉格朗日形式。