自動微分

Conal ElliotのBeautiful differentiationを見て思い出したこと。
http://conal.net/papers/beautiful-differentiation/

昔
http://www.bcl.hamilton.ie/~barak/papers/ifl2005.pdf
を見て考えたのは、ここにある定義に従って、自動微分を定義すると、例えばNewton法は

newton f x0 = iterate (\x-> x - (normal.f.lift $ x)/(d f x)) x0
  where
   normal (D  x _) = x

のような感じで定義される。しかしliftとかださい。

論文のテーマになっている

d (\ x -> x * d (\ y -> (lift x) + y) 1) 1

とかでも、入れ子が複雑になると、どこにlift付けるか訳が分からなくなりそう

というわけで自動微分の型を、

deriv :: (forall x.(Floating x)=>x->x) -> (forall a.(Floating a)=>(a->a))

にしとくと、色々エレガントになるんじゃないかと思ったんだけど

deriv (\ x -> x * deriv (\ y -> x + y) 1) 1

とか型推論エラーではねられるようになってしまった。

Newton法は、キレイになるんだけど(型を書く必要があるけれど)

newton :: (Floating a) => (forall x.(Floating x)=>x->x) -> a -> [a]
newton f x0 = iterate (\x-> x - (f x)/(deriv f x)) x0

どういう風になれば、この問題が解決できるようになるのか、よく分からない。

あと、Conal Elliotの論文にHaskellでベクトル空間をどう定義するのか的な話があった。ベクトルはリストでも配列でもないので、昔どっかの論文で、内積は = sum $ zipWith(*) as bsと書けるよねみたいな記述があって、それはダメじゃん、じゃあどうすんの？とか思って、基底の集合に係数を対応させるようなマップとして定義すると綺麗になるとこまでは行き着くというか
http://d.hatena.ne.jp/m-a-o/20070704#p1
のパクリだけど、それだと次元が2とか3とかの時はともかく、100とか1000とかなったら、基底の集合の型とかどうやって作るのか、みたいな問題が発生して結局これもどうすればよいのか、よく分からない

以下実装

{-# OPTIONS -fglasgow-exts #-}
data Num a=> Diff a = D a a

instance (Num a) => Eq (Diff a) where
  (D x x') == (D y y') = (x==y) && (x'==y')

instance (Num a) => Show (Diff a) where
  show (D x x') = show (x,x')

instance (Num a) => Num (Diff a) where
  (D x x') + (D y y') = D (x+y) (x' + y')
  (D x x') * (D y y') = D (x*y) (x*y' + x'*y)
  (D x x') - (D y y') = D (x-y) (x'-y')
  negate (D x x') = D (-x) (-x')
  abs (D x x') |(abs x /= x) = D (-x) (-x')
               |(x==0)&&(abs x' /= x') = D (-x) (-x')
               |otherwise = D x x'
  signum (D x x') = D (signum x) 0
  fromInteger x = D (fromInteger x) 0

instance (Fractional a) => Fractional (Diff a) where
  fromRational x = D (fromRational x) 0
  recip (D x x') = D (recip x) (-x'/(x*x))

instance (Floating a) => Floating (Diff a) where
  pi = D pi 0
  exp (D x x') = D (exp x) (x' * exp x)
  log (D x x') = D (log x) (x'/x)
  sin (D x x') = D (sin x) (x' * cos x)
  cos (D x x') = D (cos x) (-x' * sin x)
  sinh (D x x') = D (sinh x) (x' * cosh x)
  cosh (D x x') = D (cosh x) (x' * sinh x)
  asin (D x x') = D (asin x) (x' /(sqrt $ 1-x*x))
  acos (D x x') = D (acos x) (-x'/(sqrt $ 1-x*x))
  atan (D x x') = D (atan x) (x'/(x*x+1))
  asinh (D x x') = D (asinh x) (x'/(sqrt $ 1+x*x))
  acosh (D x x') = D (acosh x) (x'/(sqrt $ 1-x*x))
  atanh (D x x') = D (atanh x) (x' /(1-x*x))

--微分
deriv :: (forall x.(Floating x)=>x->x) -> (forall a.(Floating a)=>(a->a))
deriv f = (\x -> mt $ f (D x 1))
 where mt (D _ x) = x