OSX
Xcode
機械学習
Swift
Accelerate

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

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

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

勾配降下法とは?

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

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

y = x^2 + n

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

SS2.jpg

説明の前にコード

この式の 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で 新規プロジェクトでmacOSCommand Line Toolを選択して、

SS1.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

説明

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

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

先ほどの例は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.swiftmain.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

偏微分すると
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 )

標準正規分布をマイナスしたものです。

SS6.jpg

これの導関数は

{ 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]
見つけました

次回は

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

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