2022年度プログラミング演習A・B

第15回:最小二乗法

1. 最小二乗法とは

複数の点(データ)が与えられたときに、それらによく当てはまる関数を求めること(特に、パラメータを持つ関数について最適なパラメータを決定すること)を、統計学の用語で回帰(regression)と言います。与えられた点を通る必要はないため、前回と前々回の授業で扱った補間とは異なります。

$xy$ 平面上の $x$ 座標が相異なる $n$ 個の点 $(x_1,y_1),\ldots,(x_n,y_n)$ に対して、$m$ 個の関数 $g_1(x),\ldots,g_m(x)$ の線形結合として表される

$$ g(x)=\sum_{j=1}^m c_j g_j(x) $$

の最適なパラメータ $c_1,\ldots,c_m\in\mathbb{R}$ を決定することを考えます。例えば

$$ g_j(x)=x^{j-1}\quad (j=1\ldots,m) $$

の場合は $g(x)$ は $m-1$ 次の多項式、特に $m=2$ の場合は $g(x)$ は $1$ 次関数(直線)になり、実用上このように指定することが多いです。

(※最大次数の項の係数が $0$ になる可能性があるため、上記は厳密には「$m-1$ 次以下」です。)

あとは何をもって「パラメータが最適である」と言うかが問題となりますが、最小二乗法(least squares method)においては、残差平方和(residual sum of squares)

$$ \sum_{i=1}^n \left(g(x_i)-y_i\right)^2=\sum_{i=1}^n \left(\sum_{j=1}^m c_j g_j(x_i)-y_i\right)^2 $$

が最小となるようにパラメータ $c_1,\ldots,c_m\in\mathbb{R}$ を決定します。これは、大雑把に言えば「関数値とデータのずれを最小にしよう」という考え方です。

特に、ある関数 $f$ 上の点 $(x_1,f(x_1)),\ldots,(x_n,f(x_n))$ が与えられた場合には

$$ \sum_{i=1}^n \left(g(x_i)-f(x_i)\right)^2=\sum_{i=1}^n \left(\sum_{j=1}^m c_j g_j(x_i)-f(x_i)\right)^2 $$

が最小となるようにパラメータ $c_1,\ldots,c_m\in\mathbb{R}$ を決定することになり、このようにして得られる関数 $g(x)$ は $f(x)$ の近似と解釈することができます。

2. 理論

以下では、最小二乗法における最適なパラメータ $c_1,\ldots,c_m\in\mathbb{R}$ を理論的に求めます。

$n\times m$ 行列 $A$、$m$ 次元ベクトル $c$、$n$ 次元ベクトル $b$ を

$$ A= \begin{pmatrix} g_1(x_1) & g_2(x_1) & \cdots & g_m(x_1)\\ g_1(x_2) & g_2(x_2) & \cdots & g_m(x_2)\\ \vdots & \vdots & \vdots & \vdots\\ g_1(x_n) & g_2(x_n) & \cdots & g_m(x_n) \end{pmatrix} ,\quad c= \begin{pmatrix} c_1\\ c_2\\ \vdots\\ c_m \end{pmatrix} ,\quad b= \begin{pmatrix} y_1\\ y_2\\ \vdots\\ y_n \end{pmatrix} $$

とおくと、残差平方和について

$$ \sum_{i=1}^n \left(g(x_i)-y_i\right)^2=\sum_{i=1}^n \left(\sum_{j=1}^m c_j g_j(x_i)-y_i\right)^2=(Ac-b)^T(Ac-b) $$

が成り立ちます。そのため、関数 $r:\mathbb{R}^m\to\mathbb{R}$ を

$$ r(c)=(Ac-b)^T(Ac-b) $$

により定めれば、最小二乗法では $r(c)$ を最小にするようなベクトル $c$ を求めればよいことが分かります。

ここで、$r(c)$ の勾配ベクトルを計算すると

$$ \nabla r(c)= \begin{pmatrix} \frac{\partial r(c)}{\partial c_1}\\ \vdots\\ \frac{\partial r(c)}{\partial c_m} \end{pmatrix} =2A^T(Ac-b)=2A^TAc-2A^Tb $$

となります(証明は省略、$1$ 変数関数 $(ax-b)^2$ の微分が $2a(ax-b)$ であることを考えると理解しやすいです)。

$r(c)$ は凸関数であるため、$0_m$ を $m$ 次元の零ベクトルとして

$$ \nabla r(c)=0_m $$

を満たす $c$ を求めれば、それが $r(c)$ を最小にするベクトルになります。

この条件を変形すると連立一次方程式

$$ A^TAc=A^Tb $$

が得られ、これを正規方程式(normal equation)と呼びます。

正規方程式が唯一の解を持つかどうかは、行列 $A^TA$ の正則性に依存します。$A^TA$ を $A$ のグラム行列(Gram matrix)と呼び、次の二つは同値であることが知られています。

なお、$n\times m$ 行列 $A$ は $m$ 個の $n$ 次元ベクトルを列ベクトルとして持つため、もし $n<m$ であれば明らかに $A$ の列ベクトルは線形独立になりません。よって、最小二乗法を用いる上では、前提条件として少なくとも $n\geq m$ が必要であることに注意してください。

3. 使用例

3.1. 最小二乗法による線形回帰

$n$ 個の点 $(x_1,y_1),\ldots,(x_n,y_n)$ に対して、最適な $1$ 次関数

$$ g(x)=c_1+c_2x $$

を求めることを考えます。これを線形回帰(linear regression)と言い、求めた $1$ 次関数のグラフを回帰直線(regression line)と呼びます。

最小二乗法を用いるならば、前述の式で

$$ m=2,\quad g_1(x)=1,\quad g_2(x)=x $$

の場合を考えることになるため、

$$ A= \begin{pmatrix} 1 & x_1\\ 1 & x_2\\ \vdots & \vdots\\ 1 & x_n \end{pmatrix} ,\quad c= \begin{pmatrix} c_1\\ c_2 \end{pmatrix} ,\quad b= \begin{pmatrix} y_1\\ y_2\\ \vdots\\ y_n \end{pmatrix} $$

とおいて正規方程式

$$ A^TAc=A^Tb $$

を解けばよいことが分かります。

演習1

次の $11$ 個の点に対して、正規方程式を解いて最小二乗法による線形回帰を行い、それらの点および求めた回帰直線(範囲は $0\leq x\leq 1$)を描画してください。

$$ (0,1.09),\quad (0.1,1.17),\quad (0.2,1.42),\quad (0.3,1.54), $$$$ (0.4,1.71),\quad (0.5,1.99),\quad (0.6,2.26),\quad (0.7,2.44), $$$$ (0.8,2.55),\quad (0.9,2.84),\quad (1,3.04) $$

最小二乗法による線形回帰では、正規方程式を解かなくても $(x_1,y_1),\ldots,(x_n,y_n)$ に関する統計量を用いて

$$ c_2=\frac{s_{xy}}{s_x^2}, $$$$ c_1=\overline{y}-c_2\overline{x} $$

と計算することができます。

ここで、$\overline{x}$ は $x_1,\ldots,x_n$ の平均、$\overline{y}$ は $y_1,\ldots,y_n$ の平均、$s_x^2$ は $x_1,\ldots,x_n$ の分散、$s_{xy}$ は $x_1,\ldots,x_n$ と $y_1,\ldots,y_n$ の共分散です。

演習2

演習1と同じ $11$ 個の点に対して、統計量を用いた計算で最小二乗法による線形回帰を行い、それらの点および求めた回帰直線(範囲は $0\leq x\leq 1$)を描画してください。

ただし、統計量の計算には平均はmean(ベクトル)、分散はvar(ベクトル)、共分散はcov(ベクトル1,ベクトル2)を使用してよいです。

3.2. 最小二乗法による関数の多項式近似

関数 $f(x)$ と $n$ 個の節点 $x_1,\ldots,x_n$ が与えられたとき、最小二乗法によって $f(x)$ を近似する $m-1$ 次以下の多項式

$$ g(x)=\sum_{j=1}^m c_j x^{j-1} $$

を求めることができます(ただし、$n\geq m$)。前述のように、これは

$$ g_j(x)=x^{j-1}\quad (j=1\ldots,m) $$

の場合に相当します。

演習3(オプション)

節点

$$ x_i=-1+\frac{2(i-1)}{n-1}\quad (i=1,\ldots,n) $$

を用いて関数

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

を $m-1$ 次以下の多項式

$$ g(x)=\sum_{j=1}^m c_j x^{j-1} $$

で近似することを考えます。

$n=30,60$、$m=10,20$ を組み合わせた計 $4$ 個の場合について近似多項式を求め、元の関数 $f(x)$ とともに $-1\leq x\leq 1$ の範囲でグラフを描画してください。

第15回レポート課題

演習1~演習3に取り組んでください。