Ccmmutty logo
Commutty IT
40 min read

02 機械学習ライブラリの基礎

https://picsum.photos/seed/27a284ebf392410b99c9ecb898541c43/600/800
本章では,基礎的な機械学習手法として代表的な単回帰分析重回帰分析の仕組みを、数式を用いて説明します. ここで単回帰分析と重回帰分析を紹介することには 2 つの理由があります. 1 つ目は,回帰分析と重回帰分析の数学がニューラルネットワーク含めたディープラーニングの数学の基礎となるためです. 2 つ目は,単回帰分析のアルゴリズムを通して微分,重回帰分析のアルゴリズムを通して線形代数に関する理解を深めることができるためです.

単回帰分析

まずはじめに,単回帰分析について説明します. 機械学習手法は,教師あり学習 (supervised learning)教師なし学習 (unsupervised learning),**強化学習 (reinforcement learning)**に大別され,単回帰分析は教師あり学習に含まれます.
教師あり学習の中でも典型的な問題設定は 2 つに大別されます. 与えられた入力変数から,10100.10.1 といった実数値を予測する**回帰 (regression)と、「赤ワイン」,「白ワイン」といったカテゴリを予測する分類 (classification)**の 2 つです.
単回帰分析は回帰を行うための手法であり,1 つの入力変数から 1 つの出力変数を予測します. それに対し,重回帰分析は,複数の入力変数から 1 つの出力変数を予測します. この両手法は教師あり学習であるため,訓練の際には、入力変数 xx と目的変数 tt がペアで準備されている必要があります.

問題設定(単回帰分析)

まず,データに含まれる情報の中から何を利用し,何を予測させるかを決めます.
ここでは例として,家賃を予測する問題を考えることにします. 従って,家賃が 出力変数 yy となります.
次に, 入力変数 として何を採用するかを考えます. 家賃の予測には,部屋の広さ,駅からの距離,犯罪発生率などを考慮する必要があると思われます. ここでは部屋の広さを入力変数 xx として採用することにします. 複数の入力変数の候補がある場合に,それらを同時に扱う方法は,次の重回帰分析の説明の際に紹介します.
多くの機械学習手法は,大きく分けて次の3ステップで構成されています.
  • Step1: モデルを決める
  • Step2: 目的関数を決める
  • Step3: 最適なパラメータを求める
上記の3ステップについて,順に説明していきます.

Step1. モデルを決める(単回帰分析)

まずはじめに、入力変数 x と出力変数 y との関係をどのように定式化するかを決定します. この定式化したものをモデルもしくは数理モデルと呼びます.
どのように定式化すれば,家賃をうまく予測することができるのでしょうか. このモデル設計は現在は人手で行うのが一般的であり,機械が自動的に決めてくれるわけではありません(ただし最近ではAutoMLなどモデルを自動決定するための研究も進展してきています).
例えば,家賃と部屋の広さの組で表されるデータを 3 つ集め,「家賃」を y 軸に,「部屋の広さ」を x 軸にとってそれらをプロットしたとき,次のようになっていたとします.
家賃と部屋の広さの関係
この場合,部屋が広くなるほど,家賃が高くなるという関係が予想されます. また,この 2 変数間の関係性は直線によって表現を行うことができそうだと考えられます. そこで、2 つのパラメータ wwbb によって特徴づけられる直線の方程式
f(x;w,b)=wx+b f(x; w, b) = wx + b
によって,部屋の広さと家賃の関係を表すことを考えます. ここで,ww重み (weight)bb は**バイアス (bias)**の頭文字を採用しています.
直線式によるモデル化
単回帰分析では,このようにモデルとして直線 f(x;w,b)=wx+bf(x; w, b) = wx + b を用います. そして,2 つのパラメータ wwbb を,直線がデータによくフィットするように調整します.
パラメータで特徴づけられたモデルを用いる場合,与えられたデータセットに適合するように最適なパラメータを求めることが目標となります. 今回はデータセットとして部屋の広さ xx と家賃 tt の組からなるデータの集合を用います. 全部で NN 個のデータがあり,nn 番目のデータが (x(n),t(n))(x^{(n)}, t^{(n)}) と表されるとき,データセットは
D={(x(1),t(1)),(x(2),t(2)),,(x(N),t(N))}={(x(n),t(n)}n=1N \begin{aligned} \mathcal{D} &= \{(x^{(1)}, t^{(1)}), (x^{(2)}, t^{(2)}), \dots, (x^{(N)}, t^{(N)})\} \\ &= \{(x^{(n)}, t^{(n)}\}_{n=1}^N \end{aligned}
と表すことができます. これを用いて、新しい xx を入力すると,それに対応する tt を予測するモデルを訓練します.
ここで,データセット中の入力値とのことを**データ点(datum)**ということがあることに注意してください. データ点とは,具体的には上の説明で登場した D\mathcal{D} 中の各 (x(1),t(1))(x^{(1)},t^{(1)}) などのことです.
ここで,この後の計算を楽に進めるために,データの中心化というテクニックを紹介します. 部屋の広さと家賃は両方とも正の値であるため,各データ点をいくつかプロットすると,下図の左のグラフのようになります. 中心化では,各次元の平均が 0\boldsymbol{0} となるよう全てのデータを同量平行移動します. 中心化はしばしば前処理として採用されます. 厳密には前章で紹介したスケーリング方法の一つである標準化(正規化)がよく用いられます.
中心化処理
この処理を行うと,下図のように,バイアス bb00 とおけるため,fc(x;w)=wxcf_c(x; w) = wx_{c} のように,モデルをバイアス成分なしで表現することができるようになります. これによって,調整すべきパラメータを減らすことができます.
中心化後の直線式
データの中心化は入出力の平均をデータの全体から引くことで実現されます. つまり,
xc(n)=x(n)xˉtc(n)=t(n)tˉ \begin{aligned} x^{(n)}_{c} &= x^{(n)} - \bar{x} \\ t^{(n)}_{c} &= t^{(n)} - \bar{t} \end{aligned}
という変換を全ての nn について行います.
例えば,具体的な数値で見ると,下の表のようになります.
nx(n)x^{(n)}xˉ\bar{x}xc(n)x^{(n)}_c
113-2
2330
3532
中心化後を示す添え字の cc に関しては表現が冗長となるため,今後はこの添え字を省略し,データの中心化を事前に行っていることを前提とします. この時,モデルは
f(x;w)=wx f(x; w) = wx
となり,単回帰分析の目標は,データセット D={x(n),t(n)}n=1N\mathcal{D} = \{ x^{(n)}, t^{(n)} \}_{n=1}^{N} に基づいて,パラメータ ww適切に調整することになります.

Step2. 目的関数を決める(単回帰分析)

1章で説明したように,教師あり学習では多くの場合,目的関数を設計し,その目的関数を最小化(または最大化)することでモデルの訓練を行います.
今回は教師データと予測値が一致することが目標であり,乖離度を表す最小化すべき目的関数として教師データと予測値の二乗誤差を使います. ここで,モデルの出力が予測値となります.すなわち,y=f(x;w)y = f(x; w) です. 二乗誤差が 00 であれば t=yt = y となり予測が完全に教師データと一致していることを意味します. nn 番目の物件に対する教師データ t(n)t^{(n)} と予測値 y(n)=f(x(n);w)y^{(n)} = f(x^{(n)}; w) の二乗誤差は
(t(n)y(n))2 (t^{(n)} - y^{(n)})^2
です. 特定の物件についてだけ考慮するのではなく,データセット中の全ての物件の情報を考慮してモデルの訓練を行うために,上式で計算される各データ点における二乗誤差を全物件に対してそれぞれ計算して和をとったものを目的関数とします. すなわち,
L=(t(1)y(1))2+(t(2)y(2))2++(t(N)y(N))2=n=1N(t(n)y(n))2 \begin{aligned} \mathcal{L} &= \left( t^{(1)} - y^{(1)} \right)^2 + \left( t^{(2)} - y^{(2)} \right)^2 + \dots + \left( t^{(N)} - y^{(N)} \right)^2 \\ &= \sum^{N}_{n=1} \left( t^{(n)} - y^{(n)} \right)^2 \\ \end{aligned}
です. また,Step1で決めたモデル
y(n)=f(x(n);w)=wx(n) y^{(n)} = f(x^{(n)}; w) = wx^{(n)}
をこれに代入すると,目的関数は
L=n=1N(t(n)wx(n))2 \mathcal{L} = \sum^{N}_{n=1}\left( t^{(n)} - wx^{(n)} \right)^2
とパラメータを含んだ形式で表現することができます.このような関数を損失関数とよぶことを思い出してください.

Step3. 最適なパラメータを求める(単回帰分析)

次に目的関数を最小にするようなパラメータを求めます. 今回用いた目的関数は二次関数であるため,微分して求まる導関数を 0 とおいてそのときの x を求めれば,最小値を取る時の x が分かります. すなわち,目的関数の「傾きが0」となる点が値 00 をとる点です.
それでは,目的関数を微分しましょう.
wL=wn=1N(t(n)wx(n))2 \begin{aligned} \dfrac{\partial }{\partial w} \mathcal{L} &= \dfrac{\partial}{\partial w} { \sum^{N}_{n=1} ( t^{(n)}-wx^{(n)})^{2} } \end{aligned}
微分という操作は線形性を持っているため,和の微分は,微分の和と等しくなります. これを利用して次を得ます.
wL=n=1Nw(t(n)wx(n))2 \dfrac{\partial}{\partial w} \mathcal{L} = \sum^{N}_{n=1} \dfrac {\partial }{\partial w} \left( t^{(n)}-wx^{(n)} \right)^2
パラメータ ww による微分を表す w\dfrac {\partial }{\partial w} と総和を表す \sum が入れ替わっています. 次に和の各項を見ます.
w(t(n)wx(n))2 \dfrac {\partial }{\partial w}\left( t^{(n)}-wx^{(n)} \right)^2
t(n)wx(n)t^{(n)} - wx^{(n)}()2(\cdot)^2合成関数になっています.
u(n)=t(n)wx(n)u^{(n)} = t^{(n)} - wx^{(n)}, g(u(n);w)=(u(n))2g(u^{(n)}; w) = (u^{(n)})^2 とおくと,
w(t(n)wx(n))2=wg(u(n))=u(n)wg(u(n))u(n)=x(n)(2u(n))=2x(n)(t(n)wx(n)) \begin{aligned} \dfrac {\partial }{\partial w} \left( t^{(n)} - wx^{(n)} \right)^2 &= \dfrac {\partial }{\partial w} g(u^{(n)}) \\ &= \dfrac {\partial u^{(n)}}{\partial w} \dfrac{\partial g(u^{(n)})}{\partial u^{(n)}} \\ &= -x^{(n)} \cdot (2 u^{(n)}) \\ &= -2x^{(n)}( t^{(n)} - wx^{(n)} ) \end{aligned}
が得られます.これより,
wL=n=1Nw(t(n)wx(n))2=n=1N2x(n)(t(n)wx(n)) \begin{aligned} \dfrac{\partial }{\partial w} \mathcal{L} &= \sum^{N}_{n=1} \dfrac {\partial }{\partial w} \left( t^{(n)} - wx^{(n)} \right)^2 \\ &= -\sum^{N}_{n=1} 2x^{(n)} \left( t^{(n)} - wx^{(n)} \right) \end{aligned}
となります.この微分の値を 0 とおいて ww について解くと,
wL=02n=1Nx(n)(t(n)wx(n))=02n=1Nx(n)t(n)+2n=1Nw(x(n))2=02n=1Nx(n)t(n)+2wn=1N(x(n))2=0wn=1N(x(n))2=n=1Nx(n)t(n) \begin{aligned} \dfrac {\partial }{\partial w} \mathcal{L} &= 0 \\ -2 \sum^{N}_{n=1} x^{(n)} \left( t^{(n)} - wx^{(n)} \right) &= 0 \\ -2 \sum^{N}_{n=1} x^{(n)} t^{(n)} + 2 \sum^{N}_{n=1} w(x^{(n)})^2 &= 0 \\ -2 \sum^{N}_{n=1} x^{(n)} t^{(n)} + 2w \sum^{N}_{n=1} (x^{(n)})^2 &= 0 \\ w \sum^{N}_{n=1} (x^{(n)})^2 &= \sum^{N}_{n=1} x^{(n)} t^{(n)} \end{aligned}
より,
w=n=1Nx(n)t(n)n=1N(x(n))2 \begin{aligned} w &= \dfrac {\displaystyle \sum^{N}_{n=1} x^{(n)} t^{(n)}} {\displaystyle \sum^{N}_{n=1} (x^{(n)})^2} \end{aligned}
と求まります. この右辺に着目すると,与えられたデータセット D={x(n),t(n)}n=1N\mathcal{D} = \{x^{(n)}, t^{(n)}\}_{n=1}^{N} のみから決定されていることがわかります.
それでは,以下の数値例を使って実際にパラメータ ww を求めてみましょう.
nx(n)x^{(n)}t(n)t^{(n)}
112
223.9
336.1
まずは,データの中心化を行うために,平均を求めます.
xˉ=13(1+2+3)=2tˉ=13(2+3.9+6.1)=4 \begin{aligned} \bar{x} &= \dfrac{1}{3} (1 + 2 + 3) = 2 \\ \bar{t} &= \dfrac{1}{3}(2 + 3.9 + 6.1) = 4 \end{aligned}
各データ点からそれぞれの値を引きます.
x1=12=1x2=22=0x3=32=1t1=24=2t2=3.94=0.1t3=6.14=2.1 \begin{aligned} x_{1} &= 1 - 2 = -1 \\ x_{2} &= 2 -2 = 0 \\ x_{3} &= 3- 2 = 1\\ t_{1} &= 2 - 4 = -2\\ t_{2} &= 3.9 - 4 = -0.1\\ t_{3} &= 6.1 - 4 = 2.1 \end{aligned}
それでは,中心化後の値を用いて,最適なパラメータww を計算します.
w=n=1Nx(n)t(n)n=1N(x(n))2=x(1)t(1)+x(2)t(2)+x(3)t(3)(x(1))2+(x(2))2+(x(3))2=1×(2)+0×0.1+1×2.1(1)2+02+12=2.05 \begin{aligned} w &= \dfrac {\displaystyle \sum_{n=1}^{N} x^{(n)} t^{(n)}} {\displaystyle \sum_{n=1}^{N} (x^{(n)})^2} \\ &= \dfrac { x^{(1)} t^{(1)} + x^{(2)} t^{(2)} + x^{(3)} t^{(3)} } { (x^{(1)})^2 + (x^{(2)})^2 + (x^{(3)})^2 } \\ &= \dfrac { -1 \times (-2) + 0 \times 0.1 + 1 \times 2.1 }{ (-1)^2 + 0^2 + 1^2 } \\ &= 2.05 \end{aligned}
これで単回帰分析の学習が完了しました. この求まったパラメータを w^\hat{w} とすると,これを使用したモデル f(x;w^)f(x; \hat{w})学習済みモデルとなります.
続いて,このモデルを使って新しいサンプルに対する予測をしてみましょう. 学習したモデルを使って新たな入力データについて予測値を計算する処理を**推論 (inference)**とよびます. 例えば,新しいサンプル x(q)=1.5x^{(q)}=1.5 に対する予測値は次のように求まります,
y(q)tˉ=w^(x(q)xˉ)y(q)=w^(x(q)xˉ)+tˉ=2.05×(1.52)+4=2.975 \begin{aligned} y^{(q)} - \bar{t} &= \hat{w}(x^{(q)}-\bar{x}) \\ \Rightarrow y^{(q)} &= \hat{w}(x^{(q)}-\bar{x}) + \bar{t} \\ &= 2.05 \times (1.5 - 2) + 4 \\ &= 2.975 \end{aligned}
モデルは中心化データを用いて学習を行ったので,推論の際にも入力値・予測値それぞれに同様の操作を行う必要があることに注意しましょう.
以上が,単回帰分析の一連の手順となります.

重回帰分析

次に,多変数の入力変数を扱う重回帰分析を扱います. この重回帰分析を学ぶことで線形代数に関する知識が深まります.
重回帰分析は単回帰分析と同様に教師あり学習の一種であり,回帰を行う手法です. 問題設定はほとんど単回帰分析と同じですが,重回帰分析では入力変数の数が複数となります. つまり,複数の入力変数から出力変数を予測する機械学習手法の一つです.

問題設定(重回帰分析)

ここでは単回帰分析の場合と同様に家賃を予測する問題を考え,家賃を出力変数 yy とします. 入力変数としては,単回帰分析では考慮しきれていなかった駅からの距離や犯罪発生率なども同時に考慮します. 例えば,部屋の広さ x1x_{1}, 駅からの距離 x2x_{2}, ..., 犯罪発生率 xMx_{M} のように MM 個の入力変数があるとします(M=1M=1の場合,単回帰分析の問題に帰着されます).
単回帰分析と同様,以下の3ステップに従います.
  • Step1: モデルを決める
  • Step2: 目的関数を決める
  • Step3: 最適なパラメータを求める

Step1. モデルを決める(重回帰分析)

単回帰分析のモデルは,
f(x;w,b)=wx+b f(x; w, b) = wx + b
であり,ww を重み(weight),bb をバイアス (bias) とよびました. 重回帰分析では,この式を複数の入力変数へと拡張し,
f(x;w,b)=w1x1+w2x2++wMxM+b f({\bf x}; {\bf w}, b) = w_{1}x_{1} + w_{2}x_{2} + \dots + w_{M}x_{M} + b
のような線形結合の形で表します. ここで,太字の x,w{\bf x}, {\bf w} はそれぞれベクトルを表し,それぞれの mm 番目の要素が wm,xmw_m, x_m で表されています.
ここでは,各入力変数は出力変数に線形に影響を与えることが仮定されています. 一方,もし入力変数間に非線形な依存関係が想定される場合には,そのことを考慮したモデル化を行う必要があります. それについては後述します.
重回帰分析のモデルは総和の記号を使って整理すると,
f(x;w,b)=m=1Mwmxm+b f({\bf x}; {\bf w}, b) = \sum_{m=1}^{M} w_{m} x_{m} + b
のように書くことができます.さらにここで,x0=1x_0 = 1w0=bw_0 = bとおくと,
f(x;w,b)=w1x1+w2x2++wMxM+b=w1x1+w2x2++wMxM+w0x0=w0x0+w1x1++wMxM=m=0Mwmxm \begin{aligned} f({\bf x}; {\bf w}, b) &= w_{1}x_{1} + w_{2}x_{2} + \dots + w_{M}x_{M} + b\\ &= w_{1}x_{1} + w_{2}x_{2} + \dots + w_{M}x_{M} + w_{0}x_{0}\\ &= w_{0}x_{0} + w_{1}x_{1} + \dots + w_{M}x_{M} \\ &= \sum_{m=0}^M w_m x_m \end{aligned}
のようにバイアス bb を総和記号の中に含めることができます. w0w_0 を含む w{\bf w} を改めて w{\bf w},最初の要素 x0x_0 として 11 を付け加えた入力ベクトルを改めて x{\bf x} と定義しなおすと,このモデルはベクトルを用いて
f(x;w)=[w0w1wM][x0x1xM]=wTx=xTw \begin{aligned} f({\bf x}; {\bf w}) &= \begin{bmatrix} w_{0} & w_{1} & \dots & w_{M} \end{bmatrix} \begin{bmatrix} x_{0} \\ x_{1} \\ \vdots \\ x_{M} \end{bmatrix} \\ &= {\bf w}^{T} {\bf x} = {\bf x}^T {\bf w} \end{aligned}
と書けます. これで,ベクトルの内積を用いて表現することができました.
これで重回帰分析を行うためのモデルが決定できました. このモデルはパラメータとして M+1M+1 個の重み w{\bf w} を持っています.

Step2. 目的関数を決める(重回帰分析)

単回帰分析では,目標値 tt と予測値 yy の二乗誤差が小さいほど良い予測であるとし,その総和を目的関数にしました. 重回帰分析でも,そのような予測値 yy を求めるというのは同じであるため,同じ目的関数を使います.
ここで,nn 個目の入力値 x(n){\bf x}^{(n)} に対する予測値 f(x(n);w)f({\bf x}^{(n)}; {\bf w})y(n)y^{(n)},目標値を t(n)t^{(n)} とおくと,目的関数は
L=(t(1)y(1))2+(t(2)y(2))2++(t(N)y(N))2 \begin{aligned} \mathcal{L} &= \left( t^{(1)} - y^{(1)} \right)^2 + \left( t^{(2)} - y^{(2)} \right)^2 + \dots + \left( t^{(N)} - y^{(N)} \right)^2 \end{aligned}
となります. 単回帰分析では,これを総和記号 \sum を用いて表すことができましたが,これは以下のようにベクトル同士の内積を用いて表すこともできます.
L=(t(1)y(1))2+(t(2)y(2))2++(t(N)y(N))2=[t(1)y(1)t(2)y(2)t(N)y(N)][t(1)y(1)t(2)y(2)t(N)y(N)]=(ty)T(ty) \begin{aligned} \mathcal{L} &= \left( t^{(1)} - y^{(1)} \right)^2 + \left( t^{(2)} - y^{(2)} \right)^2 + \dots + \left( t^{(N)} - y^{(N)} \right)^2 \\ &= \begin{bmatrix} t^{(1)} - y^{(1)} & t^{(2)} - y^{(2)} & \dots & t^{(N)} - y^{(N)} \end{bmatrix} \begin{bmatrix} t^{(1)} - y^{(1)} \\ t^{(2)} - y^{(2)} \\ \vdots \\ t^{(N)} - y^{(N)} \end{bmatrix} \\ &= \left( {\bf t} - {\bf y} \right)^{\rm T} \left( {\bf t} - {\bf y} \right) \end{aligned}
また,y{\bf y}
y=[y(1)y(2)y(N)]=[(x(1))Tw(x(2))Tw(x(N))Tw]=[(x(1))T(x(2))T(x(N))T]w \begin{aligned} {\bf y} = \begin{bmatrix} y^{(1)} \\ y^{(2)} \\ \vdots \\ y^{(N)} \end{bmatrix} = \begin{bmatrix} ({\bf x}^{(1)})^{\rm T} {\bf w} \\ ({\bf x}^{(2)})^{\rm T} {\bf w} \\ \vdots \\ ({\bf x}^{(N)})^{\rm T} {\bf w} \end{bmatrix} = \begin{bmatrix} ({\bf x}^{(1)})^{\rm T} \\ ({\bf x}^{(2)})^{\rm T} \\ \vdots \\ ({\bf x}^{(N)})^{\rm T} \end{bmatrix} \boldsymbol{w} \end{aligned}
のように書くことができます. ここで,x(n){\bf x}^{(n)}mm 番目の要素を xnmx_{nm} と書くことにすると,上式はさらに
y=[x10x11x12x1Mx20x21x22x2MxN0xN1xN2xNM][w0w1w2wM]=Xw \begin{aligned} {\bf y} &= \begin{bmatrix} x_{10} & x_{11} & x_{12} & \dots & x_{1M} \\ x_{20} & x_{21} & x_{22} & \dots & x_{2M} \\ \vdots & \vdots & \vdots & \ddots \\ x_{N0} & x_{N1} & x_{N{2}} & \dots & x_{NM} \end{bmatrix} \begin{bmatrix} w_{0} \\ w_{1} \\ w_{2} \\ \vdots \\ w_{M} \end{bmatrix} \\ &= {\bf X}{\bf w} \end{aligned}
と表せます. X{\bf X} の各行が各サンプルを表しており,各列が入力変数を表しています. つまり,各サンプルごとに MM 個の入力変数(xn0x_{n0} は常に 11)を持ちます.
具体例を挙げてみます. 例えば,m=1m=1 が部屋の広さ,m=2m=2 が駅からの距離,m=3m=3 が犯罪発生率に対応する入力変数だとして,nn 個目のデータ点が部屋の広さ =50m2= 50m^{2} ,駅からの距離 =600m= 600 m ,犯罪発生率 =2= 2% であるような物件を表している場合,これは
xnT=[1506000.02] {\bf x}_{n}^{T} = \begin{bmatrix} 1 & 50 & 600 & 0.02 \end{bmatrix}
というベクトルになります. 今,3 つの入力変数を考えているので,M=3M=3 ですが,入力ベクトル xn{\bf x}_n の1つ目の要素 xn1x_{n1} はバイアス w0w_0 に対応する要素で常に値は 11 となることに注意してください.

Step3. パラメータを最適化する(重回帰分析)

それでは,目的関数を最小化するパラメータ w{\bf w}を求めてみましょう.
※ここでは,式変形を駆使しながら最適パラメータの解析的な解を求めていきますが,導出過程が少々複雑になります.導出結果は次節(§2.3)で示されているので,導出に興味のある方以外は本節はスキップしていただいて構いません.
まずは目的関数をパラメータ w{\bf w} を用いて表し直すと
L=(ty)T(ty)=(tXw)T(tXw)={tT(Xw)T}(tXw)=(tTwTXT)(tXw) \begin{aligned} \mathcal{L} &= \left( {\bf t} - {\bf y} \right)^{\rm T} \left( {\bf t} - {\bf y} \right) \\ &= \left( {\bf t} - {\bf X}{\bf w} \right)^{\rm T} \left( {\bf t} - {\bf X}{\bf w} \right) \\ &= \left\{ {\bf t}^{\rm T} - ({\bf X}{\bf w})^{\rm T} \right\} \left( {\bf t} - {\bf X}{\bf w} \right) \\ &= \left( {\bf t}^{\rm T} - {\bf w}^{\rm T}{\bf X}^{\rm T} \right) \left( {\bf t} - {\bf X}{\bf w} \right) \end{aligned}
となります. ここで,転置の公式 (AB)T=BTAT({\bf A}{\bf B})^{\rm T} = {\bf B}^{\rm T}{\bf A}^{\rm T} を使っていることに注意しましょう. さらに分配法則を使って展開すると,
L=tTttTXwwTXTt+wTXTXw \begin{aligned} \mathcal{L} &= {\bf t}^{\rm T}{\bf t} - {\bf t}^{\rm T}{\bf X}{\bf w} - {\bf w}^{\rm T}{\bf X}^{\rm T}{\bf t} + {\bf w}^{\rm T}{\bf X}^{\rm T}{\bf X}{\bf w} \\ \end{aligned}
となります. この目的関数に対しパラメータの ww についての偏微分を計算します. その前にこの式はもう少し整理することができます. はじめに,
(1)T=1 (1)^T = 1
のようにスカラは転置しても変化しません. 上式の中で出てくる tTXw{\bf t}^{\rm T}{\bf X}{\bf w} はスカラなので,
(tTXw)T=tTXw ({\bf t}^{\rm T}{\bf X}{\bf w})^{\rm T} = {\bf t}^{\rm T}{\bf X}{\bf w}
が成り立ちます. さらに,転置の公式 (ABC)T=CTBTAT({\bf A}{\bf B}{\bf C})^{\rm T} = {\bf C}^{\rm T}{\bf B}^{\rm T}{\bf A}^{\rm T} より,
(tTXw)T=wTXTt ({\bf t}^{\rm T}{\bf X}{\bf w})^{\rm T} = {\bf w}^{\rm T}{\bf X}^{\rm T}{\bf t}
も成り立ちます.これより,
(tTXw)T=tTXw=wTXTt ({\bf t}^{\rm T}{\bf X}{\bf w})^{\rm T} = {\bf t}^{\rm T}{\bf X}{\bf w} = {\bf w}^{\rm T}{\bf X}^{\rm T}{\bf t}
を導くことができます. 目的関数を L\mathcal{L} とおくと,上の式を利用して,
L=tTt2tTXw+wTXTXw \begin{aligned} \mathcal{L} = {\bf t}^{\rm T}{\bf t} - 2 {\bf t}^{\rm T}{\bf X}{\bf w} + {\bf w}^{\rm T}{\bf X}^{\rm T}{\bf X}{\bf w} \\ \end{aligned}
とまとめることができます.ここで, w{\bf w} について偏微分を行いやすくするため, w{\bf w} 以外の定数項をまとめると,
L=tTt2tTXw+wTXTXw=tTt2(XTt)Tw+wTXTXw=γ+βTw+wTAw \begin{aligned} \mathcal{L} &= {\bf t}^{\rm T}{\bf t} - 2 {\bf t}^{\rm T}{\bf X}{\bf w} + {\bf w}^{\rm T}{\bf X}^{\rm T}{\bf X}{\bf w} \\ &= {\bf t}^{\rm T}{\bf t} - 2 \left( {\bf X}^{\rm T}{\bf t} \right)^{\rm T}{\bf w} + {\bf w}^{\rm T}{\bf X}^{\rm T}{\bf X}{\bf w} \\ &= \gamma + \boldsymbol{\beta}^{\rm T}{\bf w} + {\bf w}^{\rm T}{\bf A}{\bf w} \end{aligned}
となります. 線形代数の章で学んだような w{\bf w} に関する二次形式(二次関数)で表現することができました. ここで,A=XTX, β=2XTt, γ=tTt{\bf A}= {\bf X}^{\rm T}{\bf X}, \ \boldsymbol{\beta} =-2 {\bf X}^{\rm T}{\bf t}, \ \gamma = {\bf t}^{\rm T}{\bf t} です. また,β\boldsymbol{\beta} を転置の形式にした理由は,線形代数の章で学んだベクトルで微分するための公式集にある形式に合わせるためです.
それでは,この目的関数を最小化するパラメータ w{\bf w} の求め方を考えましょう. 前述の通り,目的関数はパラメータ w{\bf w} に関しては二次関数になっています. 例えば,
w=[w1w2],A=[1234],β=[12],γ=1 \begin{aligned} {\bf w} = \begin{bmatrix} w_{1} \\ w_{2} \end{bmatrix}, {\bf A} = \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix}, \boldsymbol{\beta} = \begin{bmatrix} 1 \\ 2 \end{bmatrix}, \gamma = 1 \end{aligned}
のとき,
L=wTAw+βTw+γ=[w1w2][1234][w1w2]+[12][w1w2]+1=[w1w2][w1+2w23w1+4w2]+w1+2w2+1=w1(w1+2w2)+w2(3w1+4w2)+w1+2w2+1=w12+5w1w2+4w22+w1+2w2+1 \begin{aligned} \mathcal{L} & = {\bf w}^{\rm T}{\bf A}{\bf w} + \boldsymbol{\beta}^{\rm T}{\bf w} + \gamma \\ &= \begin{bmatrix} w_{1} & w_{2} \end{bmatrix} \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix} \begin{bmatrix} w_{1} \\ w_{2} \end{bmatrix} + \begin{bmatrix} 1 & 2 \end{bmatrix} \begin{bmatrix} w_{1} \\ w_{2} \end{bmatrix} + 1 \\ &= \begin{bmatrix} w_{1} & w_{2} \end{bmatrix} \begin{bmatrix} w_{1} + 2w_{2} \\ 3w_{1} + 4w_{2} \end{bmatrix} + w_{1} + 2 w_{2} + 1 \\ &= w_{1} \left( w_{1} + 2w_{2} \right) + w_{2} \left( 3w_{1} + 4w_{2} \right) + w_{1} + 2w_{2} + 1 \\ &= w^{2}_{1} + 5w_{1}w_{2} + 4w^{2}_{2} + w_{1} + 2w_{2} + 1 \\ \end{aligned}
となります. さらに w1,w2w_{1}, w_{2} に関してそれぞれまとめると,
L=w12+(5w2+1)w1+(4w22+2w2+1)=4w22+(5w1+2)w2+(w12+w1+1) \begin{aligned} \mathcal{L} &= w^{2}_{1} + \left( 5w_{2} + 1\right) w_{1} + \left( 4w^{2}_{2} + 2w_{2} + 1 \right) \\ &= 4w^{2}_{2} + \left( 5w_{1} + 2 \right) w_{2} + \left( w^{2}_{1} + w_{1} + 1 \right) \end{aligned}
となり,それぞれの二次関数であることがわかります.
そして,二次関数であれば,下図のような形となります.
パラメータと目的関数の関係(2次元)
これを3次元でイメージすると,下図のようになります.
パラメータと目的関数の関係(3次元)
そして,目的関数である二乗誤差の総和が最小となる点では各変数で微分した時の傾きが0となります.
目的関数が最小となる点
この例では,w1w_{1}w2w_{2} の2つのパラメータの場合で考えましたが,これは w0w_{0}, w1w_{1}, w2w_{2}, \ldots, wMw_{M} の場合でも同様に考えることができ,目的関数が最小となる点は
{w0L=0w1L=0     wML=0 \begin{cases} \dfrac {\partial }{\partial w_{0}}\mathcal{L} = 0\\ \dfrac {\partial }{\partial w_{1}}\mathcal{L} = 0\\ \ \ \ \ \ \vdots \\ \dfrac {\partial }{\partial w_{M}}\mathcal{L} = 0\\ \end{cases}
となり,これをまとめると,
[w0Lw1LwML]=[000]wL=0 \begin{aligned} \begin{bmatrix} \dfrac {\partial}{\partial w_{0}} \mathcal{L} \\ \dfrac {\partial}{\partial w_{1}} \mathcal{L} \\ \vdots \\ \dfrac {\partial}{\partial w_{M}} \mathcal{L} \\ \end{bmatrix}&=\begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ \end{bmatrix} \\ \Rightarrow \dfrac {\partial}{\partial {\bf w}} \mathcal{L} &= \boldsymbol{0} \\ \end{aligned}
のようにベクトルでの微分として表されます. あとは,上式を満たす w{\bf w} を決めます. まずは w{\bf w} が求めやすくなるように,代入と式変形を行います. (下記の計算ではベクトルでの微分をはじめとして線形代数の章で学んだ内容を活用しているため,計算途中がわからなくなった場合には,線形代数の章を再度確認しながら進めてください.)
wL=w(γ+βTw+wTAw)=0w(γ)+w(βTw)+w(wTAw)=00+β+(A+AT)w=02XTt+{XTX+(XTX)T}w=02XTt+2XTXw=0XTXw=XTt \begin{aligned} \dfrac {\partial }{\partial {\bf w}} \mathcal{L} = \dfrac {\partial }{\partial {\bf w}} \left( \gamma + \boldsymbol{\beta}^{\rm T}{\bf w} + {\bf w}^{\rm T}{\bf A}{\bf w} \right) = {\bf 0} \\ \dfrac {\partial }{\partial {\bf w}} \left( \gamma \right) + \dfrac {\partial }{\partial {\bf w}} \left( \boldsymbol{\beta}^{\rm T}{\bf w} \right) + \dfrac {\partial }{\partial {\bf w}} \left( {\bf w}^{\rm T}{\bf A}{\bf w} \right) = {\bf 0} \\ {\bf 0} + \boldsymbol{\beta} + \left( {\bf A}+{\bf A}^{\rm T} \right) {\bf w} = {\bf 0} \\ -2{\bf X}^{\rm T}{\bf t} + \left\{ {\bf X}^{\rm T}{\bf X} + \left( {\bf X}^{\rm T}{\bf X} \right)^{\rm T} \right\} {\bf w} = {\bf 0} \\ -2{\bf X}^{\rm T}{\bf t}+2{\bf X}^{\rm T}{\bf X}{\bf w}={\bf 0}\\ {\bf X}^{\rm T}{\bf X}{\bf w}={\bf X}^{\rm T}{\bf t}\\ \end{aligned}
ここで,XTX{\bf X}^{\rm T} {\bf X}に逆行列が存在すると仮定して,両辺に左側から (XTX)1\left( {\bf X}^{\rm T}{\bf X}\right)^{-1} をかけると,
(XTX)1XTXw=(XTX)1XTtIw=(XTX)1XTtw=(XTX)1XTt \begin{aligned} \left( {\bf X}^{\rm T}{\bf X}\right)^{-1}{\bf X}^{\rm T}{\bf X} {\bf w} =\left( {\bf X}^{\rm T}{\bf X}\right)^{-1}{\bf X}^{\rm T}{\bf t} \\ {\bf I}{\bf w}=\left( {\bf X}^{\rm T}{\bf X}\right)^{-1}{\bf X}^{\rm T}{\bf t} \\ {\bf w}=\left( {\bf X}^{\rm T}{\bf X}\right)^{-1}{\bf X}^{\rm T}{\bf t} \end{aligned}
となり,与えられたデータセット X,t{\bf X}, {\bf t} から,最適なパラメータ w{\bf w} が求まりました.ここで,I{\bf I} は単位行列です.また,式変形の際には,
w=XTtXTX {\bf w} = \dfrac{{\bf X}^{\rm T}{\bf t}}{{\bf X}^{\rm T}{\bf X}}
のような分数が表れないように注意してください.これは行列の計算には割り算が定義されていないためです. そのため,逆行列を使って行列積のみで計算しています.
また,もうひとつよくある間違いとして,w{\bf w} を求めるために以下のような式変形をしてしまう例が挙げられます.
XTXw=XTt(XT)1XTXw=(XT)1XTtXw=tX1Xw=X1tw=X1t \begin{aligned} {\bf X}^{\rm T}{\bf X}{\bf w} &= {\bf X}^{\rm T}{\bf t} \\ \left( {\bf X}^{\rm T} \right)^{-1} {\bf X}^{\rm T}{\bf X}{\bf w} &= \left( {\bf X}^{\rm T} \right)^{-1}{\bf X}^{\rm T}{\bf t} \\ {\bf X}{\bf w} &= {\bf t} \\ {\bf X}^{-1}{\bf X}{\bf w} &= {\bf X}^{-1}{\bf t} \\ {\bf w} &= {\bf X}^{-1}{\bf t} \end{aligned}
しかし,これは一般には成立しません. その理由は,逆行列を持つための条件の一つである正方行列であることX{\bf X} が常に満たすとは限らないためです. サンプル数 NN と入力変数の数 M+1M+1 が等しくない場合,XRN×(M+1){\bf X} \in \mathcal{R}^{N \times (M+1)} は正方行列ではなく,逆行列をもちません. 一方,XTX{\bf X}^{\rm T} {\bf X}XTXR(M+1)×(M+1){\bf X}^{\rm T}{\bf X} \in \mathcal{R}^{(M+1) \times (M+1)} であり,サンプル数 NN に依存することなく,常に正方行列となることに注目してください.(逆行列が求まるためにはもう少し厳密な条件がありますが,ここでは説明しません.)
推論の際は学習で得られたパラメータを w^\hat{\bf w} として
yq=w^Txq y_{q} = \hat{\bf w}^{\rm T}{\bf x}_{q}
のように計算することで予測値が得られます.
で表されることが分かりました.この最適なパラメータを計算するために,以下の5つを扱います.
  • ベクトルの定義
  • 行列の定義
  • 転置
  • 行列積
  • 逆行列
具体的に,以下のようなデータセットが与えられているケースを想定してみましょう.この例では,データのサンプル数NN44であり,入力データXXの変数の数は22です.そしてttは教師データとなります.
X=[123125134159], t=[1568] \boldsymbol{X} = \begin{bmatrix} 1 & 2 & 3 \\ 1 & 2 & 5 \\ 1 & 3 & 4 \\ 1 & 5 & 9 \end{bmatrix}, \ \boldsymbol{t} = \begin{bmatrix} 1 \\ 5 \\ 6 \\ 8 \end{bmatrix}
ここでX\boldsymbol{X}パラメータ w\boldsymbol{w} がバイアス b\boldsymbol{b} を包含する 形式を想定しており,従って入力データX\boldsymbol{X}の1列目には11が格納されています.
それでは実装方法について見ていきましょう.まずは,NumPyの読み込みから始めます.

NumPyによる実装

それでは重回帰分析で行われる計算をPythonを用いてコンピュータに実行させてみましょう. PythonにはNumPyとよばれる線形代数を簡単に扱えるライブラリが存在し,広く利用されています. 次の章で紹介するChainerの中でもNumPyは用いられており,様々なデータ解析・科学計算のライブラリで広く採用されているため,NumPyの使い方にある程度慣れ親しんでおくことは後々役立ちます.
以下では,Pythonの基本的な文法はすでに理解していることを前提としています. 例えば,変数(数値・文字列,リスト,タプル,辞書),制御構文(for,if),関数,クラスを理解している必要があります.
重回帰分析では,最終的に最適なパラメータ w{\bf w}
w=(XTX)1XTt {\bf w} = \left( {\bf X}^{\rm T}{\bf X}\right)^{-1}{\bf X}^{\rm T}{\bf t}
と閉じた形で求まりました. この計算を行うために,これから以下の5つをNumPyを使って行います.
  • ベクトルの定義
  • 行列の定義
  • 転置
  • 行列積
  • 逆行列
以下のようなデータセットが与えられているとします.
X=[123125134159], t=[1568] {\bf X} = \left[ \begin{matrix} 1 & 2 & 3 \\ 1 & 2 & 5 \\ 1 & 3 & 4 \\ 1 & 5 & 9 \end{matrix} \right], \ {\bf t} = \left[ \begin{matrix} 1 \\ 5 \\ 6 \\ 8 \end{matrix} \right]
データ数 N=4N = 4 であり,入力変数の数 M=2M = 2 です. それではこのデータセットにフィットするモデル f(x;w)=xTwf({\bf x}; {\bf w}) = {\bf x}^{\rm T}{\bf w} のパラメータ w{\bf w} をNumPyを使って計算してみましょう.
まずNumPyを読み込みます.
python
import numpy as np
ベクトルの定義は以下のように行います.
python
t = np.array([1, 5, 6, 8])
ベクトルを表示してみましょう.
python
print(t)
[1 5 6 8]
行列の定義も行い,表示してみましょう.
python
X = np.array([
    [1, 2, 3],
    [1, 2, 5],
    [1, 3, 4],
    [1, 5, 9]
])
python
print(X)
[[1 2 3] [1 2 5] [1 3 4] [1 5 9]]
ここではnp.arrayという関数を用いて,PythonのリストからNumPyの多次元配列の形式(np.ndarray)への変換を行っています.
次に,Xの転置を行ってみましょう.np.ndarrayで定義されている場合,.Tをつけるだけで転置することができます.
python
print(X.T)
[[1 1 1 1] [2 2 3 5] [3 5 4 9]]
縦と横が入れ替わっていることを確認できます.
行列積は以下のように np.dot によって実現できます.行列積を行う際には,一番目の行列の列数と,二番目の行列の行数が同じであることに注意して下さい.
python
XX = np.dot(X.T, X)
python
print(XX)
[[ 4 12 21] [ 12 42 73] [ 21 73 131]]
ここからさらに,XTX{\bf X}^{\rm T}{\bf X} に対する逆行列,(XTX)1\left( {\bf X}^{\rm T}{\bf X} \right)^{-1} を計算します.逆行列を求めるには,np.linalg.inv を用います.
python
XX_inv = np.linalg.inv(XX)
python
print(XX_inv)
[[ 1.76530612 -0.39795918 -0.06122449] [-0.39795918 0.84693878 -0.40816327] [-0.06122449 -0.40816327 0.24489796]]
これで重回帰分析のために必要な演算が揃いました.
最適なパラメータ (XTX)1XTt\left({\bf X}^{\rm T}{\bf X} \right)^{-1} {\bf X}^{\rm T}{\bf t} を求めると,
python
Xt = np.dot(X.T, t)
python
print(Xt)
[ 20 70 124]
python
w = np.dot(XX_inv, Xt)
python
print(w)
[-0.14285714 0.71428571 0.57142857]
このように求まりました. NumPyを使って数式と同じ計算をプログラムに書き下せば,コンピュータに具体的な数値を使った計算を高速に行わせることができます.

Scikit-learnによる機械学習アルゴリズムの実行

重回帰分析であればNumPyで比較的容易に実装することができましたが,実践的に使用する機械学習手法のアルゴリズムの多くは複雑であり,初学者が一から書くのは難しい場合も少なくありません. PythonにはScikit-learnと呼ばれる様々な機械学習手法の実装が含められたライブラリが公開されており,初学者でも簡単に扱うことができます.
ここでは重回帰分析をScikit-learnを用いて行う方法を紹介します. データセットは先程と同様に X{\bf X}t{\bf t} を使用しますが,Scikit-learnにおいては,パラメータ w{\bf w} がバイアス b{\bf b} を包含しない 形式を想定しており,入力データ X{\bf X} の1列目から 11 を取り除く必要があることに注意してください. 従って,
X=[23253459], t=[1568] {\bf X} = \begin{bmatrix} 2 & 3 \\ 2 & 5 \\ 3 & 4 \\ 5 & 9 \end{bmatrix}, \ {\bf t} = \begin{bmatrix} 1 \\ 5 \\ 6 \\ 8 \end{bmatrix}
をデータセットとして準備します.

Scikit-learn 基礎編

Scikit-learnはsklearnという名前で読み込みます.
python
import sklearn
重回帰分析を使用する場合は以下の LinearRegression というクラスを読み込みます.
python
from sklearn.linear_model import LinearRegression
なお,使い方を調べる際には,公式のリファレンスに加えて,使用例を見るのも有用です(例えば検索エンジンで「重回帰分析 Scikit-learn」のようなキーワードで検索すればたくさんの使用例が見つかります).
LinearRegression クラスを利用するために,インスタンス化を行い,model と名付けます.
python
model = LinearRegression()
これで重回帰分析を行う準備が完了しました. この model を使って,最適なパラメータを求めるには以下のようにします.
python
# データセットの定義
X = np.array([
    [2, 3],
    [2, 5],
    [3, 4],
    [5, 9]
])
t = np.array([1, 5, 6, 8])

# 最適なパラメータの計算
model.fit(X, t)
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False)
求まったパラメータの検証結果を以下のようにして確認することができます.
python
model.score(X, t)
0.6923076923076923
LinearRegression クラスのオブジェクトが持つ score() メソッドを呼ぶと,以下の式で示される,決定係数とよばれる指標が計算され,結果が返されます.
R2=1i(tiyi)2i(titˉ)2 R^{2} = 1 - \dfrac{\sum_{i}\left( t_{i} - y_{i} \right)^{2}}{\sum_{i}\left( t_{i} - \bar{t} \right)^{2}}
Scikit-learnは,よく使う機械学習手法を簡単に切り替えて使えるよう,統一されたインターフェースで様々な機械学習手法が実装されています. Scikit-learnが好まれる理由の一つには,様々なアルゴリズムが「.fit()で学習,.score()で検証」という同じインターフェースから利用できることが挙げられるでしょう.
また,アルゴリズムによって内容は多少異なりますが,パラメータも model オブジェクトに属性として格納されているため,学習後に確認することができます.
python
# パラメータw
model.coef_
array([0.71428571, 0.57142857])
python
# バイアスb
model.intercept_
-0.14285714285714235

Scikit-learn 応用編

Scikit-learnは機械学習の実装を支援する多くの機能を兼ね備えています.本節では,サンプルデータセットの使用方法,及びデータセットの分割方法について紹介していきます.
サンプルデータセットの使用
まずはじめにサンプルデータセットの取り扱いを紹介します. Scikit-learnでは,幾つかのデータセットが提供されています. 今回はその中から,米国ボストン市郊外における地域別の物件価格のデータセットを使用することにします.
このデータセットには 506 件のデータが登録されており,各サンプルには対象地域の平均物件価格と,それに紐づく情報として対象地域の平均的な物件情報(一戸あたりの部屋数,築年数,雇用施設からの距離など),人口統計情報(低所得者の割合,教師あたりの生徒数など),生活環境に関する情報(犯罪発生率など)などが含まれています. このデータセットを利用して,物件や人口統計などの情報から,平均物件価格を予測するモデルを構築してみましょう.
入力変数は全部で13種類あり,詳細は以下の通りです.
  • CRIM : 人口11人あたりの犯罪発生率
  • ZN : 25,00025,000平方フィート以上の住宅区画が占める割合
  • INDUS : 非小売業が占める面積の割合
  • CHAS : チャールズ川に関するダミー変数 (1 : 川沿い,0 : それ以外)
  • NOX : 窒素酸化物の濃度
  • RM : 住居あたりの平均部屋数
  • AGE : 1940年以前に建てられた物件の割合
  • DIS : 5つのボストン雇用施設からの重み付き距離
  • RAD : 都心部の幹線道路へのアクセス指数
  • TAX : $ 10,00010,000あたりの固定資産税の割合
  • PTRATIO : 教師1人あたりの生徒数
  • B : 黒人の比率を表す指数
  • LSTAT : 低所得者の割合
それでは, load_boston() 関数を実行して,データセットを読み込みましょう.
python
from sklearn.datasets import load_boston
python
boston = load_boston()
boston はPythonの辞書と同じインターフェースが備わっており,'data' キーに入力値が,'target' キーに目標値が格納されています.
python
X = boston['data']
t = boston['target']
python
print(X)

Discussion

コメントにはログインが必要です。