# Juliaによる数値計算

## 浮動小数点の演算

- **浮動小数点**:
    - 小数をコンピュータ上で表現するための数
    - 固定長の仮数部と指数部をビット2進数により表現する
    - 2進数で表現するため、ほとんどの小数は無限小数となり、有限のリソースで計算を行うコンピュータ上では誤差が生じる

簡単な例として、0.001 を 1000 個加算してみると以下のようになる


```julia
a = 0

for i = 1:1000
    a += 0.001
end

println(a)
```

    1.0000000000000007


数学的には 0.001 を 1000 個加算すれば 1 になるはずだが、上記のように $7.0\times10^{-16}$ の誤差が生じてしまっている

これが2進数表現で小数を表現する場合に生じてしまう誤差であり、**計算機イプシロン** と呼ばれる

Julia において、計算機イプシロンは `eps` 関数を用いて取得することが出来る


```julia
# Float64型浮動小数点における計算機イプシロン
eps(Float64)
```




    2.220446049250313e-16



上記のように浮動小数点では誤差が生じるが、その誤差は非常に小さく、実用上は近似計算してしまって問題ないことがほとんどである

しかしながら、このような内部表現による誤差があることは、常に意識しておかなければ予想外の問題が発生しうるため注意が必要である

ここではまず、条件分岐に関する注意点を示す

浮動小数点に関する条件式では数値誤差を十分に考慮しなければならず、特に数値比較で `==` のような完全一致比較を行うのは悪手であることが多い


```julia
s = 0
# s を 0.0 から 0.1 ずつ加算していく
for i = 1:20
    # s が 1.0 になったら終了
    ## => 浮動小数点の内部表現により s が正確に 1.0 になることはないため、この条件式は永遠に true にならない
    if s == 1.0
        break
    end
    println(s)
    s += 0.1
end
```

    0
    0.1
    0.2
    0.30000000000000004
    0.4
    0.5
    0.6
    0.7
    0.7999999999999999
    0.8999999999999999
    0.9999999999999999
    1.0999999999999999
    1.2
    1.3
    1.4000000000000001
    1.5000000000000002
    1.6000000000000003
    1.7000000000000004
    1.8000000000000005
    1.9000000000000006


上記のように、浮動小数点の比較に `==` や `!=` のような完全一致演算子を用いるのは意図と異なる結果になることが多い

こういった場合には、十分に小さな正の数を用意して、誤差がその数値以内になるかどうかをチェックするという手法がよく用いられる

例えば、Juliaにおいては概ね17桁程度の精度の浮動小数点を扱う事ができるため、$10^{-10}$ のような十分に小さな値を設定して数値比較を行えば、意図通りの結果が得られるはずである


```julia
# 十分に小さな正の値
## ϵ: \epsilon <Tab>
ϵ = 1e-10

s = 0

# s を 0.0 から 0.1 ずつ加算していく
for i = 1:20
    # s が 1.0（付近）になったら終了したい
    ## => |s-1.0| が ϵ より小さくなったら終了
    if abs(s - 1.0) < ϵ
        break
    end
    println(s)
    s += 0.1
end
```

    0
    0.1
    0.2
    0.30000000000000004
    0.4
    0.5
    0.6
    0.7
    0.7999999999999999
    0.8999999999999999


上記では $\epsilon$ を $10^{-10}$ と設定したが、この値をどのようにして決めるのが良いのかは、解く問題によるため一概には言えない

理論的に値が決められることもあれば、実験的にうまくいく値を採用することもある

### 演算による桁落ち

次に2次方程式の解を求める関数を考えてみる

2次方程式 $ax^2 + bx + c = 0 \quad (a\neq0)$ の解は

$$
x = \frac{-b\pm\sqrt{b^2 - 4ac}}{2a}
$$

で与えられる

これを計算する関数を定義してみる


ここでは話を簡単にするため、解が実数になる場合、つまり $b^2-4ac \geq 0$ の場合のみを考えることにする


```julia
# 2次方程式の解を求める関数
qeq(a, b, c) = begin
    d = sqrt(b^2 - 4*a*c)
    # 2つの解をTupleで返す
    ((-b + d) / 2a, (-b - d) / 2a)
end

"""
test: x^2 + 5x + 6 = 0
    a = 1
    b = 5
    c = 6
"""
# answer: (x + 2)(x + 3) = 0 より
## x = -2, -3 になるはず
qeq(1, 5, 6)
```




    (-2.0, -3.0)



上記の通り、理論値通り正確に計算できた

続いて、2次方程式 $x^2 + 1.000000001x + 0.000000001 = 0$ について考えてみる


```julia
"""
test: x^2 + 1.000000001x + 0.000000001 = 0
    a = 1
    b = 1.000000001
    c = 0.000000001
"""
# answer: (x + 1)(x + 0.000000001) = 0 より
## x = -1, -0.000000001 になるはず
qeq(1, 1.000000001, 0.000000001)
```




    (-1.0000000272292198e-9, -1.0)



上記の通り、理論値と異なる値が算出された

これは浮動小数点の **桁落ち** と呼ばれる現象により、誤差が生じているためである

今回の場合、分子の $-b + \sqrt{b^2 - 4ac}$ が 0 に非常に近いため、桁落ちが生じている

実際に $\sqrt{b^2 - 4ac}$ を計算してみると以下のようになる


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




    0.999999999



上記の結果と $-b = -1.000000001$ を加算すると、小数点第8位までは消えてしまう

つまり、もともとは約17桁まで信用できた数字が、9桁分までの数字が消えてしまうことにより、信用できる桁数が残り8桁になってしまうということである

これが桁落ちと呼ばれる現象である


```julia
-1.000000001 + sqrt(1.000000001^2 - 4 * 1 * 0.000000001)
```




    -2.0000000544584395e-9



このように、絶対値の近い減算処理では有効桁数を大きく失う可能性があるため危険である

#### 桁落ちの回避方法
桁落ちを回避するためには、絶対値が近い数値同士の減算処理が発生しそうな計算を回避できるように、数式を同値な式で置き換えれば良い

上記2次方程式の場合、$\sqrt{}$ の計算結果は必ず正になるはずであるため、絶対値の近い数値同士の減算処理が発生する可能性は $b$ の符号によって分けて考える必要がある

$b\geq0$ のときは、分子の $\pm$ が $+$ のときに、近い数値同士の減算処理が発生しうる

逆に $b\gt0$ のときは、分子の $\pm$ が $-$ のときに、近い数値同士の減算処理が発生しうる

これらをまとめて表現するために符号関数 $\rm{sign}(x)$ という関数を導入する

$$
\rm{sign}(x) = \left\{
\begin{align}
-1\quad&(x\lt0) \\
0\quad&(x=0) \\
1\quad&(x\gt0) \\
\end{align}
\right.
$$

符号関数を用いて、直接計算しても問題ない（絶対値の近い減算処理が**発生しない**）式を考えると、以下のように表される

$$
\frac{-b-\rm{sgin}(b)\sqrt{b^2 - 4ac}}{2a}
$$

次に直接計算すると問題のある（絶対値の近い減算処理が発生しうる）方の解を計算する方法を考える

ここで、2次方程式 $ax^2 + bx + c = 0$ の解を $\alpha, \beta$ と置くと、次の式が成り立つことが知られている

$$
\alpha\beta = \frac{c}{a}
$$

このことは、以下のように実際に計算してみるとすぐに分かる

$$
\begin{align}
\alpha\beta &= \frac{-b + \sqrt{b^2 - 4ac}}{2a}\times\frac{-b - \sqrt{b^2 - 4ac}}{2a} \\
&= \frac{(-b)^2 - (\sqrt{b^2 - 4ac})^2}{4a^2} \\
&= \frac{b^2 - (b^2 - 4ac)}{4a^2} \\
&= \frac{c}{a}
\end{align}
$$

上記を利用することで、直接計算に問題のある方の解 $\beta$ を、直接計算しても問題ない方の解 $\alpha$ で表すことが出来るようになる

つまり、

$$
\begin{align}
\alpha &= \frac{-b - \rm{sign}(b)\sqrt{b^2 - 4ac}}{2a} \\
\beta &= \frac{c}{a\alpha}
\end{align}
$$

という形で計算すれば桁落ちを回避できるはずである


```julia
# 2次方程式の解を求める関数（桁落ち回避版）
qeq(a, b, c) = begin
    # Julia には sign 関数は事前定義されているためそれを利用
    alpha = (-b - sign(b) * sqrt(b^2 - 4*a*c)) / 2a
    beta = c / (a * alpha)
    (alpha, beta)
end

# test: x^2 + 1.000000001x + 0.000000001 = 0
qeq(1, 1.000000001, 0.000000001)
```




    (-1.0, -1.0e-9)



上記の通り、前回よりも正確な値を計算することが出来た

このように、誤差を比較的小さく抑えられる方法を **数値的に安定な解法** と呼ぶ

### 数値範囲の考慮
次に、ソフトプラス関数と呼ばれる次のような関数を考える

$$
\rm{softplus}(x) = \ln(1+e^x)
$$

この関数は $x$ が十分に大きい場合、$\rm{softplus}(x) \approx x$ で近似される

実際 $e^x$ は $x$ が大きくなると急激に大きくなるため、$1 + e^x \approx e^x$ と近似できるようになり、$\rm{softplus}(x) \approx \ln e^x = x$ となる

まずは、この関数を定義そのままで記述してみる


```julia
softplus(x) = log(1 + exp(x))
    
# test
softplus(-1) |> println
softplus(0) |> println
softplus(1000) |> println
```

    0.31326168751822286
    0.6931471805599453
    Inf


上記のように $x=-1$ と $x=0$ のときは順調に計算出来ていたが、$x=1000$ のときは `Inf`（無限大）となってしまった

これは、Julia で扱える浮動小数点の範囲を超えて大きな値になったことを意味する

しかしながら、前述の考察によると $\rm{softplus}(1000) \approx 1000$ になるはずであるため、浮動小数点型で十分に計算できる範囲のように思われる

それでも `Inf` になってしまうのは、`exp(1000)` の計算時点で浮動小数点型の範囲を超えてしまうためである


```julia
exp(1000)
```




    Inf



ここでの問題は、数学的に考えれば計算結果はさほど大きくない（Julia の浮動小数点型として十分に取り扱い可能）にも関わらず、計算過程で `Inf` が発生してその後の計算も全て `Inf` になってしまうという点である

この問題を克服するには、元の式を同値な式に変形し、計算過程で `Inf` が発生しないようにすることである

試しに以下のように式変形してみる

$$
\begin{align}
\ln(1 + e^x) &= \ln\{e^x(x^{-x} + 1)\} \\
&= \ln{e^x} + \ln(e^{-x} + 1) \\
&= x + \ln(e^{-x} + 1)
\end{align}
$$

上記のように変形すれば、$x$ が十分に大きくなると $e^{-x}\right0$ となるため `Inf` が発生する問題は回避できる

しかし、今度は $x$ が十分に小さくなる場合に $e^{-x}\right\inf$ となってしまい、この場合は元の式の方が良かったことになる

したがって、どこかの $x$ を境目にして2つの式を分けるのが良さそうである

ここでは $x=0$ を境目にして、以下のように式変形してみる

$$
\begin{align}
\rm{softplus}(x) &= \left\{
\begin{array}{cc}
\ln(1 + e^x) & (x\lt0) \\
x + \ln(1 + e^{-x}) & (x\geq0)
\end{array}
\right. \\
&= \max(0,x) + \ln(1 + e^{-|x|})
\end{align}
$$

ここで $\max(a,b)$ とは $a$, $b$ の内の大きい方を返す関数である

上式に従って、`softplus` 関数を修正してみる


```julia
softplus(x) = max(0, x) + log(1 + exp(-abs(x)))

# test
softplus(-1) |> println
softplus(0) |> println
softplus(1000) |> println
softplus(-1000) |> println
```

    0.31326168751822286
    0.6931471805599453
    1000.0
    0.0


上記の通り、今度は $x=1000$ の場合も $x=-1000$ の場合も問題なく計算できている

以上、コンピュータによる数値計算で気をつけるべき点を、例を挙げながら見てきたが、ポイントをまとめると以下のようになる

- 原則として、浮動小数点型の比較には `==` や `!=` のような完全一致演算子を用いない
    - 十分な小さな正の数を設定して、「比較対象との差がその数値以内になる」という条件式を用いる方法がよく採用される
- 数値が近いものの減算処理は、桁落ち現象により有効桁数が失われるため注意
    - 数値が近いものの減算処理が入らないように、同値な式に変形して計算する
- 計算過程で `Inf` や `-Inf` が生じてしまう計算に注意
    - 浮動小数点型で扱える範囲の数値しか現れないように式変形して計算する


```julia

```

f15063486d6f66130d73f4442984b318448680cb
ipynb
Jupyter Notebook
02_machine-learning/03-01_float.ipynb
amenoyoya/julia_ml-tuto
9c0be0923ea00ca4d1d51c0c6f61f6f2748232be
MIT
02_machine-learning/03-01_float.ipynb
amenoyoya/julia_ml-tuto
9c0be0923ea00ca4d1d51c0c6f61f6f2748232be
MIT
02_machine-learning/03-01_float.ipynb
amenoyoya/julia_ml-tuto
9c0be0923ea00ca4d1d51c0c6f61f6f2748232be
MIT
Qwen/Qwen-72B
1. YES
2. YES