微分の近似について

flickrholdr_550_300_math
  • ????????????????????

すこしマニアックな話題を。ある関数の微分をコンピュータ上で計算するための近似式として最も基本的なものは以下のものでしょう。差分(商)とか呼ばれるものですね。
\begin{align*}
f'(x)\simeq\frac{f(x+\Delta x)-f(x)}{\Delta x}
\end{align*}
もうすこし精度をよくしたい場合に次のような式もよく利用されます。
\begin{align*}f'(x)\simeq\frac{f(x+\Delta x)-f(x-\Delta x)}{2\Delta x}\end{align*}
しかし複素数を使ったこんな近似の仕方もあるようです。
\begin{align*}f'(x)\simeq\frac{\text{Im}(f(x+i\Delta x))}{\Delta x}\end{align*}
直感的にはよくわからないこの式がうまくいくことを理論と数値実験でみてみよう、というのがこの記事の目的です。

差分商の裏側

差分商というのは微分の定義式
\begin{align*}f'(x)=\lim_{\Delta x\rightarrow 0}\frac{f(x+\Delta x)-f(x)}{\Delta x}\end{align*}
という式の極限をとらずに有限の小さな値でとめてしまって近似する、ということがすぐにわかります。したがって近似はうまくいきそうなものです。実際、これは理論的に確かめることができます。テイラー展開をしてみます。
\begin{align*}f(x+\Delta x)=f(x)+f'(x)\Delta x+\frac{f^{(2)}(x)}{2}\Delta x^2+\frac{f^{(3)}(x)}{2}\Delta x^3+\cdots\end{align*}
この式を変形して三角不等式を適用すると差分商の誤差を評価できます。
\begin{align*}\bigg\|\frac{f(x+\Delta x)-f(x)}{\Delta x}-f'(x)\bigg\|\le\bigg\|\frac{f^{(2)}(x)}{2}\Delta x\bigg\|+\cdots=O(\Delta x)\end{align*}
この式より差分商の刻み幅を小さくすると誤差がそれに比例して小さくなることがわかりました(これをオーダー1とかいいます)。
二番目の公式についても同様の計算をしてみると刻み幅の二乗に比例して誤差が小さくなる、つまりオーダー2であることがわかります(通常、刻み幅はゼロに近づけるものなので1乗よりは2乗、2乗よりは3乗に比例するほうが誤差の評価としてはよいということになります)。

問題の式の裏側

さて、問題の複素数の混じった公式についてです。ややこしそうですが、これも同様のロジックによって誤差の評価ができてしまいます。
まず、さきほどと同様ににテーラー展開を行います。しかし、今回加える摂動は純虚数の摂動です。
\begin{align*}f(x+i\Delta x)=\bigg[f(x)-\frac{f^{(2)}(x)}{2}\Delta x^2+\cdots\bigg]+i\bigg[f'(x)\Delta x-\frac{f^{(3)}(x)}{6}\Delta x^3+\cdots\bigg]\end{align*}
この式の両辺の虚部をとると誤差評価できます。
\begin{align*}\bigg\|\frac{\text{Im}(f(x+i\Delta x))}{\Delta x}-f'(x)\bigg\|\le\bigg\|\frac{f^{(3)}(x)}{6}\Delta x^2\bigg\|+\cdots=O(\Delta x^2)\end{align*}
複素関数についてあまり学んだことがない人はこんな疑問を持つ人がいるかもしれません。いま f という関数は実数を引数に取る関数なのに、なぜ複素数で定義できるのか?という疑問です。こんな疑問をもった人は「複素関数論」「解析接続」「解析関数」などのキーワードで検索してみてください。実際計算する段階では多くの言語では exp や sin などの関数は複素数にも対応した実装がされているので問題ないでしょう。

数値実験

簡単な数値実験をして見ましょう。
\begin{align*}f(x)=\cos(x)\qquad f'(x)=-\sin(x)\end{align*}
という関数の x=pi/4 での誤差を数値実験で評価してみましょう。これは以下のような簡単な R のスクリプトを使います(これは無視してかまいません)。

# Settings
x0 = pi/4.0                  # 微分する座標
f = function(x){ cos(x) }    # 微分される関数
df = function(x){ -sin(x) }  # 真の導関数

dx=cumprod(rep(1.003,10900))*1.0e-15  # 摂動の大きさをいろいろ変化させる

# Ordinal difference
df1=f(x0+dx)-f(x0)
err1=abs(df1/dx-df(x0))

# Lyness difference
df2=f(x0+(0+1i)*dx)
err2=abs(Im(df2)/dx-df(x0))

# Plot
plot(dx,err1,type="n",log="x",ylab="Error of d(cos)/dx at x=1")
lines(dx,err1,col="darkorange")
lines(dx,err2,col="darkblue")
legend("topright",legend=c("Ordinal","Lyness"),lw=2,col=c("darkorange","darkblue"))

# 対数プロット
plot(dx,err1,log="xy",type="l",ylim=c(1.0e-16,1e-1),ylab="log of error",col="darkorange")
lines(dx,err2,col="darkblue")
legend("bottomright",legend=c("Ordinal","Lyness"),lw=2,col=c("darkorange","darkblue"))

# 対数プロットの傾きを求める
thd=1e-6
lm(log(err1[dx>thd])~log(dx[dx>thd])) # 傾き 0.9982 : 理論通り!
lm(log(err2[dx>thd])~log(dx[dx>thd])) # 傾き 2.000  : 理論通り!

刻み幅をいろいろ変化させたときの誤差をプロットすると次のようになります。Lyness というほうが複素数のほうです。
すこしわかりにくいですが、Lyness のほうが普通の差分商(Ordinal)よりも(どの刻み幅においても)誤差が小さいですね。もうすこしみやすくするために縦軸を対数をとってみたグラフが次のグラフです。

かなりはっきりと違いがわかります。重要なのは2点

  • 刻み幅を減らしていくと、ともに誤差は減少していくが、差分商のほうが誤差の改善度合いが小さい
  • 刻み幅を減らしていくと、差分商のほうは途中から誤差が増加し始める

1点目はさきほど行った誤差のオーダーの議論からすぐにわかります。誤差オーダーの式の両辺の対数をとってみてください。対数グラフを書くとオーダーが直線の傾きとしてもとまります。傾きをフィッティングによって求めるとオーダーに一致する1、2という値が求まります(Rスクリプトの一番最後を見てみてください)。

2点目については桁落ち誤差が関係しています。差分商では同じような数の引き算が分母に現れますが、これが「桁落ち誤差」と呼ばれる誤差を発生させてしまいます(説明するのは大変なので気になる人は検索してみてください)。いっぽうの複素数の公式では引き算が存在しないため桁落ち誤差は存在しません。ということで誤差をマシンイプシロンと呼ばれるコンピュータにとっての特別な量の付近まで減らすことができるということです。

まとめ

複素数の関数は計算負荷が大きくなったり、場合によってはプログラミングが面倒になるかもしれません。また、誤差のオーダーについてはもっとよい公式はいくらでも存在します。しかし、ついつい忘れがちな桁落ち誤差を気にしなくて良いという点は◎です。

はてなブックマーク - 微分の近似について
Pocket