最急降下法とマクローリン展開を用いた円周率の推定手法について

はじめに

円周率の近似値を推定する方法は、多数知られている。具体的には、級数を用いた解析学的な手法や、幾何学を用いたオーソドックスな方法が挙げられる。今回は、$cos \theta$は$0<\theta <2 \pi$の区間内では、$\theta=\pi$のとき最小値$-1$を取ることをもとにして円周率の推定を行う。ただし、目的関数の極小値を探索するために、再急降下法というアルゴリズムを用いた。

以下の画像は、初期値を0.5とした探索結果である。このように、試行回数を増やすほど、円周率$3.14 \cdot \cdot \cdot$に漸近していくことが分かる。

最急降下法による最小値探索

導入

以下の余弦関数を、$0<\theta<2\pi$の範囲で探索し、最小値-1を取る$\theta$ の値つまり、$\pi$の値を推定する。つまり、

cos \theta = -1 (0<\theta<2 \pi)

の解を探索する。

余弦のマクローリン展開について

余弦のマクローリン展開は以下のサイトにより、次のように表すことができる。

cos \theta = \sum_{n=0}^{\infty}{(-1)^n \frac{\theta^{2n}}{(2n)!}}

今回は、10次の多項式として余弦を近似する。

最急降下法

簡潔に最急降下法アルゴリズムを説明しておこう。$y=f(x)$の最小値を求めることを考える。
解空間の最小値を取る極値が一つ解空間の谷が一つの場合に有効な最小値の推定手法である。$a$に初期値を与え、以下のような式を用いて$a$を更新していく。つまり、$a_k$を算出する。

a_{k+1}=a_{k}-\alpha(\frac{\delta }{\delta a}f(a_{k}))

ただし$\alpha,\beta$は$(a,b)$の学習率である。
これは、ボールが坂に下り落ちる様に似ている。偏微分係数は各方向の坂の傾きと捉えることができる。

プログラム

上記の考えをPythonプログラムによって反映させると以下の通りになる。

python talor.py

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
#再急降下法の学習率
alpha=0.01
#微小分
delta=1.0*10**-5
#マクローリン展開において近似する多項式の次数
n=10
#マクローリン展開
def cos_near(x):
  num_ary=[]
  val_ary=[]
  #   # valは(2n)の階乗である
  for m in range(2*n):
    val=1
    num=2*m
    for i in range(1,num+1):
      val=val*i
    num_ary.append(num)   
    val_ary.append(val)
  B=0
  for i in range(len(num_ary)):
    A=(-1)**i*(x**num_ary[i]/val_ary[i])
    B=B+A  
  return B
#試行回数
i=0
#各値をしまう配列
x_ary=[]
y_ary=[]
i_ary=[]
err_ary=[]
#初期値
x=0.5
#再急降下法のメインループ
while cos_near(x)>-1+delta:
    #微分の計算
    dfdx=(cos_near(x+delta)-cos_near(x))/delta
    #再急降下法によるパラメータの更新
    x=x-alpha*dfdx
    y=cos_near(x)
    err=(abs(x-math.pi)/math.pi)*100
    #反復回数と目的関数の値を格納
    i_ary.append(i)
    x_ary.append(x)
    y_ary.append(y)
    err_ary.append(err)
    i=i+1
##結果の出力
plt.xlabel('反復回数')
plt.ylabel('円周率の近似値')
plt.plot(i_ary,x_ary)
plt.savefig('円周率の近似値を最急降下法を用いて求めるアルゴリズム.png')
plt.show()
##結果の出力
plt.xlabel('反復回数')
plt.ylabel('円周率の近似値の誤差[%]')
plt.plot(i_ary,err_ary)
plt.savefig('円周率の近似値を最急降下法を用いて求めるアルゴリズム_誤差.png')
plt.show()

xx=np.linspace(0,2*math.pi,100)
yy=cos_near(xx)
##結果の出力
plt.plot(xx,yy)
plt.scatter(x_ary,y_ary)
plt.savefig('円周率の近似値を最急降下法を用いて求めるアルゴリズム_余弦波.png')
plt.show()

結果

試行回数による円周率の推定値の推移と誤差の推移、余弦カーブ上における推定点の推移のグラフを示す。

試行回数による円周率の推定値の推移と誤差の推移

以下の画像のように、試行回数を増やすほど、円周率の真値に漸近されていることが分かる。

誤差の推移

誤差は、試行回数を増やすにつれて減少傾向にあるが、特に中盤での減少度が激しかった。

余弦カーブ上における推定点の推移

余弦カーブ上の点の推移は、以下のようになった。

このように、最小値を探索することができていることが分かる。

まとめ

今回は、三角関数のマクローリン展開と最急降下法を用いて、円周率を推定してみた。結果、かなり良い精度で円周率を推定することができた。円周率の推定方法は大学入試の積分問題や東大入試でもだされたこともあり非常に奥が深い。なので、今回のアルゴリズムをしっかりと学習することで、基礎力の強化を図りたい。