からっぽのしょこ

読んだら書く!書いたら読む!同じ事は二度調べ(たく)ない

ステップ27:sin関数のテイラー展開【ゼロつく3のノート(数学)】

はじめに

 『ゼロから作るDeep Learning 3』の初学者向け攻略ノートです。『ゼロつく3』学習の補助となるように適宜解説を加えていきます。本と一緒に読んでください。

 本で登場する数学的な内容をもう少し深掘りして解説していきます。

 この記事は、主にステップ27「テイラー展開の微分」を補足する内容です。
 $\sin$関数のマクローリン展開とテイラー展開をPythonで実装します。

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

・テイラー展開

 テイラー展開とは、任意の関数$f(x)$を次の多項式で近似する方法です。

$$ f(x) = f(a) + f'(a) (x - a) + \frac{1}{2!} f''(a) (x - a)^2 + \frac{1}{3!} f'''(a) (x - a)^3 + \cdots \tag{27.1} $$

 特に、$a = 0$のときをマクローリン展開と呼びます。

$$ f(x) = f(0) + f'(0) x + \frac{1}{2!} f''(0) x^2 + \frac{1}{3!} f'''(0) x^3 + \cdots \tag{27.2} $$

 $f'(x)$は$f(x)$の$x$での微分を表し、$f''(x)$は$f'(x)$を更に$x$で微分したもので2階微分と呼びます。同様に、$f'''(x)$なら3階微分です。また、$f'(a)$は$f'(x)$に$x = a$を代入したものです。

 ここでは、理論的な内容は省略して、プログラムによる計算とグラフを通して理解を深めていきます。

 次のライブラリを利用します。

# 利用するライブラリ
import numpy as np
import math
import matplotlib.pyplot as plt
import matplotlib.animation as animation

 階乗$x!$の計算にmathfactorial()を使います。matplotlibanimationは、アニメーション(gif画像)作成用のモジュールです。処理がややこしいので手元で実行するかはお好みでどうぞ。

・$\sin$関数のマクローリン展開の確認

 まずは、簡単な例として$\sin$関数のマクローリン展開を行います。

 $\sin (x)$をグラフで確認しましょう。$\sin (x)$はnp.sin()で計算できます。

# x軸の値を設定
x = np.arange(-10.0, 10.0, 0.01)

# sin(x)を計算
sin_x = np.sin(x)

# 確認
print(x[:10])
print(sin_x[:10])
[-10.    -9.99  -9.98  -9.97  -9.96  -9.95  -9.94  -9.93  -9.92  -9.91]
[0.54402111 0.53560333 0.527132   0.51860795 0.51003204 0.50140513
 0.49272808 0.48400175 0.47522703 0.46640478]


 グラフを描きます。

# 作図
plt.figure(figsize=(12, 3)) # 画像サイズ
plt.plot(x, sin_x, label='sin(x)') # sin(x)
plt.xlabel('x') # x軸ラベル
plt.ylabel('f(x)') # y軸ラベル
plt.grid() # グリッド線
plt.legend() # 凡例
plt.show()

$\sin$

 $\sin$関数自体にはここでは深入りしません。

 ただし、$\sin (x)$を微分すると$\cos$関数が登場するので、$\cos (x)$もグラフで見ておきましょう。$\cos (x)$はnp.cos()で計算できます。

# cos(x)を計算
cos_x = np.cos(x)

# 確認
print(x[:10])
print(cos_x[:10])
[-10.    -9.99  -9.98  -9.97  -9.96  -9.95  -9.94  -9.93  -9.92  -9.91]
[-0.83907153 -0.8444697  -0.84978342 -0.85501216 -0.8601554  -0.86521263
 -0.87018334 -0.87506703 -0.87986321 -0.88457141]


 $\sin(x)$と$\cos(x)$を重ねて作図します。

# 作図
plt.figure(figsize=(12, 3)) # 画像サイズ
plt.plot(x, sin_x, label='sin(x)') # sin(x)
plt.plot(x, cos_x, label='cos(x)') # cos(x)
plt.xlabel('x') # x軸ラベル
plt.ylabel('f(x)') # y軸ラベル
plt.grid() # グリッド線
plt.legend() # 凡例
plt.show()

$\sin$と$\cos$

 $\sin(0) = 0$、$\cos(0) = 1$ということだけ覚えておけば今回の話は大丈夫です。

 では、マクローリン展開($a = 0$の場合)について考えていきます。値が0のオブジェクトaを作成しておきます。

# aを指定
a = 0.0


 無限に続く式(27.2)に関して、1つ目の項だけで$f(x)$を近似することを考えます。

$$ f(x) \approx f(0) $$

 これを0次近似と呼びます。$\approx$は近似を表す記号です。両辺は等しくないので$=$では繋げません。
 $f(x) = \sin(x)$とすると、$\sin(0) = 0$なので

$$ \sin(x) \approx \sin(0) = 0 $$

となります。

 実際に計算してみましょう。

# 0次近似を計算
approx_x0 = np.sin(0)
print(approx_x0)
0.0


 $\sin (x)$のグラフと重ねて確認しましょう。引数のprop={"family":"MS Gothic"}matplotlibで作成するグラフに日本語を表示するためのおまじないです。

# 作図
plt.figure(figsize=(10, 5))
plt.plot(x, sin_x, label='sin(x)') # 元の関数
plt.scatter(a, approx_x0, color='orange', label='0次近似') # 0次近似
plt.grid()
plt.legend(prop={"family":"MS Gothic"})
plt.show()

0次近似

 $x = a = 0$のとき、$\sin(x) = \sin(0) = 0$となります。つまり、$x = 0$の点では等しいですが、近似しているとは言えませんね。

 続いて、式(27.2)に関して2つ目の項まで考えてみます。

$$ f(x) \approx f(0) + f'(0) x $$

 これを1次近似と呼びます。
 $\sin(x)$の1階微分は$\cos(x)$です。$\cos(0) = 1$なので

$$ \begin{aligned} \sin(x) &\approx \sin(0) + \cos(0) x \\ &= 0 + 1 x \\ &= x \end{aligned} $$

となります。

 実際に計算してグラフで確認します。

# 1次近似を計算
approx_x1 = approx_x0 + np.cos(a) * (x - a)
#approx_x1 = approx_x0 + np.cos(a) * x
#approx_x1 = x

 上から、テイラー展開の式(2.71)、マクローリン展開の式(2.72)、今求めた式の3パターンで計算しています。当然、どれも同じ結果になります。

 グラフを描きます。

# 作図
plt.figure(figsize=(10, 5))
plt.plot(x, sin_x, label='sin(x)') # 元の関数
plt.plot(x, approx_x1, label='1次近似') # 1次近似
plt.grid()
plt.legend(prop={"family":"MS Gothic"})
#plt.xlim(-1.0, 1.0) # x軸の表示範囲
plt.ylim(-4.0, 4.0) # y軸の表示範囲
plt.show()

1次近似

 1次関数(直線)で近似しようとしています。

 グラフを表示する範囲を狭めて確認しましょう。

1次近似

 $x = 0$の点を中心にわずかな範囲で近似できています。

 同様に、3つ目の項まで含めてみましょう。

$$ f(x) \approx f(0) + f'(0) x + \frac{1}{2!} f''(0) x^2 $$

 これは2次近似ですね。
 $\sin(x)$の2階微分は$- \sin(x)$です。$\sin(0) = 0$なので

$$ \begin{aligned} \sin(x) &\approx \sin(0) + \cos(0) x - \frac{1}{2!} \sin(0) x^2 \\ &= 0 + 1 x - 0 \\ &= x \end{aligned} $$

となります。1次近似と同じ式になりました。グラフも同じなので省略します。

 3次近似を計算します。

$$ f(x) \approx f(0) + f'(0) x + \frac{1}{2!} f''(0) x^2 + \frac{1}{3!} f'''(0) x^3 $$

 $\sin(x)$の3階微分は$- \cos(x)$です。$\cos(0) = 1$なので

$$ \begin{aligned} \sin(x) &\approx \sin(0) + \cos(0) x - \frac{1}{2!} \sin(0) x^2 - \frac{1}{3!} \cos(0) x^3 \\ &= 0 + x - 0 - \frac{1}{3!} x^3 \\ &= x - \frac{x^3}{3!} \end{aligned} $$

となります。

 同様にやってみます。

# 3次近似を計算
approx_x3 = approx_x1 - 1.0 / math.factorial(3) * np.cos(a) * (x - a)**3
#approx_x3 = approx_x1 - 1.0 / math.factorial(3) * np.cos(a) * x**3
#approx_x3 = approx_x1 - x**3 / math.factorial(3)
# 作図
plt.figure(figsize=(10, 5))
plt.plot(x, sin_x, label='sin(x)') # 元の関数
plt.plot(x, approx_x1, linestyle='dotted', label='1次近似') # 1次近似
plt.plot(x, approx_x3, label='3次近似') # 3次近似
plt.grid()
plt.legend(prop={"family":"MS Gothic"})
#plt.xlim(-2.0, 2.0)
plt.ylim(-4.0, 4.0)
plt.show()

3次近似

 1次近似よりも少し広い範囲で近似できているのが分かります。

 4次近似を計算します。

$$ f(x) \approx f(0) + f'(0) x + \frac{1}{2!} f''(0) x^2 + \frac{1}{3!} f'''(0) x^3 + \frac{1}{4!} f''''(0) x^4 $$

 $\sin(x)$の4階微分は$\sin(x)$です。$\sin(0) = 0$なので

$$ \begin{aligned} \sin(x) &\approx \sin(0) + \cos(0) x - \frac{1}{2!} \sin(0) x^2 - \frac{1}{3!} \cos(0) x^3 + \frac{1}{4!} \sin(0) x^4 \\ &= 0 + x - 0 - \frac{1}{3!} x^3 + 0 \\ &= x - \frac{x^3}{3!} \end{aligned} $$

となります。これも3次近似と同じ結果になりました。

 5次近似を計算します。

$$ f(x) \approx f(0) + f'(0) x + \frac{1}{2!} f''(0) x^2 + \frac{1}{3!} f'''(0) x^3 + \frac{1}{4!} f''''(0) x^4 + \frac{1}{5!} f'''''(0) x^5 $$

 $\sin(x)$の5階微分は$\cos(x)$です。$\cos(0) = 1$なので

$$ \begin{aligned} \sin(x) &\approx \sin(0) + \cos(0) x - \frac{1}{2!} \sin(0) x^2 - \frac{1}{3!} \cos(0) x^3 + \frac{1}{4!} \sin(0) x^4 + \frac{1}{5!} \cos(0) x^5 \\ &= 0 + x - 0 - \frac{1}{3!} x^3 + 0 + \frac{1}{5!} x^5 \\ &= x - \frac{x^3}{3!} + \frac{x^5}{5!} \end{aligned} $$

となります。

# 5次近似を計算
approx_x5 = approx_x3 + 1.0 / math.factorial(5) * np.cos(a) * (x - a)**5
#approx_x5 = approx_x3 + 1.0 / math.factorial(5) * np.cos(a) * x**5
#approx_x5 = approx_x3 + x**5 / math.factorial(5)
# 作図
plt.figure(figsize=(10, 5))
plt.plot(x, sin_x, label='sin(x)') # 元の関数
plt.plot(x, approx_x1, linestyle='dotted', label='1次近似') # 1次近似
plt.plot(x, approx_x3, linestyle='dotted', label='3次近似') # 3次近似
plt.plot(x, approx_x5, label='5次近似') # 5次近似
plt.grid()
plt.legend(prop={"family":"MS Gothic"})
#plt.xlim(-2.5, 2.5)
plt.ylim(-4.0, 4.0)
plt.show()

5次近似

 3次近似のときよりも更に広い範囲で近似できました。

・sin関数のマクローリン展開の実装

 ここまでで、$\sin (x)$のマクローリン展開に関して次のパターンが見えたと思います。

$$ \sin(x) = \frac{x^1}{1!} - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \frac{x^9}{9!} - \cdots $$

 この計算は、総和の記号$\sum$を使って次の式で表記できます。

$$ \sin(x) = \frac{x}{1!} - \frac{x^3}{3!} + \frac{x^5}{5!} - \cdots = \sum_{i=1}^{\infty} (-1)^i \frac{x^{2i+1}}{(2 i + 1)!} \tag{27.3} $$

 この式が$\sin$関数のマクローリン展開です。$(-1)^i$は、$i$が偶数のとき$+1$になり、奇数のとき$-1$になります。また、$2i + 1$は常に奇数になります。
 近似に用いる要素数が多いほど近似精度が良くなり、無限個のとき等しくなります。

 $\sin (x)$のマクローリン展開による$n$次近似(無限個ではなく$n$個の要素による近似)は次の式になります。

$$ \sin(x) \approx \sum_{i=1}^n (-1)^i \frac{x^{2i+1}}{(2 i + 1)!} $$

 この計算を関数として定義します。

# マクローリン展開によるsin(x)の近似式を定義
def sin_maclaurin(x, n):
    # 初期化
    y = 0
    
    # n次近似を計算
    for i in range(n):
        c = (-1)**i / math.factorial(2 * i + 1)
        t = c * x**(2 * i + 1)
        y = y + t
    return y


 作成した関数を使って、先ほどまでの計算を行ってみましょう。

# 近似する次元数を指定
n = 9

# x軸の値を設定
x = np.arange(-10.0, 10.0, 0.01)

# sin(x)を計算
sin_y = np.sin(x)

# n次近似を計算
approx_xn = sin_maclaurin(x, n)

# 作図
plt.figure(figsize=(10, 5)) # 画像サイズ
plt.plot(x, sin_y, label='sin(x)') # 元の関数
plt.plot(x, approx_xn, label=str(n) + '次近似') # n次近似
plt.grid() # グリッド線
plt.legend(prop={"family":"MS Gothic"}) # 凡例
plt.ylim(-4.0, 4.0) # y軸の表示範囲
plt.show()

マクローリン展開による$\sin$の9次近似


 $n$が大きいほど近似できる範囲が大きくなる(近似精度がよくなる)のでした。それをアニメーションで確認します。

・コード(クリックで展開)

# 近似する最大次元数を指定
n_max = 20

# x軸の値を設定
x = np.arange(-15.0, 15.0, 0.01)

# y軸の値を計算
sin_x = np.sin(x)

# 画像サイズを指定
fig = plt.figure(figsize=(10, 5))

# 作図処理を関数として定義
def update(n):
    # n次近似を計算
    approx_xn = sin_maclaurin(x, n + 1)
    
    # 前フレームのグラフを初期化
    plt.cla()
    
    # 作図
    plt.plot(x, sin_x, label='sin(x)') # 元の関数
    plt.plot(x, approx_xn, label=str(n + 1) + '次近似') # n次近似
    plt.grid() # グリッド線
    plt.legend(prop={"family":"MS Gothic"}, loc='upper right') # 凡例
    plt.ylim(-4.0, 4.0) # y軸の表示範囲

# gif画像を作成
approx_anime = animation.FuncAnimation(fig, update, frames=n_max, interval=200)

# gif画像を保存
approx_anime.save("step27_sin_maclaurin.gif")

マクローリン展開による$\sin$の近似

 $n$が大きくなるにしたがって近似できていくのを確認できます。

 $n$を無限に大きくすると等しくなるとのことでした。しかし、無限には計算できませんし、無限でなくとも処理は少ないほど嬉しいものです。そこで次は、何次近似まで求めるかを考えます。

・どこまで計算するか

 $\frac{x^{2i+1}}{(2 i + 1)!}$について掘り下げてみましょう。ただしここでは簡単に、$\frac{x^i}{i!}$として扱います。

 $x = 5$において、$i = 3$と$i = 9$のときを計算してみます。

$$ \begin{aligned} \frac{5^1}{1!} &= 5 \\ \frac{5^3}{3!} &= \frac{ 5 * 5 * 5 }{ 3 * 2 * 1 } = 20.833 \\ \frac{5^9}{9!} &= \frac{ 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 }{ 9 * 8 * 7 * 6 * 5 * 4 * 3 * 2 * 1 } = 5.382 \end{aligned} $$
# 値を指定
x = 5
i = 9

# i番目の項の計算
y = x**i / math.factorial(i)
print(y)
5.3822889109347445

 $i$の値によって増減しそうです。

 $i = 0, 1, 2, \cdots, n$の範囲で計算して、グラフにしましょう。math.factorial()はスカラしか処理できないので、リスト内包表記で処理します。

# 値を指定
x = 5
n = 20

# y軸の値を計算
y = np.array([x**i / math.factorial(i) for i in range(n + 1)])

# 作図
plt.figure(figsize=(10, 5)) # 画像サイズ
plt.plot(np.arange(n + 1), y, label='x=' + str(x))
plt.xlabel('i') 
plt.ylabel('y')
plt.title('$y = \\frac{x^i}{i!}$', fontsize=20)
plt.grid() # グリッド線
plt.legend() # 凡例
plt.show()

$i$番目の項の推移

 x軸は$i$です。$x > i$の間は値が増加しますが、$x < i$だと減少していくのが分かります。

 $x$の幅を持たせた上で、$i$の変化による$y$の変化をアニメーションで確認してみましょう。

・コード(クリックで展開)

# x軸の値を設定
x = np.arange(-15.0, 15.0, 0.01)

# 次数の最大値を指定
i_max = 50

# 画像サイズを指定
fig = plt.figure(figsize=(10, 5))

# 作図処理を関数として定義
def update(i):
    # i番目の項を計算
    y = x**i / math.factorial(i)
    
    # 前フレームのグラフを初期化
    plt.cla()
    
    # 作図
    plt.plot(x, y, label='i=' + str(i))
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('$y = \\frac{x^i}{i!}$', fontsize=20)
    plt.grid() # グリッド線
    plt.legend(loc='upper right') # 凡例
    plt.ylim(-10.0, 10.0) # y軸の表示範囲

# gif画像を作成
approx_anime = animation.FuncAnimation(fig, update, frames=i_max + 1, interval=200)
approx_anime.save("step27_factorial.gif")

$i$番目の項の推移

 $i$が大きくなるにしたがって、$y$がほとんど0になり、またその影響を受ける$x$の範囲が広がるのが分かります。$y$がほとんど0とは、近似にほとんど影響しないということです。つまり、必要な$x$の範囲で$y$の影響が小さくなると処理を打ち切って(計算を止めて)もよさそうです。

 本では、tの値が指定した値(閾値)threshodlを下回ると計算を止めるように実装しています。

・$\sin$関数のテイラー展開の実装

 最後に、$\sin (x)$のテイラー展開を行います。

 テイラー展開の式(27.1)も$\sum$を使って次の式で表記できます。

$$ f(x) = \frac{1}{0!} f(a) (x - a)^0 + \frac{1}{1!} f'(a) (x - a)^1 + \frac{1}{2!} f''(a) (x - a)^2 + \cdots = \sum_{i=1}^{\infty} \frac{f^{(i)}(a)}{i!} (x - a)^i $$

 $i$階微分($f'(x)$のダッシュの数が$i$個)を$f^{(i)}(x)$で表すことにします($f'''(x)$なら$f^{(3)}(x)$です)。また、0の階乗$0!$と0乗$x^0$は1と定義されています。$0! = 1,\ x^0 = 1$なので、式(2.71)の1つ目と2つ目の項を変形すると右辺が成り立つのが分かりやすと思います。
 無限個の要素を足し合わせると元の関数と等しくなります。

 $\sin (x)$のテイラー展開による$n$次近似($n$個の要素による近似)は次の式になります。

$$ \sin (x) \approx \sum_{i=1}^n \frac{\{\sin(a)\}^{(i)}}{i!} (x - a)^i $$

 また、$\sin (x)$の$i$階微分は次になります。

$$ \{\sin (x)\}^{(i)} = \begin{cases} \cos (x) &\quad i = 1, 5, 9, \cdots \\ - \sin (x) &\quad i = 2, 6, 10, \cdots \\ - \cos (x) &\quad i = 3, 7, 11, \cdots \\ \sin (x) &\quad i = 4, 8, 12, \cdots \end{cases} $$

 $i$を4で割ったときに、余りが1なら$\cos (x)$、2なら$- \sin (x)$、3なら$- \cos (x)$、0なら$\sin (x)$です。

 この計算を関数として定義します。割り算の余り(剰余)は%で計算できます。

# テイラー展開によるsin(x)の近似式を定義
def sin_taylor(x, a, n):
    # 初期化
    y = 0
    
    # n次近似を計算
    for i in range(n):
        if i % 4 == 1:
            c = np.cos(a) / math.factorial(i)
        elif i % 4 == 2:
            c = - np.sin(a) / math.factorial(i)
        elif i % 4 == 3:
            c = - np.cos(a) / math.factorial(i)
        elif i % 4 == 0:
            c = np.sin(a) / math.factorial(i)
        t = c * (x - a)**i
        y = y + t
    return y

 細かい違いですが、sin_maclaurin()では$i = 0$の計算は行われませんでしたが、sin_taylor()では行います。

 早速計算してみましょう。

・コード(クリックで展開)

 マクローリン展開の$n$次近似の処理コードと、オブジェクトaを作成してsin_maclaurin()sin_taylor()に置き換えたところだけが異なります。

# 値を指定
a = 5.0
n = 20

# x軸の値を設定
x = np.arange(-15.0, 15.0, 0.01)

# sin(x)を計算
sin_y = np.sin(x)

# n次近似を計算
approx_xn = sin_taylor(x, a, n)

# 作図
plt.figure(figsize=(10, 5)) # 画像サイズ
plt.plot(x, sin_y, label='sin(x)') # 元の関数
plt.plot(x, approx_xn, label=str(n) + '次近似') # n次近似
plt.grid() # グリッド線
plt.legend(prop={"family":"MS Gothic"}) # 凡例
plt.ylim(-4.0, 4.0) # y軸の表示範囲
plt.show()

テイラー展開による$\sin$の20次近似

 $x = 5$を中心に近似されていきます。

・コード(クリックで展開)

 アニメーションのコードについても、aを作成してsin_taylor()を使うところだけが異なります。

# 値を指定
n_max = 50
a = 5

# x軸の値を設定
x = np.arange(-15.0, 15.0, 0.01)

# y軸の値を計算
sin_x = np.sin(x)

# 画像サイズを指定
fig = plt.figure(figsize=(10, 5))

# 作図処理を関数として定義
def update(n):
    # n次近似を計算
    approx_xn = sin_taylor(x, a, n + 1)
    
    # 前フレームのグラフを初期化
    plt.cla()
    
    # 作図
    plt.plot(x, sin_x, label='sin(x)') # 元の関数
    plt.plot(x, approx_xn, label=str(n + 1) + '次近似') # n次近似
    plt.grid() # グリッド線
    plt.legend(prop={"family":"MS Gothic"}, loc='upper right') # 凡例
    plt.ylim(-4.0, 4.0) # y軸の表示範囲

# gif画像を作成
approx_anime = animation.FuncAnimation(fig, update, frames=n_max, interval=200)

# gif画像を保存
approx_anime.save("step27_sin_taylor.gif")

テイラー展開による$\sin$の近似


 以上で、$\sin$関数を用いてテイラー展開について確認しました。ステップ29にて、テイラー展開を用いて最適化を行います。

参考文献

おわりに

 ゼロつく3巻ではこの記事が1つ目の記事です。3巻での理論的な内容は、1・2巻の知識(記事)で何とかなるかなぁ、攻略ノートを作らなくてもいいかなぁと思いながら進めていました。が、テイラー展開については他の本で登場したときにもふわっとした理解で読み飛ばしたので、ここらでしっかり理解しておこうというのもあり、丁寧にやってみることにしました。折角書いたので誰かの役に立てばとブログにて共有しておきます。
 しかしこうも寄り道してると本編が全然進みません。この先の内容も必要に応じて書いたり書かなかったりするつもりです。

 しかしホント数学ネタは上げたくないです。いや他のネタなら気にならないかと言ったらそんなこともないですが。だからこそしっかり調べて理解しようという動機になるのですが。今回は間をとって理論的な内容を書くのは避けつつ解説を試みました。なぜ式(27.1)で近似できるのかという話も面白いので、時間があれば調べてみるといいですよ。本編の進捗に影響しますけどね。

 この記事の投稿4時間前にアップされたカバー動画です。勉強のお供にぜひ♪

【次ステップの内容】

www.anarchive-beta.com