『機械学習のエッセンス(http://isbn.sbcr.jp/93965/)』のPythonサンプルをJuliaで書き換えてみる。（第04章01数値計算の基本）

はじめに

『機械学習のエッセンス　実装しながら学ぶPython、数学、アルゴリズム』を本屋で見つけて買ってみました。この本は機械学習で使用されているアルゴリズムをPythonの基本から丁寧に説明しているので、ここでのPythonサンプルをJuliaで書き換えてみたらいい文法の勉強になると思いました。

第04章以降でいくつか試してみます。
※著作権には詳しくないので、念のためPythonのコードや本に書いてある解説は載せていません。

環境

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.0.3 (2018-12-18)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

浮動小数点の演算

誤差を生む演算の例

• Pythonと同じように書いてみたらエラーになった。

julia> s = 0
0

julia> for i in 0:999
       s += 0.001
       end
ERROR: UndefVarError: s not defined
Stacktrace:
 [1] top-level scope at ./REPL[40]:2 [inlined]
 [2] top-level scope at ./none:0

• 変数sのスコープをglobalに修正。

julia> s = 0
0

julia> for i in 0:999
       global s += 0.001
       end

julia> s
1.0000000000000007

誤差により無限ループとなる例

julia> s = 0
0

julia> while s != 1.
         print(s)
         global s += 0.1
       end

(無限に続く)

誤差を生まない演算方法

• Pythonと同じように書いてみたらエラーになった。

julia> esp = 1e-10
1.0e-10

julia> s = 0
0

julia> while abs(s - 1.) > eps
         println(s)
         global s += 0.1
       end
ERROR: MethodError: no method matching isless(::typeof(eps), ::Float64)
Closest candidates are:
  isless(::Float64, ::Float64) at float.jl:459
  isless(::Missing, ::Any) at missing.jl:66
  isless(::AbstractFloat, ::AbstractFloat) at operators.jl:148
  ...
Stacktrace:
 [1] <(::Function, ::Float64) at ./operators.jl:260
 [2] >(::Float64, ::Function) at ./operators.jl:286
 [3] top-level scope at ./none:0

Juliaには既にeps(Machine epsilon)という関数が用意されているようです。
https://docs.julialang.org/en/v1/manual/integers-and-floating-point-numbers/index.html

julia> eps
eps (generic function with 9 methods)

• eps関数を利用して実行する。

julia> s = 0
0

julia> while abs(s - 1.) > eps()
         println(s)
         global s += 0.1
       end
0
0.1
0.2
0.30000000000000004
0.4
0.5
0.6
0.7
0.7999999999999999
0.8999999999999999

julia> 

演算による桁落ち

qeq.jl

function qeq(a, b, c)
  d = sqrt(b^2 - 4 * a * c)
  ((-b + d) / (2 * a), (-b - d) / (2 * a))
end

答えが正しい場合

julia> include("qeq.jl")
qeq (generic function with 1 method)

julia> qeq(1, 5, 6)
(-2.0, -3.0)

間違っている場合

julia> qeq(1, 1.000000001, 0.000000001)
(-1.0000000272292198e-9, -1.0)

$\sqrt{b^2 - 4ac}$の計算

julia> sqrt(1.000000001^2 - 4 * 1 * 0.000000001)
0.999999999

$\alpha = \frac{-b - sign(b)\sqrt{b^2 - 4ac}}{2a}$
$\beta = \frac{c}{a\alpha}$
の実装

qeq2.jl

function qeq2(a, b, c)
  alpha = (-b - sign(b) * sqrt(b^2 - 4 * a * c)) / (2 * a)
  beta = c / (a * alpha)
  (alpha, beta)
end

実行

julia> include("qeq2.jl")
qeq2 (generic function with 1 method)

julia> qeq2(1, 5 ,6)
(-3.0, -2.0)

julia> qeq2(1, 1.000000001, 0.000000001)
(-1.0, -1.0e-9)

数値範囲の考慮

$softplus(x) = \log(1 + e^x)$
の実装

function softplus(x)
  log(1 + exp(x))
end

実行

julia> include("softplus.jl")
softplus (generic function with 1 method)

julia> softplus(-1)
0.31326168751822286

julia> softplus(0)
0.6931471805599453

julia> softplus(1000)
Inf

$e^{1000}$の計算

julia> exp(1000)
Inf

$softplus2(x) = max(0, x) + log(1 + e^{-|x|})$
の実装

softplus2.jl

function softplus2(x)
  max(0, x) + log(1 + exp(-abs(x)))
end

改良した関数softplus2を使う

julia> include("softplus2.jl")
softplus2 (generic function with 1 method)

julia> softplus2(-1)
0.31326168751822286

julia> softplus2(0)
0.6931471805599453

julia> softplus2(1000)
1000.0

julia> softplus2(-1000)
0.0