はじめに
テイラー展開は、複雑な関数を累乗の多項式関数の無限和で表現するという技法である。したがって、元の関数から近似関数を引いたものを2乗したものを積分した誤差関数は極めて0に小さくなるはずである。そこで、今回はこのような誤差関数が最小になる係数を最急降下法で推定し、係数がテイラー展開の理論値と比較しあまりズレていないということを確認する。ただし、題材の関数として$exp(-x)$を用いた。
誤差関数の低下の推移
最終的な係数の比較
問題設定
I_n=\int_0^{n}\{ \exp(-x)-\sum_{k=0}^{n} a_k x^k \}^2\exp(-x) dx
と表されるとき、$I_{\infty}$を最小にするのは、
a_k=\frac{(-1)^k}{k!}
となるはず。これを示す。
ただし、$1\le k\le n$である。
解答
まず、2乗項を計算する。
(\exp(-x)-\sum_{k=0}^{n}a_k x^k)^2=\exp(-2x)-2\exp(-x)\sum_{k=0}^n a_k x^k+(\sum_{k=0}^n a_k x^k)^2
ここで、
(\sum_{k=0}^n a_k x^k)^2=\sum_{k=0}^n \sum_{m=0}^n a_k a_m x^{k+m}
I_n=\int_0^{n}\{ \exp(-3x)+\sum_{k=0}^n \sum_{m=0}^n a_k a_m x^{k+m}\exp(-x)-2\sum_{k=0}^{n}{ a_k x^k \exp(-2x) }\} dx
したがって、積分を前に出すと以下のとおりである。
I_n=\int_0^{n} \exp(-3x) dx+\sum_{k=0}^n \sum_{m=0}^n a_k a_m \int_0^{n}x^{k+m}\exp(-x)dx-2\sum_{k=0}^{n}a_k\int_0^{n}x^{k}\exp(-2x)dx
ゆえに、
I_n=\frac{1}{3}+\sum_{k=0}^{n}(\sum_{m=0}^nc_k a_k a_m - 2b_k a_k )
と表すことができる。
ただし、
b_k= \int_0^{n}x^{k}\exp(-2x)dx=\frac{1}{2^{k+1}}\int_0^{n}(2x)^{k}\exp(-2x)dx
ここで、$n\to \infty$のとき、ガンマ関数の性質より、$b_k \to \frac{k!}{2^{k+1}}$となる。
一方で、$c_k$は以下のように表されるの。
c_k= \int_0^{n}x^{k+m}\exp(-x)dx\to (k+m)!
ここで、$I_n$を$a_k$で偏微分することを考える。
\frac{\delta I_n}{\delta a_k}=-2b_k+2(\sum_{m=0}^n a_k c_k)
これが0に等しいとき$I_nは$最小値を取る。
したがって、
\frac{\delta I_n}{\delta a_k}=-\frac{k!}{2^k}+2(\sum_{m=0}^n a_k (k+m)!) =0
ここで、上記の連立方程式の解の1つが、
a_k=\frac{(-1)^k}{k!}
であると仮定する。$n\to \infty$であるとき、
-\frac{k!}{2^k}+2(\sum_{m=0}^n a_k (k+m)!) =-\frac{k!}{2^k}+2(\sum_{m=0}^\infty (\frac{(-1)^k}{k!}) (k+m)!) =-\frac{k!}{2^k}+2(k! \sum_{m=0}^{\infty}{_{k+m}C_{k}})=-\frac{k!}{2^k}+\frac{2 k!}{2^(k+1)}=0
ここで、
(\sum_{m=0}^\infty (\frac{(-1)^k}{k!}) (k+m)!) =k! \sum_{m=0}^\infty{_{k+m}C_{k}}=k! \sum_{m=0}^n{_{k+m}C_{m}}=k! \{ 1-(-1)\}^{-(n+1)}=k! 2^{-(k+1)}
という負の二項定理なるものを用いた。
したがって、
a_k=\frac{(-1)^k}{k!}
である。
プログラム
上記の議論をプログラムに反映すると、以下のようになる。
ただし、Pythonによる数値積分（シンプソン法）は以下の記事を参考にした。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
#a :下限
a=0
#b :上限
b=1
#m :分割数
m=100
n = int(m/2)
#N次式を作成
N=10
#係数の初期値を0として作成する。
a_ary=np.zeros(N+1)
# 積分する関数の定義
def f(x,a_ary):
#f_nearは多項式関数
f_near=0
for i in range(len(a_ary)):
f_near=f_near+(a_ary[i]*x**(i))
err=(np.exp(-x)-f_near)**2*np.exp(-x)
return err
def Simpson(f,a,b,m,a_ary):
df = (b-a)/(2*n)
h = (b-a)/m
sum = 0
for i in range(1,n+1):
x1 = a + df*(2*i -2)
x2 = a + df*(2*i -1)
x3 = a + df*(2*i)
f1 = f(x1, a_ary)
f2 = f(x2, a_ary)
f3 = f(x3, a_ary)
sum += h*(f1 + 4*f2 + f3)/3
return sum
hh=1.0*10**(-5)
# 積算回数
iterations=100
# 学習率
alpha=0.1
# 損失関数
error_pre=[]
# 面積分が0.001以下になったときに計算を終了する。
while Simpson(f,a,b,m,a_ary) > 0.001:
for p in range(0,len(a_ary)):
a_ary2=np.copy(a_ary)
a_ary2[p]=a_ary[p]+hh
d_err_a=(Simpson(f,a,b,m,a_ary2)-Simpson(f,a,b,m,a_ary))/hh
a_ary[p]=a_ary[p]-alpha*d_err_a
error_pre.append(Simpson(f,a,b,m,a_ary))
# 結果の表示
plt.plot(error_pre)
plt.xlabel("Iterations")
plt.ylabel("Error")
plt.title("Error Reduction")
plt.grid()
plt.savefig("error_reduction_talor.png")
plt.show()
#階乗計算
def factorial(num):
num2=1
for i in range(1,num+1):
num2=num2*i
return num2
y_ary=[]
xx = np.arange(len(a_ary))
width = 0.35
for i in range(len(a_ary)):
y_ary.append((-1)**i/factorial(i))
#計算値
plt.bar(xx, a_ary, width,color='blue',label="計算値")
#理論値
plt.plot(xx, y_ary, color='red', label='(-1)^n/n!（理論値）')
plt.xlabel("次数")
plt.ylabel("係数")
plt.title("係数のヒストグラム")
plt.legend()
plt.savefig("histogram_a_talor.png")
plt.show()
これを実行すると以下のようなグラフが出力される。
誤差関数の低下の推移
最終的な係数の比較
青い棒グラフは、最終的な試行において生成された係数の推定値である。
このように、初期値を0としているため、精度はそこまで高くはないが正負はほぼ一致している。
まとめ
今回は、$\exp(-x)$をテイラー展開したときの各整数次における係数を推定するプログラムを作成した。結果、やはり損失関数を最小にするためには、テイラー展開で算出できる係数を使う必要性があるということがわかった。
また、理論計算で抽象度が上がってしまったため、GeminiやChatGptなどのAIを用いて試行錯誤した。