勾配降下法（Gradient Descent）を Swift で動かしながら理解する最小のコード

すみません、最小は煽りかもしれません・・・汗）

去年の正月休みに、[Coursera の Stanford 大学の Machine Learning](https://www.coursera.org/learn/machine-learning) コースを[取ってみました](https://www.coursera.org/account/accomplishments/verify/EDZYAHHNGCZA)。Andrew Ng 先生がとても丁寧に教えてくれるほんとにいいコースだと思います。終わった時は泣きそうでした。おすすめです。復習をかねて教わった勾配降下法を Swift で書きながら理解しなおしてみたいと思います。( XCode Version 10.0 Beta 4 でチェックしています )

# 勾配降下法とは？

ネットにとてもたくさん説明が載っていますが　[このページがとてもイメージをつかみやすい](http://yaju3d.hatenablog.jp/entry/2017/08/27/233459) ので参照してみてください。
http://yaju3d.hatenablog.jp/entry/2017/08/27/233459


ボウルにピンポン球を入れたら最終的に底にとまるように、下に凸な関数に投げ入れたピンポン球がとまったところの`x`を求めるイメージです。

# `y = x^2 + n`

例として`y = x^2 - 2`をグラフ化してみます。

![SS2.jpg](https://qiita-image-store.s3.amazonaws.com/0/73493/6aa7a141-4fa4-349c-9c53-09b7a54190b2.jpeg)
## 説明の前にコード


この式の `y` が最小値の時の `x` を求めるコードの実行部分は最低限だとこんな感じにすればよさそうです。説明は実行結果の後に書きます。

```GradientDescent.swift
func
勾配降下法(
	導関数		: ( Double ) -> Double
,	初期値		: Double
,	学習率		: Double
,	コールバック	: ( Double ) -> Bool
) {
	var	w = 初期値
	repeat { w = w - 学習率 * 導関数( w ) } while コールバック( w )
}
```

```main.swift
var	i = 0
勾配降下法(
	導関数: { p in 2.0 * p }	//	x^2 + n の導関数
,	初期値: 1.0
,	学習率: 0.1
 ) { p in					//	コールバック
	print( i, p )
	i += 1
	return i < 40			//	とりあえず 40 回で終わり
}
```

プログラムの名前とかパラメータを説明のため日本語にしています。普段からそうしてるわけではありません・・・汗）

`XCode`で 新規プロジェクトで`macOS`の`Command Line Tool`を選択して、

![SS1.png](https://qiita-image-store.s3.amazonaws.com/0/73493/70a9b679-758d-b2f5-3984-b9fb6ece70c5.png)

ソースを上の2つのファイルにして実行してみてください。


実行結果は以下のようになります。徐々に0に近づいているのがわかります。

```
0 [0.8]
1 [0.64]
2 [0.512]
3 [0.4096]
4 [0.32768]
5 [0.26214400000000004]
6 [0.20971520000000005]
7 [0.16777216000000003]
8 [0.13421772800000004]
9 [0.10737418240000003]

途中略

30 [0.0009903520314283045]
31 [0.0007922816251426436]
32 [0.000633825300114115]
33 [0.0005070602400912919]
34 [0.0004056481920730336]
35 [0.00032451855365842687]
36 [0.0002596148429267415]
37 [0.0002076918743413932]
38 [0.00016615349947311455]
39 [0.00013292279957849163]
```

可視化するアプリを作ってみたのでそれで見てみます。

![ezgif-5-8d97b2bf5a.gif](https://qiita-image-store.s3.amazonaws.com/0/73493/93dae432-16ee-01b8-1bf3-9062815e7c33.gif)

## 説明

数学の記号を極力用いないでざっくりイメージで書きます。用語が数学の規定どおりでないのは承知の上です。責めないでください。

```GradientDescent.swift
func
勾配降下法(
	導関数		: ( Double ) -> Double
,	初期値		: Double
,	学習率		: Double
,	コールバック	: ( Double ) -> Bool
) {
	var	w = 初期値
	repeat { w = w - 学習率 * 導関数( w ) } while コールバック( w )
}
```

* `初期値`とは関数の中に投げ入れられたピンポン球が着地した位置のイメージです。ただしバウンドしません。
* `学習率`とはピンポン球のスピードのイメージです。ただし底に行くほど進み方は遅くなります。

* `x^2`を微分すると`2x`になります。

	Swiftで書けば：
`func F( p: Double ) -> Double { return p * p }`
を微分すると
`func F( p: Double ) -> Double { return 2 * p }`
になります。この関数を`導関数`と呼びます。

	ラムダで書けば：
`{ p in p * p }`
を微分すると 
`{ p in 2 * p }`
になります。

* 導関数の値は接線の傾きを表します。例えば` x^2 + n `の関数の` x = 1 `での接線の傾きは` 2 * x `、すなわち` 2 `です。

* 接線の傾きの符号の逆方向にピンポン球は進みます。

* コールバックに現在のピンポン球の位置を渡し、コールバックが` false `を返すまで計算し続けます。


# 多変数化

## `z = x^2 + y^2 + n`

`z = x^2 + y^2 - 10`のグラフ

![ezgif-5-fd6d0a335a.gif](https://qiita-image-store.s3.amazonaws.com/0/73493/8bdd9b78-278f-19ba-48b7-257327ee3d1b.gif)

先ほどの例は1変数でしたがこの式は変数が複数あります。なので、変数を`Array`で渡します。それに応じて他のものも変更します。また`Array`同士の引き算と掛け算（数学でいうベクトルの積とは違います）とスカラー値と`Array`の掛け算を定義します。

## まず演算から。

以下のファイルを作成してください。`Accelerate`フレームワークを使っているので、アップルオンリーです。

```MiniJPMatrix.swift
import Accelerate

func
- ( _ l: [ Double ], _ r: [ Double ] ) -> [ Double ] {
	guard l.count == r.count else { fatalError() }
	var	v = [ Double ]( repeating: 0, count: l.count )
	vDSP_vsubD( r, 1, l, 1, &v, 1, vDSP_Length( v.count ) )
	return v
}

func
* ( _ l: [ Double ], _ r: [ Double ] ) -> [ Double ] {
	guard l.count == r.count else { fatalError() }
	var	v = [ Double ]( repeating: 0, count: l.count )
	vDSP_vmulD( l, 1, r, 1, &v, 1, vDSP_Length( v.count ) )
	return v
}

func
* ( _ s: Double, _ p: [ Double ] ) -> [ Double ] {
	var	v = [ Double ]( repeating: 0, count: p.count )
	var	w = s
	vDSP_vsmulD( p, 1, &w, &v, 1, vDSP_Length( v.count ) )
	return v
}
```

これにより以下のようなことができるようになります。

`[ 4.0, 3.0 ] - [ 1.0, 2.0 ] -> [ 3.0, 1.0 ]`
`[ 1.0, 2.0 ] * [ 3.0, 4.0 ] -> [ 3.0, 8.0 ]`
`         2.0 * [ 3.0, 4.0 ] -> [ 6.0, 8.0 ]`

## パラメータを`Array`にする

`GradientDescent.swift`と`main.swift`を変更します。

```GradientDescent.swift
func
勾配降下法(
	導関数		: ( [ Double ] ) -> [ Double ]
,	初期値		: [ Double ]
,	学習率		: [ Double ]
,	コールバック	: ( [ Double ] ) -> Bool
) {
	var	w = 初期値
	repeat { w = w - 学習率 * 導関数( w ) } while コールバック( w )
}
```

```main.swift
var	i = 0
勾配降下法(
	導関数: { p in [ 2.0, 2.0 ] * p }	//	x^2 + y^2 + n の導関数
,	初期値: [ 1.0, 1.0 ]
,	学習率: [ 0.1, 0.1 ]
 ) { p in					//	コールバック
	print( i, p )
	i += 1
	return i < 40			//	とりあえず 40 回で終わり
}
```

Andrew 先生のコースでは学習率はスカラー値でしたが、各変数ごとに学習率を変えれるようにこれも`Array`にしてあります。

## 実行結果

```
0 [0.8, 0.8]
1 [0.64, 0.64]
2 [0.512, 0.512]
3 [0.4096, 0.4096]
4 [0.32768, 0.32768]
5 [0.26214400000000004, 0.26214400000000004]
6 [0.20971520000000005, 0.20971520000000005]
7 [0.16777216000000003, 0.16777216000000003]
8 [0.13421772800000004, 0.13421772800000004]
9 [0.10737418240000003, 0.10737418240000003]

途中略

30 [0.0009903520314283045, 0.0009903520314283045]
31 [0.0007922816251426436, 0.0007922816251426436]
32 [0.000633825300114115, 0.000633825300114115]
33 [0.0005070602400912919, 0.0005070602400912919]
34 [0.0004056481920730336, 0.0004056481920730336]
35 [0.00032451855365842687, 0.00032451855365842687]
36 [0.0002596148429267415, 0.0002596148429267415]
37 [0.0002076918743413932, 0.0002076918743413932]
38 [0.00016615349947311455, 0.00016615349947311455]
39 [0.00013292279957849163, 0.00013292279957849163]
```

## 多変数の導関数

多変数の式の導関数を求めるには、次のように偏微分といわれることをします。

例えば
`x^2 + y^2 + n`
を x に関して微分すると
`2 * x`
を y に関して微分すると
`2 * y`
となります。

上のソースでは、`導関数`に`{ p in [ 2.0, 2.0 ] * p }`を渡しています。

この`p`には
`[ x, y ]`
が渡ってくるので、
`[ 2.0, 2.0 ] * [ x, y ]`
すなわち
`[ 2.0 * x, 2.0 * y ]`
を返します。

# さまざまな関数とその導関数

下に凸な関数の導関数さえあれば極値の`x`を求められるのでやってみます。

## `z = 3x^2 + 5y^2 + 6xy`

![ezgif-5-9c703a0cf8.gif](https://qiita-image-store.s3.amazonaws.com/0/73493/ea105ad1-010d-64d8-a824-9bbd14c844e8.gif)

偏微分すると
`x: 6x + 6y`
`y: 10y + 6x`

なので`導関数`パラメータに以下の関数を渡します。

```
{ p in
	let	x = p[ 0 ]
	let	y = p[ 1 ]
	return [
		6.0 * x + 6 * y
	,	10.0 * y + 6 * x
	]
}
```
### 結果

```
0 [-0.20000000000000018, -0.6000000000000001]
1 [0.28, 0.1200000000000001]
2 [0.039999999999999925, -0.16800000000000004]
3 [0.1168, -0.023999999999999938]
4 [0.06111999999999996, -0.07008]
5 [0.06649599999999999, -0.03667199999999997]
6 [0.04860159999999998, -0.039897599999999984]
7 [0.04337919999999998, -0.029160959999999993]
8 [0.03484825599999999, -0.026027519999999988]
9 [0.02955581439999999, -0.02090895359999999]

途中略

30 [0.000624413447362416, -0.0004500517582240578]
31 [0.000519796433879401, -0.0003746480684174496]
32 [0.00043270741460223015, -0.0003118778603276406]
33 [0.0003602096820374764, -0.00025962444876133806]
34 [0.0002998585420717934, -0.00021612580922248588]
35 [0.0002496189023622089, -0.00017991512524307603]
36 [0.0002077966360907292, -0.00014977134141732533]
37 [0.0001729814592866869, -0.0001246779816544375]
38 [0.00014399937270733726, -0.00010378887557201211]
39 [0.00011987307442614217, -8.639962362440236e-05]
```

## `y = -1 / sqrt( 2 * .pi ) * exp( - x^2 / 2 )`

[標準正規分布](https://ja.wikipedia.org/wiki/正規分布)をマイナスしたものです。

![SS6.jpg](https://qiita-image-store.s3.amazonaws.com/0/73493/a1e75f86-4809-8286-68a3-71a7035987c9.jpeg)

これの導関数は

```
{ p in
	return [ p[ 0 ] / sqrt( 2 * .pi ) * exp( -( p[ 0 ] * p[ 0 ] ) / 2 ) ]
}
```

なので

```main.swift
import	Foundation

var	i = 0
勾配降下法(
	導関数: { p in
		[ p[ 0 ] / sqrt( 2 * .pi ) * exp( -( p[ 0 ] * p[ 0 ] ) / 2 ) ]
	}
,	初期値: [ 1.0 ]
,	学習率: [ 0.2 ]
 ) { p in					//	コールバック
	print( i, p )
	i += 1
	return i < 40			//	とりあえず 40 回で終わり
}
```

### 実行結果

```
0 [0.9516058550961713]
1 [0.9033268092100339]
2 [0.8553984061164297]
3 [0.8080592550971047]
4 [0.7615440585887202]
5 [0.716076598552253]
6 [0.6718632012597742]
7 [0.6290871573562863]
8 [0.5879044666066391]
9 [0.5484411247642597]

途中略

30 [0.10415888579607616]
31 [0.09589316862989072]
32 [0.08827709816616808]
33 [0.08126099582784446]
34 [0.07479867816390755]
35 [0.06884727899439798]
36 [0.06336706426678124]
37 [0.05832124467255749]
38 [0.05367578977518854]
39 [0.04939924637363279]
```

# 終了条件

今までの例ではとりあえず40回回して終了、と乱暴にやってきましたが、実際にはどうするかの一例を紹介します。

まず`MiniJPMatrix.swift`に以下の関数を付け加えます。

```
func
ユークリッド距離の2乗( _ p: [ Double ] ) -> Double {
	var	v = 0.0
	vDSP_svesqD( p, 1, &v, vDSP_Length( p.count ) )
	return v
}
```
ユークリッド距離は直角3角形の斜辺を求める`sqrt( x^2 + y^2 )`ってやつです。それの2乗なので、x^2 + y^2 を返します。
ちなみにユークリッド距離は3次元だと`sqrt( x^2 + y^2 + z^2 )`となります。

`main.swift`を以下のようにします。

```main.swift
import	Foundation

var	反復数マックス = 1000

var	しきい値の2乗 = 1E-8

var	i = 0
var	r = [ 1.0 ]


勾配降下法(
	導関数: { p in
		[ p[ 0 ] / sqrt( 2 * .pi ) * exp( -( p[ 0 ] * p[ 0 ] ) / 2 ) ]
	}
,	初期値: r
,	学習率: [ 0.2 ]
 ) { p in					//	コールバック

	if i >= 反復数マックス {
		print( "規定回数に達しました" )
		return false
	} else if ユークリッド距離の2乗( p - r ) < しきい値の2乗 {
		print( "見つけました" )
		return false
	} else {
		for w in p {
			if w == .infinity || w == -.infinity {
				print( "発散しました" )
				return false
			}
		}
	}

	print( i, p )
	i += 1
	r = p
	return true
}
```

前回コールバックが呼ばれた時の値を保存しておいて、新たな値とのユークリッド距離がしきい値以下なら見つかったと判断します。
上の例ではユークリッド距離の2乗がしきい値（1E-4)の2乗(1E-8)以下なら見つかったと判断しています。

またどこかの値が`.infinity`に達してしまったら無限に発散していると判断します。

## 実行結果

```
0 [0.9516058550961713]
1 [0.9033268092100339]
2 [0.8553984061164297]
3 [0.8080592550971047]
4 [0.7615440585887202]
5 [0.716076598552253]
6 [0.6718632012597742]
7 [0.6290871573562863]
8 [0.5879044666066391]
9 [0.5484411247642597]

途中略

80 [0.0016346121656212512]
81 [0.0015041891588787229]
82 [0.001384172364013263]
83 [0.001273731493938356]
84 [0.0011721025070172348]
見つけました
```

# 次回は

次のヒマな時があれば、次のことをやりたいと思っています。

* この記事に加筆、特に終了条件、学習率など
* 視覚化の記事を書く
* 機械学習への応用の記事を書く