最急降下法によるフーリエ級数展開の数値的推定

はじめに

以前の記事[1]では、テイラー展開、すなわち関数を累乗の多項式の有限和で近似する問題を、係数を未知パラメータとした二乗誤差の最小化問題として定式化し、最急降下法によって数値的に係数を求めた。

フーリエ級数展開も、近似対象を三角関数の有限和で表現するという点でテイラー展開と本質的に同じ構造を持つ。そこで今回は、有限区間 $(0,\ \pi/4)$ 上で定義された関数 $y=x$ を対象に、正弦関数の有限和

f_N(x)=\sum_{n=0}^{N} a_n \sin\!\left(\frac{n\pi x}{L}\right),\qquad L=\frac{\pi}{4}

への展開係数 $a_0,\dots,a_N$ を、二乗誤差を目的関数とした最急降下法により数値的に求めた。

理論値が解析的に導出できる問題であるため、数値解と理論解を直接比較できる点が本稿の主眼である。

問題設定

区間 $(0,L)$ 上で定義された関数 $y=x$ に対し、次数 $N$ の正弦級数 $f_N(x)$ による近似を考える。近似の良さを測る目的関数として、区間上の二乗誤差積分

I(a_0,\dots,a_N)=\int_0^{L}\left\{x-\sum_{n=0}^{N} a_n \sin\!\left(\frac{n\pi x}{L}\right)\right\}^2 dx

を導入し、これを最小化する係数ベクトル $(a_0,\dots,a_N)$ を求める問題として定式化する。なお $n=0$ の項は $\sin(0)=0$ であるから、係数 $a_0$ は目的関数に一切寄与しない恒等的に無効な変数である。

理論的な解

区間 $(0,L)$ における関数系 ${\sin(n\pi x/L)}$ は、次の直交性

\int_0^L \sin\!\left(\frac{m\pi x}{L}\right)\sin\!\left(\frac{n\pi x}{L}\right) dx = \frac{L}{2}\delta_{mn}

を満たす。

この直交性より、$I$ を最小化する係数は、通常のフーリエ正弦係数の公式

a_n = \frac{2}{L}\int_0^L x\sin\!\left(\frac{n\pi x}{L}\right) dx

で解析的に与えられる。部分積分によりこれを計算すると

a_n = \frac{2L(-1)^{n+1}}{n\pi}\qquad (n\ge 1)

を得る。$L=\pi/4$ の場合、これは $a_n=(-1)^{n+1}/(2n)$ に一致する。今回はこの理論値を、最急降下法による数値解の妥当性を確認するための基準として用いる。

数値解法(最急降下法)

$I(a_0,\dots,a_N)$ を各係数 $a_p$ について中心差分ではなく前進差分で数値微分し、

\frac{\partial I}{\partial a_p}\approx \frac{I(a_0,\dots,a_p+h,\dots,a_N)-I(a_0,\dots,a_p,\dots,a_N)}{h}

を得たうえで、学習率 $\alpha$ による更新則

a_p \leftarrow a_p-\alpha\frac{\partial I}{\partial a_p}

をすべての係数に対して反復適用した。積分 $I$ 自体は解析的に評価せず、シンプソン則による数値積分で評価している。これにより、目的関数の形を陽に微分せずとも、任意の関数形に対して同じ枠組みで最急降下法を適用できる。

プログラム

実装したPythonコードは以下の通りである。

furie.py

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math

#a :下限
a=0
#b :上限
b=math.pi/4
#L :区間の長さ(直交基底 sin(n*pi*x/L) の周期を決める)
L=b-a
#m :分割数
m=1000
n = int(m/2)

#N次式を作成
N=10
#係数の初期値をランダムに作成する。
a_ary= np.random.randn( N+1)

# 積分する関数の定義
def f(x,a_ary):
  #f_nearはsin関数の和
  f_near=0
  for i in range(len(a_ary)):
     f_near=f_near+(a_ary[i]*np.sin(i*np.pi*x/L))
  err=(x-f_near)**2
  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=1000
# 学習率
alpha=0.01
# 損失関数
error_pre=[]
# 面積分が0.01以下、またはiterations回数に達したときに計算を終了する。
it=0
while Simpson(f,a,b,m,a_ary) > 0.01 and it < iterations:
  it+=1
  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))
  print(Simpson(f,a,b,m,a_ary))

# 結果の表示
plt.figure()
plt.plot(error_pre)
plt.xlabel("Iterations")
plt.ylabel("Error")
plt.title("Error Reduction")
plt.grid()
plt.savefig("error_reduction_furie.png")
plt.show()

y_ary=[]
xx = np.arange(len(a_ary))
width = 0.35

#理論値: (0,L)上でy=xを正弦級数展開したときの係数 b_n = 2L(-1)^(n+1)/(n*pi)
y_ary.append(0)
for i in range(1,len(a_ary)):
  y_ary.append(2*L*(-1)**(i+1)/(i*math.pi))

plt.figure()
#計算値
plt.bar(xx, a_ary, width,color='blue',label="計算値")
#理論値
plt.plot(xx, y_ary, color='red', label='2L(-1)^(n+1)/(n*pi)（理論値）')
plt.xlabel("次数")
plt.ylabel("係数")
plt.title("係数のヒストグラム")
plt.legend()
plt.savefig("histogram_a_furie.png")
plt.show()

num=1000
xx_ary=np.linspace(a,b,num)
yy_ary=np.zeros(num)

for k in range(num):
  for i in range(len(a_ary)):
    yy_ary[k]+=(a_ary[i]*np.sin(i*np.pi*xx_ary[k]/L))

plt.figure()
plt.plot(xx_ary,xx_ary,color='red',label='y=x（理論値）')
plt.plot(xx_ary,yy_ary,color='blue',label='フーリエ級数近似（計算値）')
plt.xlabel("x")
plt.ylabel("y")
plt.title("フーリエ級数によるy=xの近似")
plt.legend()
plt.savefig("fit_curve_furie.png")
plt.show()

区間長 $L$ に対して基底関数の角周波数を $n\pi/L$ に取っている点が要点である。基底の周波数を区間長に合わせて正しく取らなければ、基底同士が直交せず最急降下法が正しく収束しない。

結果

グラフ1:誤差の収束過程

反復回数に対する二乗誤差積分 $I$ の推移。

反復開始時は係数の初期値をランダムに与えたため誤差は $I\approx3.6$ 程度であったが、反復を重ねるごとに単調に減少し、およそ550回の反復で終了条件である $I<0.01$ に達した。

グラフ2:係数の理論値と計算値の比較

青棒が最急降下法による計算値、赤線が理論値 $a_n=2L(-1)^{n+1}/(n\pi)$。

$n=1$ 以上の係数については、計算値と理論値がほぼ完全に一致した。$n=0$ の棒だけは理論値(定義上0)から大きく外れているが、これは $\sin(0)=0$ より $a_0$ が目的関数に一切影響しない変数であり、初期のランダム値のまま更新されずに残っているためで、近似の精度とは無関係である。

グラフ3:近似曲線と理論直線の比較

赤線が $y=x$、青線が求まった係数によるフーリエ級数近似曲線。

区間の中央付近では両者はよく一致しているが、右端 $x\to L=\pi/4$ に近づくにつれて近似曲線が理論直線から外れ、振動しながら大きくオーバーシュートしている様子が確認できる。

考察

実際にプログラムを動かしたところ、係数そのものは理論値とほぼ一致したにもかかわらず、グラフ3に示したように近似曲線は区間の右端付近で理論値から明らかに外れる部分があった。これは有限次数 $N$ の三角関数の和でしか関数を表現できないことに起因する誤差であり、より具体的には次の理由による。

正弦級数 $\sum a_n\sin(n\pi x/L)$ は、$x=0$ と $x=L$ で必ず値0を取る。ところが近似対象である $y=x$ は $x=0$ では $0$ に一致するものの、$x=L=\pi/4$ では $\pi/4\neq 0$ である。すなわち、この正弦級数が表現しているのは $y=x$ を奇関数として周期 $2L$ で拡張した関数であり、その拡張は $x=L$ の点で $\pi/4$ から $-\pi/4$ へ飛び移る不連続なジャンプを持つ。有限次数の三角関数の和で不連続点(あるいはその近傍)を表現しようとすると、次数をいくら上げても一定の割合のオーバーシュートが解消されないギブス現象が生じることが知られており、グラフ3で見られた右端付近の振動的なずれは、まさにこのギブス現象として説明できる。

したがって、今回観測された「理論値からの外れ」は、最急降下法の収束不足や数値積分の誤差によるものではなく、有限個の三角関数の和という表現能力そのものの限界、より正確には近似対象の周期拡張が不連続点を持つことに起因する本質的な誤差であるといえる。この誤差は次数 $N$ を増やしても $x=L$ 近傍では大きさが漸近的に一定値に収束し完全には解消されない一方、不連続点から離れた区間内部では次数を増やすほど改善されると考えられる。

まとめ

区間 $(0,\pi/4)$ 上で定義された $y=x$ を対象に、正弦関数の有限和による近似係数を、二乗誤差を目的関数とした最急降下法で数値的に求めた。基底関数の角周波数を区間長 $L$ に対応する $n\pi/L$ に正しく取ることで、数値解は理論値である $a_n=2L(-1)^{n+1}/(n\pi)$ とほぼ完全に一致した。一方で、近似曲線は区間端 $x=L$ の近傍でギブス現象による理論値からのずれが残ることを確認した。これは最急降下法という数値解法固有の誤差ではなく、有限次数の三角関数の和で不連続な周期拡張を表現しようとすることに起因する本質的な誤差であると考えられる。

参考文献

1. テイラー展開の係数を最急降下法で求める https://qiita.com/arairuca/items/32d9d286fb46d0c7c0b2