Julia1.0のDeepLearningフレームワークFlux.jl入門

#概要
Fluxの公式ドキュメントhttps://fluxml.ai/Flux.jl/stable/
に従ってFluxチュートリアルを行っていく。

#Fluxとは
FluxとはJuliaのDLフレームワーク

#インストール
`Julia1.0`かそれ以降で動作する。Pkgモード(`]`を入力する)で以下のように入力する。
`delete`キーを押すとjuliaモードに戻る。

```shell
(v1.0) pkg> add Flux
```

CUDAを使いたいときは以下も入力。(この記事では使わない。使い方はhttps://fluxml.ai/Flux.jl/stable/gpu.htmlを参考)

```shell
(v1.0) pkg> add CuArrays
```

#勾配を求める

##一変数関数f(x)

f(x)の導関数f´(x)を求める。 
(f´(x)の`´`はmacの場合`option`と`E`キーで入力できる。)


```julia
julia> using Flux.Tracker
julia> f(x) = 3x^2 + 2x + 1 
julia> f´(x) = Tracker.gradient(f, x)[1] #f´(x)=6x+2と定義
julia> f´(2) #f´(2)を求める
14.0 (tracked)
julia> f´´(x) = Tracker.gradient(f´, x)[1] # f´´(x)=6
julia> f´´(1)
6.0 (tracked)
```

##多変数関数
###勾配とは
n変数関数f(x1,x2,x3,・・・,xn)について
n次元ベクトル

```math
\nabla = \Biggl({\frac{\partial f}{\partial x_1}} ,..., {\frac{\partial f}{\partial x_n}}\Biggr)
```
を勾配と言う。n=1のときはただの導関数。

##勾配を求める

```math
f(W, b, x) = W \times x + b
```

この時の勾配は

```math
\nabla f(W, b, x) = \Bigl(x , 1 , W\Bigr)
```

W = 2, b = 3, x = 4のとき

```math
\nabla f(2, 3, 4) = \Bigl(4, 1, 2\Bigr)
```
これをjuliaで書くと以下のようになる。

```julia
julia> f(W, b, x) = W * x + b
julia> Tracker.gradient(f, 2, 3, 4)
(4.0 (tracked), 1.0, 2.0 (tracked))
```

##パラメータの数が多いとき
機械学習のモデルのパラメータが多いときは`param`というものを使って勾配を求める。
先ほどの例を`param`を使って書いてみる。

```julia
julia> W = param(2) #W=2を代入
julia> b = param(3) #W=3を代入
julia> f(x) = W * x + b
julia> params = Params([W, b])
julia> grads = Tracker.gradient(() -> f(4), params)
julia> grads[W] 
4.0
julia> grads[b]
1.0
```

この例ではパラメータが少ないので違いがあまり分からないが、
`param`を使うことによって、実際パラメータが多いとき、非常に扱いやすくなる。

#Simple Models
入力xからyを予測する簡単な線形回帰の問題を考えて見る。

2×5行列Wと2×1行列bを乱数で作成し,同様に入力として5成分ベクトルx,2成分ベクトルyも乱数で作成し、その時の損失関数

```math
loss(x, y) = \sum_{n=1}^{2}\Bigl(y_i-ŷ_i\Bigr)^2
```
を最小にする。

```julia
julia> W = rand(2, 5)  #Wを2×5行列を乱数で作成
 0.930514  0.396266  0.787115  0.682745  0.781531
 0.456005  0.84561   0.20131   0.703145  0.139113

julia> b = rand(2) #bを2×1行列を乱数で作成
 0.3862844664160474
 0.3715160899698007
julia> predict(x) = W*x .+ b
julia> function loss(x, y) #損失関数loss(x, y)を定義
         ŷ = predict(x)
         sum((y .- ŷ).^2)
       end
julia> x, y = rand(5), rand(2) # Dummy data
([0.98456, 0.6722, 0.629101, 0.497996, 0.513054], [0.180554, 0.920561])
julia> loss(x,y) 
7.9207555228077045
```

勾配を求めるため、損失関数lossを偏微分する。

```julia
julia> W = param(W)
 0.930514  0.396266  0.787115  0.682745  0.781531
 0.456005  0.84561   0.20131   0.703145  0.139113

julia> b = param(b)
 0.3862844664160474
 0.3715160899698007

julia> #損失関数lossを偏微分する
julia> gs = Tracker.gradient(() -> loss(x, y), Params([W, b]))
```

最急降下法で損失関数を最小にするWを得る。
update!(W, αΔ)は勾配gsの値によって以下のようにWの値を更新する。
αは更新する際の更新する度合いを定める定数である。

```math
Δ = gs\\
W = W - αΔ
```


```julia
julia> using Flux.Tracker: update!

julia> Δ = gs[W] #勾配を求める
 5.16775  3.52824  3.30202  2.61388  2.69291
 2.00165  1.36661  1.27898  1.01244  1.04306

julia> update!(W, -0.1Δ)
 0.413739  0.0434426  0.456913   0.421358  0.51224  
 0.25584   0.708949   0.0734118  0.601901  0.0348072

julia> loss(x, y)
2.2616927513479204 (tracked)
```

これをloss(x, y) < εとなるまで繰り返せばよい。
[参考]もし実際にデータを持っているなら実際にモデルを訓練してみよう
https://fluxml.ai/Flux.jl/stable/training/training.html

#Building Layers
実際は線形回帰よりもっと複雑なモデルを使うことが多い。
例えばsigmoid関数のような非線形活性化関数を使いニューラルネットを構成する。

```julia
julia> W1 = param(rand(3,5))
 0.17122    0.665763  0.0395654  0.533161  0.382145
 0.541056   0.781     0.212816   0.364648  0.335324
 0.0847882  0.737511  0.403253   0.684704  0.141392
julia> b1 = param(rand(3))
 0.8287821688331789
 0.5433551383280499
 0.2721052113521929
julia> layer1(x) = W1 * x .+ b1

julia> W2 = param(rand(2, 3))
 0.639634  0.702379  0.192444
 0.108319  0.146129  0.509777
julia> b2 = param(rand(2))
 0.31566080298798305
 0.820353365740117  
julia> layer2(x) = W2 * x .+ b2

julia> model(x) = layer2(σ.(layer1(x)))

julia> model(rand(5)) 
 1.593096701001898 
 1.4430508944882239
```
このままだと保守性がなく、変更修正、層の追加などをするときに、かなり扱いにくいので, 以下のように関数を使って書き直す。

```julia
julia> function linear(in, out)
         W = param(randn(out, in))
         b = param(randn(out))
         x -> W * x .+ b
       end

julia> linear1 = linear(5, 3) 
julia> linear2 = linear(3, 2)

julia> model(x) = linear2(σ.(linear1(x)))

julia> model(rand(5))
 -0.20356715657677094
  1.9012304757200134 
```

また以下のように構造体を使って実装することもできる。

```julia
julia> struct Affine
         W
         b
       end
julia> Affine(in::Integer, out::Integer) = 
         Affine(param(randn(out, in)), param(randn(out)))
julia> (m::Affine)(x) = m.W * x .+ m.b
julia> L1 = Affine(5, 3)
julia> L2 = Affine(3, 2)
julia> model2(x) = L2(σ.(L1(x)))
```

これらはFluxのDense layerの実装である。
Fluxは多くのレイヤーがあるが、全て簡単に実装することができる。

#Stacking It Up!

`Dense`は`Dense(input, output, 活性化関数)`のように使う。
3つ目の引数はなくてもよい。
上の例のL1は`Dense(5, 3, σ)`, L2は`Dense(3, 2)`と書き換えられる。


以下のようなモデルを書くことはとても多い。

```julia
layer1 = Dense(10, 5, σ)
# ...
model(x) = layer3(layer2(layer1(x)))
```

もっと直感的に書くと

```julia
julia> using Flux
julia> layers = [Dense(10, 5, σ), Dense(5, 2), softmax]     
julia> model(x) = foldl((x, m) -> m(x), layers, init = x)
julia> model(rand(10))# => 2-element vector
```

これはまた以下のようにもっと簡単に書ける。

```julia
julia> model2 = Chain(
         Dense(10, 5, σ),
         Dense(5, 2),
         softmax)

julia> model2(rand(10)) # => 2-element vector
```

もっとシンプルに書くこともできる。

```julia
julia> m = Dense(5, 2) ∘ Dense(10, 5, σ)
julia> m(rand(10)) # => 2-element vector
```

`Chain`はどんなjuliaの関数でも動くことができるところが優れている。

```julia
julia> m = Chain(x -> x^2, x -> x+1)
julia> m(5)
26
```