今週のホタル

今週も、こども自然公園へホタルを見に行きました。







６月になってから、気温が低くて、空気も乾燥しているせいか、思っていたよりホタルの数が増えていないし、元気もないみたいです。

---

話題は、変わって

最近、金融業界でも自動微分っていうのが、話題になっているようです。
海外では、２年ぐらい前から（もっと前から？）話題になっていたようで、昨年のRisk.netには、AAD vs GPUs: banks turn to maths trick as chips lose appealなんていう記事がでていたり、EUのHPCFinance2014では、CREDIT SUISSEが自動微分を使ったらモンテカルロの計算がとっても速くなりましたという発表REAL TIME COUNTERPARTY CREDIT RISK MANAGEMENT WITH ADJOINT ALGORITHMIC DIFFERENTIATION (AAD) をしているようです。

他にも、Googleの囲碁AI『AlphaGo』で使われているTensorflowや、ハミルトニアンモンテカルロ法のソフトウェア　Stan、PDEのFEniCS dolfin-adjointなど、さまざまなソフトで、自動微分が利用されているそうです。

自動微分、いろいろ方法があるそうですが、オペレータオーバーローディングができるプログラミング言語や型の制約が緩い言語では、二重数を使うと、わりと簡単に実装できるようです。

Pythonで簡単なプログラムを作ってみました。

例として、３次関数
f(x) = 2x^3 - 7x^2 + 8
を、

def func1(x):
  y = x.__class__(2.0) * x * x * x - x.__class__(7.0) * x * x + x.__class__(8.0)
  return y

こんなかんじで定義して
二重数

class Dual:
  def __init__(self, x, dx=None):
    self.x = x
    if dx is None:
      self.dx = 0.0
    else:
      self.dx = dx
  def __add__(self, other):
    return Dual(self.x + other.x, self.dx * 1 + other.dx * 1)
  def __sub__(self, other):
    return Dual(self.x - other.x, self.dx * 1 + other.dx * (-1))
  def __mul__(self, other):
    return Dual(self.x * other.x, self.dx * other.x + other.dx * self.x)
  def __div__(self, other):
    return Dual(self.x / other.x, self.dx * (1.0 / other.x) + other.dx * (-self.x / other.x / other.x))
  def var(self, x):
    self.x = x
    self.dx = 1.0

を使って、計算
>>> x = Dual(-2.0,1.0)
>>> y = func1(x)
>>> y.x
-36.0
>>> y.dx
52.0

できているっぽい、のでグラフにしてみます。

from numpy import *
import pylab as plt
x = linspace(-2.0, 4.0, 50)
y1 = []
y2 = []
for i in x:
  n = func1(Dual(i, 1.0))
  y1.append(n.x)
  y2.append(n.dx)

plt.title("f(x) = 2x^3 - 7x^2 + 8")
plt.plot(x, y1, color="b", label = "f(x)")
plt.plot(x, y2, color="r", label = "f'(x)")
plt.legend() 
plt.grid(color='g', linestyle='--')

plt.show()

closed-form だけじゃなく、ニュートン法を使った関数などでも、できるそうです。

平方根の例

def my_sqrt(x):
  y = x.__class__(1.0)
  d = x.__class__(2.0)
  for i in range(0,100):
    y = (y + x/y) / d
  return y

>>> x = Dual(2.0,1.0)
>>> y = my_sqrt(x)
>>> y.x
1.414213562373095
>>> y.dx
0.35355339059327373

これもグラフにしてみました。

x = linspace(1.0, 10.0, 50)
y1 = []
y2 = []
for i in x:
  n = my_sqrt(Dual(i, 1.0))
  y1.append(n.x)
  y2.append(n.dx)

plt.title('f(x) = sqrt(x)')
plt.plot(x, y1, color="b", label="f(x)")
plt.plot(x, y2, color="r", label="f'(x)")
plt.legend()
plt.grid(color='g', linestyle='--')

plt.show()