上節(jié)課主要介紹了特征值與特征向量的概念凡桥，低階矩陣的特征值可以通過列出特征方程求解蟀伸，高階矩陣則可以通過冪法與反冪法迭代求解出最大特征值與最小特征值（模），要求出矩陣的全部特征值則需要借助矩陣的 QR分解來將矩陣相似化為一個上三角矩陣缅刽，相似化過程不改變矩陣的特征值啊掏，因此轉化后的上三角矩陣的對角線元素即為原矩陣的特征值。本節(jié)課主要介紹曲線擬合（機器學習中的線性回歸算法）與插值算法衰猛。

1. 擬合與插值（Fitting, Interpolation）

• 曲線擬合是給定一個數(shù)據(jù)集迟蜜，構建一個數(shù)學方程能夠最好的匹配這些數(shù)據(jù)點
• 插值則是構建一個多項式方程來完美地契合給定的數(shù)據(jù)集，當數(shù)據(jù)點比較少的時候啡省，一個簡單的多項式方程就可以實現(xiàn)插值娜睛，而數(shù)據(jù)點比較多的時候髓霞，更多地采用的是樣條插值，在各個區(qū)間構建不同的多項式方程

Fitting and Interpolation.png

2. 線性曲線擬合

線性曲線擬合是使用如下線性方程（一元一階多項式畦戒，對應機器學習線性回歸算法中只有一個特征的情況）來擬合給定的數(shù)據(jù)集： 若數(shù)據(jù)集中只有兩個點方库，則該線性方程可以準確求解出來，通過這兩個點障斋；當數(shù)據(jù)集超過兩個點時薪捍，該線性方程無法通過所有的數(shù)據(jù)點，此時需要找出擬合數(shù)據(jù)點最好的那個方程配喳。

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 常見的幾種非線性方程

（這里同線性擬合類似识樱，只考慮單變量的情況，對應線性回歸算法中的一個特征）： 對非線性擬合震束，很容易將其改寫為線性形式： 將數(shù)據(jù)集進行相應的特征變換后怜庸，采用線性最小二乘法求解相應的參數(shù)。然而驴一，對于給定的數(shù)據(jù)集休雌，如何判斷使用哪一種非線性方程進行擬合呢？一個可行的方案是首先通過繪制數(shù)據(jù)集中的點（當特征很多時肝断，很難可視化）杈曲，確定是選擇線性還是非線性進行擬合驰凛，如果選擇非線性擬合，具體選擇哪一種非線性方程則要取決于其背后的物理現(xiàn)象和數(shù)學原理担扑。同樣的恰响，可以通過以特定的方式繪制數(shù)據(jù)點并檢查這些點是否符合直線來預測所提出的非線性函數(shù)是否具有很好的擬合能力。
MATLAB實現(xiàn)：

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 二次及高階多項式擬合

多項式是指具有如下形式的函數(shù)（單變量）： 稱作多項式的階數(shù)涌献。如果給定數(shù)據(jù)集中有個點胚宦，則可以確定個從1階到階的擬合多項式，其中階多項式可以實現(xiàn)通過每一個數(shù)據(jù)點燕垃，但并不是多項式的階數(shù)越高越好枢劝，這取決于實際問題，當數(shù)據(jù)集點本身就具有噪聲污染時卜壕，提高擬合的階數(shù)并沒有意義您旁。（同機器學習中的過擬合，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 函數(shù) 返回在 x 處計算的 n 次多項式的值秦驯。輸入?yún)?shù) p 是長為 n+1 的向量尺碰，
其元素是按要計算的多項式降冪排序的系數(shù)。

    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

上述多項式擬合可以擴展到更一般的情況： 定義總誤差為 令偏導為 改為線性系統(tǒng)形式译隘，有 矩陣形式為
MATLAB實現(xiàn)：

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. 多項式插值

如前所述葱蝗，插值是讓擬合的曲線通過所有給定的數(shù)據(jù)點（總誤差為0，機器學習中的為）细燎，對個點两曼，總是能找到一個階的多項式函數(shù)通過所有的點： 矩陣形式

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

image.png

考慮兩個點, ，其一階拉格朗日多項式具有如下形式： 令該函數(shù)通過這兩個點玻驻，得 因此多項式為： 類似地悼凑，對三個點，可得其二階拉格朗日具有形式 因此璧瞬，拉格朗日多項式的通解形式為： MATLAB實現(xiàn)：

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 牛頓多項式插值

上述拉格朗日多項式插值有一個缺陷户辫，當數(shù)據(jù)集增加一個插值點時，每個系數(shù)都需要重新計算嗤锉，而牛頓多項式插值則不需要重新計算：

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 線性樣條插值

當數(shù)據(jù)點比較少時垫挨，使用低階多項式就可以實現(xiàn)插值韩肝，然而當數(shù)據(jù)點比較多時，插值的多項式需要很高的階數(shù)九榔，這會造成插值點不準確（同回歸中預測不準確）哀峻。因此，這個時候可以采用區(qū)間插值的方法哲泊，即每兩個或三個相鄰點之間剩蟀，構建具有不同參數(shù)的插值多項式（階數(shù)一致）。

線性插值就是每兩個相鄰的點之間構建各自的線性函數(shù)進行插值切威。此時喻旷，第一個點與第二個點之間的插值函數(shù)(拉格朗日插值公式)為：，同樣地牢屋，很容易得到第個點與第個點之間的插值函數(shù)為：

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 二次樣條插值

給定個數(shù)據(jù)點遍尺，二次樣條插值是在個區(qū)間內(nèi)各自構建一個二次曲線插值數(shù)據(jù)： 可以看到每個多項式有個待定系數(shù)，共有個待定系數(shù)，因此求解規(guī)則包括：

• 每個二次曲線在區(qū)間兩端點與的值為： 共構建出個方程
• 每兩個相鄰區(qū)間的節(jié)點處耳幢，對應多項式的斜率相同（保證導數(shù)連續(xù)）： 共構建個方程
• 上述條件共構建了個方程斑粱，想要求解個參數(shù)，還少一個條件鼓择，因此三幻，通常假設第一個點或者最后一個點的二階導數(shù)為零： 即第一個或最后一個多項式為線性多項式
  
  image.png

5.3 三次樣條插值

對個數(shù)據(jù)點，三次樣條插值是在個區(qū)間內(nèi)各自構建一個三次多項式進行插值呐能。三次樣條插值包含標準形式與拉格朗日形式念搬。考慮如下標準形式： 共有個待定參數(shù)摆出。求解規(guī)則：

• 每個三次曲線在區(qū)間兩端點與的值為： 共構建出個方程
• 每兩個相鄰區(qū)間的節(jié)點處朗徊，對應多項式的斜率相同（保證一階導數(shù)連續(xù)）： 共構建個方程
• 每兩個相鄰區(qū)間的節(jié)點處，對應多項式的二階導數(shù)相同（當從一條曲線轉換為下一條曲線時偎漫，斜率的變化是連續(xù)的）： 共構建個方程
• 上述條件共構建了個方程爷恳，還差兩個條件，通常要求第一個點和最后一個點的二階導數(shù)為零(自然三次樣條插值象踊，natural cubic splines)：

image.png

• 從三次多項式的二階導數(shù)出發(fā)，三次多項式的二階導數(shù)是一個線性函數(shù)节芥，采用拉格朗日多項式插值形式（點）： 積分可得： 再次積分在刺，有 令 解得
  
  因此 定義 ，則 進一步利用 头镊，可得： 注意是待定數(shù)蚣驼，然而這里只有個方程，還需要兩個條件相艇，利用標準形式中的最后一個條件颖杏，即第一個點和最后一個點的二階導數(shù)為零，可得坛芽。

MATLAB實現(xiàn)：

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. 總結

本節(jié)課主要介紹了曲線的擬合與插值留储，與機器學習中的回歸算法相似。給定一組數(shù)據(jù)集咙轩，擬合是找到最好的曲線來匹配數(shù)據(jù)集中的點（線性回歸获讳，非線性回歸），插值則是找到通過所有數(shù)據(jù)點的曲線活喊。線性回歸是最簡單的曲線擬合算法丐膝，非線性回歸可以通過特征變換轉化為線性回歸，而其中的高階多項式擬合還可以進一步擴展到更一般的情形钾菊，求解方法均是基于最優(yōu)點處的導數(shù)為零得到帅矗。對于插值，任意個點都可以找到一個階的多項式函數(shù)通過所有數(shù)據(jù)點煞烫，求解方法有拉格朗日法和牛頓法浑此， 拉格朗日法計算簡單，但增加數(shù)據(jù)點時滞详，所有參數(shù)均需重新計算尤勋，而牛頓法則不需要重新計算，但計算較為復雜茵宪。另一方面最冰，當數(shù)據(jù)點比較少時，多項式插值較為簡單稀火，但對于數(shù)據(jù)點較多的情況暖哨，多項式的階數(shù)需要很高，計算復雜，并且插值時并不可靠（過擬合）篇裁。樣條插值是在每個區(qū)間內(nèi)構建各自的多項式插值函數(shù)沛慢，一般有線性樣條插值、二次樣條插值以及三次樣條插值达布。其中為求解簡便团甲，三次樣條插值除了標準形式外，還有一個便于求解的拉格朗日形式黍聂。