4 Numerical Computation¶

In [87]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns
sns.set_context('notebook')

解析的ではなく，反復プロセスによって解の見積もりを更新する方法によって数学的問題を解くアルゴリズムを紹介．

「ゼロから作るディープラーニング」第4章も参照．

4.1 Overflow and Underflow¶

計算機上で連続的な計算を実行する際の基本的な困難さは、有限長のビットパターンで無限に多くの実数を表現する必要があることである．これは，ほとんどすべての実数について，計算機内の数値を表すときに近似誤差が発生することを意味する．多くの場合，これはちょうど丸め誤差である．丸め誤差の累積を最小限に抑えるように設計されていないと，実際に多くの操作で混乱し，理論上で動作するアルゴリズムが失敗する可能性がある．

丸め誤差の例:
- underflow (アンダーフロー): ゼロに近い数値をゼロに丸めると発生する．
- overflow (オーバーフロー): 大きな数値を $\infty$ または $-\infty$ として近似するときに発生する．

例: softmax¶

$\mathrm{softmax}(\mathbf{x})_i = \frac{\exp(x_i)}{\sum_{j=1}^n \exp(x_j)} \tag{dl 4.1}$

$x_i$ がある定数 $c$ に等しいとする．

In [57]:

c = 10
x = np.full((1, 5), c)
x

/Users/ryo-t/.pyenv/versions/anaconda3-4.3.0/lib/python3.6/site-packages/numpy/core/numeric.py:301: FutureWarning: in the future, full((1, 5), 10) will return an array of dtype('int64')
  format(shape, fill_value, array(fill_value).dtype), FutureWarning)

Out[57]:

array([[ 10.,  10.,  10.,  10.,  10.]])

In [58]:

np.exp(x)

Out[58]:

array([[ 22026.46579481,  22026.46579481,  22026.46579481,  22026.46579481,
         22026.46579481]])

In [59]:

# softmax(x)
np.exp(x) / np.exp(x).sum()

Out[59]:

array([[ 0.2,  0.2,  0.2,  0.2,  0.2]])

🙆

$c$ がデカイと $\exp(c)$ がオーバーフローし，softmax は $\infty/\infty$ で undefined．

In [60]:

c = 10e2
x = np.full((1, 5), c)
x

Out[60]:

array([[ 1000.,  1000.,  1000.,  1000.,  1000.]])

In [61]:

np.exp(x)

/Users/ryo-t/.pyenv/versions/anaconda3-4.3.0/lib/python3.6/site-packages/ipykernel/__main__.py:1: RuntimeWarning: overflow encountered in exp
  if __name__ == '__main__':

Out[61]:

array([[ inf,  inf,  inf,  inf,  inf]])

In [62]:

# softmax(x)
np.exp(x) / np.exp(x).sum()

/Users/ryo-t/.pyenv/versions/anaconda3-4.3.0/lib/python3.6/site-packages/ipykernel/__main__.py:2: RuntimeWarning: overflow encountered in exp
  from ipykernel import kernelapp as app
/Users/ryo-t/.pyenv/versions/anaconda3-4.3.0/lib/python3.6/site-packages/ipykernel/__main__.py:2: RuntimeWarning: invalid value encountered in true_divide
  from ipykernel import kernelapp as app

Out[62]:

array([[ nan,  nan,  nan,  nan,  nan]])

🙀

$c$ が小さいと $\exp(c)$ がアンダーフローし，softmax は $0/0$ で undefined．

In [63]:

c = -10e2
x = np.full((1, 5), c)
x

Out[63]:

array([[-1000., -1000., -1000., -1000., -1000.]])

In [64]:

np.exp(x)

Out[64]:

array([[ 0.,  0.,  0.,  0.,  0.]])

In [65]:

# softmax(x)
np.exp(x) / np.exp(x).sum()

/Users/ryo-t/.pyenv/versions/anaconda3-4.3.0/lib/python3.6/site-packages/ipykernel/__main__.py:2: RuntimeWarning: invalid value encountered in true_divide
  from ipykernel import kernelapp as app

Out[65]:

array([[ nan,  nan,  nan,  nan,  nan]])

🙀

解決策: $\mathrm{softmax}(\mathbf{z})$ , where $\mathbf{z} = \mathbf{x} - \mathrm{max}_i x_i$ を評価¶

In [158]:

c = 10e2
x = np.full((1, 5), c)
x

Out[158]:

array([[ 1000.,  1000.,  1000.,  1000.,  1000.]])

In [159]:

z = x - x.max()
z

Out[159]:

array([[ 0.,  0.,  0.,  0.,  0.]])

In [160]:

# softmax(z)
np.exp(z) / np.exp(z).sum()

Out[160]:

array([[ 0.2,  0.2,  0.2,  0.2,  0.2]])

🙆

In [168]:

c = -10e2
x = np.full((1, 5), c)
x

Out[168]:

array([[-1000., -1000., -1000., -1000., -1000.]])

In [169]:

z = x - x.max()
z

Out[169]:

array([[ 0.,  0.,  0.,  0.,  0.]])

In [170]:

# softmax(z)
np.exp(z) / np.exp(z).sum()

Out[170]:

array([[ 0.2,  0.2,  0.2,  0.2,  0.2]])

🙆

In [161]:

x = np.array([1, 2, 3, 4, 5])
z = x - x.max()
z

Out[161]:

array([-4, -3, -2, -1,  0])

In [167]:

np.exp(x) / np.exp(x).sum(), np.exp(z) / np.exp(z).sum()

Out[167]:

(array([ 0.01165623,  0.03168492,  0.08612854,  0.23412166,  0.63640865]),
 array([ 0.01165623,  0.03168492,  0.08612854,  0.23412166,  0.63640865]))

🙆

普通はこういうトリックが必要だが，この本で説明されているさまざまなアルゴリズムの実装に関わるすべての数値的な考慮事項は明示的に記述されていない．Theano など低レベルライブラリの開発者は，deep learning アルゴリズムを実装する際に数値問題を念頭に置くべき．でもこの本のほとんどの読者は，低レベルライブラリに頼ってればOK牧場．

4.2 Poor Conditioning¶

conditioning: 入力のわずかな変化に対して，関数がどれほど急速に変化するかを指す．入力の丸め誤差が出力に大きな変化をもたらす可能性があるため，入力がわずかに動いたときに急激に変化する関数は科学的計算には問題がある．

関数 $f(\mathbf{x}) = \mathbf{A}^{-1}\mathbf{x}$ を考える． $\mathbf{A} \in \mathbb{R}^{n \times n}$ が固有値を持つとき，condition number は

$\max_{i,j} \left| \frac{\lambda_i}{\lambda_j} \right| \tag{dl 4.2}$

In [172]:

a = np.array([[1, 0, -1], [0, 1, 0], [1, 0, 1]])
a

Out[172]:

array([[ 1,  0, -1],
       [ 0,  1,  0],
       [ 1,  0,  1]])

In [173]:

np.linalg.cond(a)

Out[173]:

1.4142135623730951

In [178]:

a = np.diag([1, 10, 100])
a

Out[178]:

array([[  1,   0,   0],
       [  0,  10,   0],
       [  0,   0, 100]])

In [179]:

np.linalg.cond(a)

Out[179]:

100.0

この値がデカイと，逆行列は入力の誤差を増幅する．実際には，誤差は逆行列自体の数値誤差によってさらに複合化される．

4.3 Gradient-Based Optimization¶

ほとんどの deep learning アルゴリズムには，ある種の最適化が必要．最適化とは， $x$ を変更することによって，ある関数 $f(x)$ を最小化または最大化するタスクのこと． $f(x)$ を最小化するという点では，ほとんどの最適化の問題に言及する．

最小化または最大化したい関数は目的関数 (objective function) または最適化基準 (criterion) と呼ばれる．最小化するときには，コスト関数 (cost function)，損失関数 (loss function)，またはエラー関数 (error function) とも呼ばれる．本書では、これらの用語を同じ意味で使用していますが、一部の機械学習の出版物では、これらの用語のいくつかに特別な意味が割り当てられています。

読者はすでに微積分に精通していると仮定し，微積分の概念が最適化にどのように関係しているかを簡単に解説する．

数値微分¶

$\frac{df(x)}{dx} = \lim_{h \to 0} \frac{f(x+h) - f(x-h)}{2h}$

In [69]:

def numerical_diff(f, x):
    h = 1e-4
    return (f(x+h) - f(x-h)) / (2*h)

$f_1(x) = 0.01x^2 + 0.1x$

In [101]:

def function_1(x):
    return 0.01*x**2 + 0.1*x

In [137]:

def tangent_line(f, x):
    d = numerical_diff(f, x)
    y = f(x) - d*x
    return lambda t: d*t + y

In [138]:

x = np.arange(0.0, 20.0, 0.1)
y = function_1(x)
plt.xlabel("x")
plt.ylabel("f(x)")

tf = tangent_line(function_1, 5)
y2 = tf(x)

plt.plot(x, y)
plt.plot(x, y2)

Out[138]:

[<matplotlib.lines.Line2D at 0x11ba486a0>]

偏微分¶

$f_2(x_0,x_1) = x_0^2 + x_1^2 \tag{zero 4.6}$

In [113]:

def function_2(x):
    if x.ndim == 1:
        return np.sum(x**2)
    else:
        return np.sum(x**2, axis=1)

$\frac{\partial{f_2}}{\partial{x_0}}$

In [120]:

def function_tmp1(x0):
    return x0**2 + 0**2  # x_1^2 は定数なので何でもいいし結果も変わらない

numerical_diff(function_tmp1, 3.0)

Out[120]:

6.000000000012662

$\frac{\partial{f_2}}{\partial{x_1}}$

In [123]:

def function_tmp1(x1):
    return 0**2 + x1**2  # x_0^2 は定数なので何でもいいし結果も変わらない

numerical_diff(function_tmp1, 4.0)

Out[123]:

7.999999999999119

勾配 (gradient)¶

$x_0$ と $x_1$ の偏微分をまとめて計算したいとする． $(x_0,x_1)$ の両方の偏微分をまとめて， $(\frac{\partial{f}}{\partial{x_0}}, \frac{\partial{f}}{\partial{x_1}})$ として計算することを考える． $(\frac{\partial{f}}{\partial{x_0}}, \frac{\partial{f}}{\partial{x_1}})$ のように，すべての変数の偏微分をベクトルとしてまとめたものを勾配 (gradient) と呼ぶ．

In [132]:

def numerical_gradient(f, x):
    h = 1e-4 # 0.0001
    grad = np.zeros_like(x) # x と同じ形状の零ベクトル
    
    for idx in range(x.size):
        tmp_val = x[idx]
        # f(x+h) の計算
        x[idx] = float(tmp_val) + h
        fxh1 = f(x)
        
        # f(x-h) の計算
        x[idx] = tmp_val - h 
        fxh2 = f(x)
        
        grad[idx] = (fxh1 - fxh2) / (2*h)
        x[idx] = tmp_val # 値を元に戻す
        
    return grad

点 $(3,4)$ での勾配:

In [133]:

numerical_gradient(function_2, np.array([3.0, 4.0]))

Out[133]:

array([ 6.,  8.])

勾配をベクトルだと思っていろんな点で求めて図示すると:

In [136]:

def _numerical_gradient_no_batch(f, x):
    h = 1e-4 # 0.0001
    grad = np.zeros_like(x) # x と同じ形状の零ベクトル
    
    for idx in range(x.size):
        tmp_val = x[idx]
        # f(x+h) の計算
        x[idx] = float(tmp_val) + h
        fxh1 = f(x)
        
        # f(x-h) の計算
        x[idx] = tmp_val - h 
        fxh2 = f(x)
        
        grad[idx] = (fxh1 - fxh2) / (2*h)
        x[idx] = tmp_val # 値を元に戻す
        
    return grad


def numerical_gradient(f, X):
    if X.ndim == 1:
        return _numerical_gradient_no_batch(f, X)
    else:
        grad = np.zeros_like(X)
        
        for idx, x in enumerate(X):
            grad[idx] = _numerical_gradient_no_batch(f, x)
        
        return grad

x0 = np.arange(-2, 2.5, 0.25)
x1 = np.arange(-2, 2.5, 0.25)
X, Y = np.meshgrid(x0, x1)

X = X.flatten()
Y = Y.flatten()

grad = numerical_gradient(function_2, np.array([X, Y]))

plt.figure()
plt.quiver(X, Y, -grad[0], -grad[1], angles="xy", color="#666666")
plt.xlim([-2, 2])
plt.ylim([-2, 2])
plt.xlabel('x0')
plt.ylabel('x1')

Out[136]:

<matplotlib.text.Text at 0x11baed860>

各地点において関数の値を最も減らす方向を示すのが勾配．

勾配法¶

機械学習の問題の多くは，学習の際に最適なパラメータを探索する．最適なパラメータとは，損失関数が最小値を取るときのパラメータの値．勾配をうまく利用して関数の最小値を探すのが勾配法．

ただし，各地点において関数の値を最も減らす方向を示すのが勾配であるため，勾配が指す先が本当に関数の最小値なのかどうかは保証されない．

$\begin{split} x_0 &= x_0 - \eta\frac{\partial{f}}{\partial{x_0}}\\ x_1 &= x_1 - \eta\frac{\partial{f}}{\partial{x_1}} \end{split} \tag{zero 4.7}$

より一般的に書くと

$\mathbf{x}' = \mathbf{x} - \epsilon\nabla_\mathbf{x}f(\mathbf{x}) \tag{dl 4.5}$

$\eta$ ( $\epsilon$ ) は学習率 (learning rate) と呼ばれる．ニューラルネットワークの学習においては，学習率の値を変更しながら，正しく学習できているかどうか，確認作業を行うのが一般的．

In [145]:

def gradient_descent(f, init_x, lr=0.01, step_num=100):
    x = init_x
    x_history = []

    for i in range(step_num):
        x_history.append(x.copy())

        grad = numerical_gradient(f, x)
        x -= lr * grad

    return x, np.array(x_history)

初期値 $(-3.0, 4.0)$ ，学習率 $\eta = 0.1$ :

In [151]:

init_x = np.array([-3.0, 4.0])    

lr = 0.1
step_num = 20
x, x_history = gradient_descent(function_2, init_x, lr=lr, step_num=step_num)

plt.plot([-5, 5], [0,0], '--b')
plt.plot([0,0], [-5, 5], '--b')
plt.plot(x_history[:,0], x_history[:,1], 'o')

plt.xlim(-3.5, 3.5)
plt.ylim(-4.5, 4.5)
plt.xlabel("X0")
plt.ylabel("X1")

Out[151]:

<matplotlib.text.Text at 0x11c55b128>

🙆

初期値 $(-3.0, 4.0)$ ，学習率 $\eta = 0.01$ :

In [152]:

init_x = np.array([-3.0, 4.0])    

lr = 0.01
step_num = 20
x, x_history = gradient_descent(function_2, init_x, lr=lr, step_num=step_num)

plt.plot([-5, 5], [0,0], '--b')
plt.plot([0,0], [-5, 5], '--b')
plt.plot(x_history[:,0], x_history[:,1], 'o')

plt.xlim(-3.5, 3.5)
plt.ylim(-4.5, 4.5)
plt.xlabel("X0")
plt.ylabel("X1")

Out[152]:

<matplotlib.text.Text at 0x11c6d4438>

小さすぎた 🙀

初期値 $(-3.0, 4.0)$ ，学習率 $\eta = 1.0$ :

In [153]:

init_x = np.array([-3.0, 4.0])    

lr = 1.0
step_num = 20
x, x_history = gradient_descent(function_2, init_x, lr=lr, step_num=step_num)

plt.plot([-5, 5], [0,0], '--b')
plt.plot([0,0], [-5, 5], '--b')
plt.plot(x_history[:,0], x_history[:,1], 'o')

plt.xlim(-3.5, 3.5)
plt.ylim(-4.5, 4.5)
plt.xlabel("X0")
plt.ylabel("X1")

Out[153]:

<matplotlib.text.Text at 0x11c81f080>

デカすぎた 🙀

line search では幾つかの $\epsilon$ に対して $f(\mathbf{x} - \epsilon\nabla_\mathbf{x}f(\mathbf{x}))$ を計算して，目的関数が最小となるような $\epsilon$ を選ぶ．

4.3.1 Beyond the Gradient: Jacobian and Hessian Matrices¶

4.4 Constrained Optimization¶

4.5 Example: Linear Least Squares¶

In [254]:

A = np.eye(2)
b = np.array([1, 1]).reshape((2, 1))

def f(x):
    if x.shape[1] == 1:
        return np.linalg.norm(A.dot(x) - b)**2 / 2
    return np.linalg.norm(A.dot(x) - b, axis=0)**2 / 2

In [275]:

step_size = 0.01
tolerance = 0.01

x = np.array([[2.0], [2.0]])
history = []

grad = A.T.dot(A).dot(x) - A.T.dot(b)
while np.linalg.norm(grad) > tolerance:
    history.append(x.copy())
    x -= step_size * grad
    grad = A.T.dot(A).dot(x) - A.T.dot(b)

In [279]:

len(history)

Out[279]:

In [278]:

history[-1]

Out[278]:

(array([[ 1.00712059],
        [ 1.00712059]]), 493)

x = [1, 1] で f(x) が最小なのでOK牧場．