最小二乗法

補間関数を使って対象関数の近似を考える時、補間の節点が増やしても補間関数の誤差が小さくならないことが分かりました。

今日の授業では、「最小二乗法」という節点における新たな関数近似の方法を解説します。

1. 復習： ルンゲ現象

数の多い節点で多項式による関数の補間を行う時、多項式の次数が高く、補間関数の誤差が大きくなる問題があります。 カール・ルンゲが、ある関数を多項式補間で近似したときの誤差を調べていて発見したので、この現象がルンゲ現象と呼ぼられる。

区間$[-1,1]$で対象関数$f$を考えます。

$$ f(x) = \frac{1}{1+25x^2} $$

節点を以下の等分割を使用します。 $$ x_i = -1+ (i-1)2/n , \quad (i=0,1,2,\cdots 2n+1) $$ 以下のグラフで$n=5$の時の補間関数を描いています。

演習１

$n=6,7,8,9,10$のとき、$f$の補間関数と誤差を確認しなさい。

In [8]:
function value = myf(x)
    value =  1 ./(1+25 .* x.^2);
end

function value = Lagrange(xlist,i,x)
    n = length(xlist);
    denominator = 1;
    numerator = 1;
    for k=1:i-1
        numerator = numerator .* (x-xlist(k));
        denominator = denominator .* (xlist(i) - xlist(k));
    end
    for k=i+1:n
        numerator = numerator .* (x-xlist(k));
        denominator = denominator .* (xlist(i) - xlist(k));
    end
    value=numerator ./ denominator;    
end

n=5;
xlist= -1:2/n:1;

x=-1:0.01:1;
f_interpolation =0;
for k=1:(n+1)
  f_interpolation = f_interpolation + Lagrange(xlist, k, x)*myf(xlist(k));
end

hold on
plot(x, myf(x),'r-')
plot(x,f_interpolation,'-')
grid on
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 -0.2 0 0.2 0.4 0.6 0.8 1 -1 -0.5 0 0.5 1 gnuplot_plot_1a gnuplot_plot_2a

2 最小二乗法のモデル

2.1 関数近似に使用される最小二乗法

与えられる関数$f(x)$に対して、基底$\{\phi_1, \phi_2, \cdots, \phi_n\}$の一次結合で現れる関数$p(x)$を利用して$f$の近似を考えます。 $$ f(x) \approx p(x) = \sum_{k=1}^n c_k \phi_i(x)\:. $$

節点$x_1,\cdots,x_N$における以下の誤差を最小化する問題を考えます。 $$ \min_{p} \sum_{i=1}^N ( p(x_i) - y_i )^2 $$ ここで、$y_i =f(x_i)$.

上記の誤差で考える近似関数の特徴

  • $N$が大きい時、各節点で$p$と$f$が一致しなくでも良いです。$N$点での誤差のバランスを取って、関数の近似を求めます。
  • 大きい$n$を使用するとき、対象関数$f$の良い近似ができます。

実際の計算では、$p_i$を多項式にすることができます。

2.2 実験のデータ処理に使用される最小二乗法

実験で測定した結果の処理には最小二乗法が多く使用されています。実験で測定した結果$(x_i,y_i)$には測定誤差があるので、 測定した結果$(x_i,y_i)$に近い直線または曲線で近似することが要求されます。

以下の例では、赤い線分は理想的なデータの分布（一次関数）であり、青い点は誤差の混入のある測定値です。 与えられる青い点のリストから線分を求める問題があります。

このとき、基底$\{\phi_i\}$を$\{1,x\}$とすると、$p(x)=c_1 + c_2 x$という形を持っています。$p$を求めるために、 以下の最小二乗法のモデルが考えられます。

$$ \min_{p} \sum_{i=1}^N ( p(x_i) - y_i )^2 $$

In [60]:
x=0:0.1:1;
y = 1 +  2.*x + (rand(1,length(x))*0.2 - 0.1);

hold on
plot(x,y,'b+')
f = 1 +  2.*x ;
plot(x,f,'r-')
grid on
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 0.5 1 1.5 2 2.5 3 0 0.2 0.4 0.6 0.8 1 gnuplot_plot_1a gnuplot_plot_2a

3. 最小二乗法の解の計算法

最小化問題を行列の「言葉」で検討してみます。$N\times n$の行列$A$を導入します。

$$ A_{ij}=\phi_j(x_i),~~ (i=1,\cdots,N; j=1,\cdots, n) \:. $$

即ち、

$$ A= \left( \begin{array}{ccccc} \phi_1(x_1) & \phi_2(x_1) & \phi_3(x_1) & \cdots & \phi_n (x_1)\\ \phi_1(x_2) & \phi_2(x_2) & \phi_3(x_2) & \cdots & \phi_n (x_2)\\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \phi_1(x_N) & \phi_2(x_N) & \phi_3(x_N) & \cdots & \phi_n (x_N) \end{array} \right) $$

$N\times 1$ベクトル$b$と$n\times 1$ベクトル$c$を定義します。 $$ b= \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_N \end{array} \right) ,~ c= \left( \begin{array}{c} c_1 \\ c_2 \\ \vdots \\ c_n \end{array} \right) $$

$N\times 1$ ベクトル$P$を定義します。 $$ P_i = p(x_i) $$ このとき、以下の式の変形がわかります。 $$ \sum_{i=1}^N ( p(x_i) - y_i )^2 = (P-b)^T\cdot (P-b) $$

さらに、$Ac=P$を利用して、最小化問題の目的関数が以下の式になることがわかります。

$$ \sum_{i=1}^N ( p(x_i) - y_i )^2 = \|Ac - b\|^2 $$

目的関数の勾配

$g(c):=\|Ac - b\|^2$を定義します。$g(c)$の展開式を確認します。

特に、$c^TA^T b$がスカラーなので、$ c^TA^T b = b^TA c $を注意してください。

$$ g(c)= (Ac - b)^T (Ac - b) = c^TA^TAc - c^TA^T b - b^T Ac + b^Tb = c^TA^TAc - 2 b^T Ac + b^T b $$

よって、$g(c)$の$c$に関する勾配が次に式になります。

$$ \frac{\partial g(c)}{\partial c} = \left( \frac{\partial g(c)}{\partial c_1}, \frac{\partial g(c)}{\partial c_2}, \cdots , \frac{\partial g(c)}{\partial c_n} \right) =2c^TA^TA - 2 b^T A $$

極限値をとる$c$を求める

$g(c)$が最小値を取る時、$g(c)$の$c$に関する勾配がゼロとなります。このとき、

$$ 2c^TA^TA = 2 b^T A,\mbox{ 即ち } A^TAc = A^T b \:. $$

よって、係数ベクトル$c$が以下の式で求められます。

$$ c= (A^TA)^{-1} (A^Tb) \:. $$

演習２

以下の実験の測定結果$(x_i,y_i)$のリストに対して、最小二乗法によって、$p=c_1+c_2x$の形の近似関数を求めてください。求めた近似関数のグラフも描いてください。

In [67]:
data = [0.00000   1.09915
   0.10000   1.17397
   0.20000   1.42845
   0.30000   1.54711
   0.40000   1.70561
   0.50000   1.99224
   0.60000   2.26987
   0.70000   2.44653
   0.80000   2.55409
   0.90000   2.84210
   1.00000   3.04617];
xlist = data(:,1);
ylist = data(:,2);

plot(xlist,ylist,'b+')
grid on
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 1 1.5 2 2.5 3 3.5 0 0.2 0.4 0.6 0.8 1 gnuplot_plot_1a

演習3

$\phi_1,\phi_2,\cdots, \phi_n$を$(n-1)$次までの多項式関数とします。即ち、

$$ \phi_1=1, \phi_2=x, \phi_3 = x^2, \cdots, \phi_n = x^{n-1} \:. $$

このとき、$p(x)$は以下の形を持っています。

$$ p(x) = \sum_{i=1}^{n} c_i x^{i-1} $$

節点を$[-1,1]$の$N-1$等分割点とする。即ち、

$$ x_i=-1 + 2(i-1)/(N-1), \quad i=1,2,\cdots, N \: . $$

最小二乗法による以下の関数の多項式近似を求めてください。 $$ f=\frac{1}{1+25x^2} \:. $$

特に、$N=30,60,90$, $n=10,20,30$を使ってください。

In [68]:
%与えられるc対して、多項式のグラフを描く関数
function draw_polynomial(c)
  x = -1:0.01:1;
  y = 0;
  n = (length(c)-1);
  for i=0:n 
    y = y + c(i+1) * x.^i;
  end
   plot(x,y,'b-')
   grid on
end
In [ ]:
N=30;
n=10;
xlist= -1:2/(N-1):1;

A=zeros(N,n);
b=zeros(N,1);

for i=1:N
  for j=1:n
     A(i,j)= ??
  end
  b(i) = myf( xlist(i) );
end

c=??

求めた$c$に対応する多項式を描きます。

In [ ]:
hold on
x=-1:0.01:1;
plot(x, myf(x),'r-')
draw_polynomial(c)
grid on