Julia の紹介
機械学習のための Julia

2015/11/28 機械学習 名古屋 第2回勉強会
後藤 俊介 ( @antimon2 )

自己紹介

  • 名前:後藤 俊介
  • 所属:株式会社コスモルート クラウドR&Dグループ 数学班
  • 言語:Ruby, JavaScript, Python, Julia, Haxe, …
  • twitter: @antimon2
  • Facebook: antimon2
  • GitHub: antimon2

Julia とは?

  • The Julia Language
  • 2015/10/04 に v0.4.0 がリリース(2015/11/27 現在の最新は v0.4.1)
  • Python/Ruby/R 等の「いいとこどり」言語(詳細後述)
  • 動作が速い!(LLVM JIT コンパイル)

Julia の特長

  • Rのように中身がぐちゃぐちゃでなく、
  • Rubyのように遅くなく、
  • Lispのように原始的またはエレファントでなく、
  • Prologのように変態的なところはなく、
  • Javaのように硬すぎることはなく、
  • Haskellのように抽象的すぎない

ほどよい言語である

引用:http://www.slideshare.net/Nikoriks/julia-28059489

Julia の目指すもの:

  • C のように高速だけど、
    Ruby のような動的型付言語である
  • Lisp のように同じ文法で書けるマクロがあって、しかも
    Matlab のような直感的な数式表現もできる
  • Python のように総合的なプログラミングができて、
    R のように統計処理も得意で、
    Perl のように文字列処理もできて、
    Matlab のように線形代数もできて、
    shell のように複数のプログラムを組み合わせることもできる
  • 超初心者にも習得は簡単で、
    超上級者の満足にも応えられる
  • インタラクティブにも動作して、コンパイルもできる

Why We Created Julia から抜粋・私訳)

機械学習への適用

自分で実装編。

例1:単純な線形回帰

1変数(or 多くても10変数程度)、データ数100件程度
⇒ 行列作ってバックスラッシュ演算子で高速に求解可能。

In [1]:
# 手作りサンプルデータ
x = rand(100) .* 100
y = randn(100) .* 5 .+ x .* 0.7 .+ 10

# 行列作成
X = [ones(100, 1) x;]

# パラメータ求解(Normal Equation 利用)
# θ = X \ y
θ = (X' * X) \ (X' * y);
In [2]:
# 結果プロット
using Gadfly  # ← Julia 定番のデータプロットライブラリ

x_linspace = linspace(0,100)
f(x) = θ[1] .+ x .* θ[2]
plot(layer(x=x, y=y, Geom.point), layer(x=x_linspace, y=f(x_linspace), Geom.line))
Out[2]:
x -150 -100 -50 0 50 100 150 200 250 -100 -95 -90 -85 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 165 170 175 180 185 190 195 200 -100 0 100 200 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 -100 -80 -60 -40 -20 0 20 40 60 80 100 120 140 160 180 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 -100 0 100 200 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 y

例2:簡単な K-means

本日のハンズオンで利用したサンプルデータを利用。

In [3]:
# データ読込
data = open(readdlm, "./CodeIQ_data.txt", "r")

K=3
X=data;
In [4]:
# クラスタ重心初期化
function initCentroids(X, K)
    randidx = randperm(size(X, 1))
    X[randidx[1:K], :]
end
Out[4]:
initCentroids (generic function with 1 method)
In [6]:
# クラスタ割り付けステップ
function findClosestCentroids(X, centroids)
    K = size(centroids, 1)
    map(1:size(X, 1)) do i
        x = vec(X[i, :])
        ts = [x - vec(centroids[j,:]) for j=1:K]
        ls = [dot(t, t) for t=ts]
        indmin(ls)
    end
end
Out[6]:
findClosestCentroids (generic function with 1 method)
In [8]:
# 重心移動ステップ
function computeCentroids(X, idxs, K)
    m, n = size(X)  # m 不使用
    centroids = zeros(K, n)
    for j = 1:K
        centroids[j, :] = mean(X[idxs .== j, :], 1) # ← idxs が j(=1,2,…,K)の行のみ抽出して、列ごとに平均値を算出
    end
    centroids
end
Out[8]:
computeCentroids (generic function with 1 method)
In [10]:
# ループ実施
function runKMeans(X, initialCentroids, max_iters=10)
    m, n = size(X)
    K = size(initialCentroids, 1)
    centroids = initialCentroids
    prev_centroids = centroids
    idxs = zeros(m, 1)
    
    for i = 1:max_iters
        idxs = findClosestCentroids(X, centroids)
        centroids = computeCentroids(X, idxs, K)
        if centroids == prev_centroids
            break
        end
        prev_centroids = centroids
    end
    
    return centroids, idxs
end
Out[10]:
runKMeans (generic function with 2 methods)
In [11]:
# 結果表示
function dispResult(X, idxs)
    plotdata = [X [string(i) for i=idxs];]
    plot(x=plotdata[:,1], y=plotdata[:,2], color=plotdata[:,3])
end
Out[11]:
dispResult (generic function with 1 method)
In [12]:
# 実行
initial_centroids = initCentroids(X, K)
centroids, idxs = runKMeans(X, initial_centroids)
dispResult(X, idxs)
Out[12]:
x -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 -20 0 20 40 -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 3 2 1 Color -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 -20 0 20 40 -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 y

機械学習への適用 その2

追加パッケージ利用編。

機械学習ライブラリ(主なもの)

  • MLBase(機械学習の基本アルゴリズムを集めたパッケージ。他↓のいくつかのライブラリで required されている)
  • Regression.jl(○○回帰(線形回帰、リッジ回帰、ロジスティック回帰等))
  • LIBSVM(SVMパッケージ。libsvm の Julia ラッパー)
  • Clustering.jl(K-means その他クラスタリング)
  • DecisionTree.jl(決定木、ランダムフォレスト)
  • BackpropNeuralNet(ニューラルネット)

例1:Regression.jl で 線形回帰

先ほど利用した手作りサンプルデータを利用

In [ ]:
# パッケージ読込
# Pkg.add("Regression")
using Regression
# ※ 2015/11/27 現在、大量の Warning が出力されます(>_<) ので出力隠しておきます
In [14]:
# 行列作成
X = [ones(100, 1) x;];
In [15]:
# 例1-1: linearreg + Regression.solve で求解
ret = Regression.solve(linearreg(X', y); options=Regression.Options(verbosity=:iter))
 Iter
Out[15]:
RiskMinSolution:
- sol:       (2,) Array{Float64,1}
- fval:      1148.2793332367219
- niters:    8
- converged: true
        f.value       f.change         g.norm           step
==================================================================
    0     1.2168e+05
    1     1.9749e+04    -1.0193e+05     1.0584e+05     1.9073e-06
    2     1.6630e+04    -3.1193e+03     9.9224e+04     6.2500e-02
    3     2.1258e+03    -1.4504e+04     2.4939e+04     5.0000e-01
    4     1.3927e+03    -7.3318e+02     1.2470e+04     5.0000e-01
    5     1.2094e+03    -1.8329e+02     6.2348e+03     5.0000e-01
    6     1.1636e+03    -4.5823e+01     3.1174e+03     5.0000e-01
    7     1.1521e+03    -1.1456e+01     1.5587e+03     5.0000e-01
    8     1.1483e+03    -3.8186e+00     2.3874e-12     1.0000e+00
Converged with 8 iterations @ f.value = 1.1483e+03
In [16]:
# 結果プロット
θ1 = ret.sol

x_linspace = linspace(0,100)
f(x) = θ1[1] .+ x .* θ1[2]
plot(layer(x=x, y=y, Geom.point), layer(x=x_linspace, y=f(x_linspace), Geom.line))
Out[16]:
x -150 -100 -50 0 50 100 150 200 250 -100 -95 -90 -85 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 165 170 175 180 185 190 195 200 -100 0 100 200 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 -100 -80 -60 -40 -20 0 20 40 60 80 100 120 140 160 180 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 -100 0 100 200 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 y
In [17]:
# 例1-2: llsq で求解
θ2 = llsq(X, y)

x_linspace = linspace(0,100)
f(x) = θ2[1] .+ x .* θ2[2]
plot(layer(x=x, y=y, Geom.point), layer(x=x_linspace, y=f(x_linspace), Geom.line))
Out[17]:
x -150 -100 -50 0 50 100 150 200 250 -100 -95 -90 -85 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 165 170 175 180 185 190 195 200 -100 0 100 200 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 -100 -80 -60 -40 -20 0 20 40 60 80 100 120 140 160 180 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 -100 0 100 200 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 y

例2:Clustering.jl で K-means

本日のハンズオンで利用したサンプルデータを利用。

In [18]:
# パッケージ読込
# Pkg.add("Clustering")
using Clustering
In [19]:
# 学習(data は先ほど読み込んだデータ)
X = data
K = 3

model = kmeans(X', K; maxiter=10, display=:iter)
  
Out[19]:
Clustering.KmeansResult{Float64}(2x3 Array{Float64,2}:
  3.91875  14.21    10.1009 
 11.1456   13.9303   2.90441,[3,1,2,3,3,1,3,3,1,2  …  3,1,3,3,2,1,1,3,2,1],[8.57266,6.92678,7.7371,13.2633,0.0355967,5.11937,5.06906,3.04261,2.91053,8.28019  …  0.765279,1.4695,0.88455,9.90113,7.78251,2.4756,1.9285,3.75704,14.6554,15.5151],[32,30,34],[32.0,30.0,34.0],545.1195459313719,2,true)
Iters               objv        objv-change | affected 
-------------------------------------------------------------
      0       6.472333e+02
      1       5.451195e+02      -1.021138e+02 |        0
      2       5.451195e+02       0.000000e+00 |        0
K-means converged with 2 iterations (objv = 545.1195459313719)
In [20]:
# 結果確認(dispResult 関数は先ほど定義したもの(Clustering パッケージには含まれません))
dispResult(X, model.assignments)
Out[20]:
x -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 -20 0 20 40 -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 3 1 2 Color -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 -20 0 20 40 -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 y