プログラミング」カテゴリーアーカイブ

Idrisで遊んでみた (2) — 自然数に関する命題の証明

前回は等号の証明に使える公理や規則(Idrisの組み込み関数)をいくつか紹介した。

これに加え、Idrisの標準ライブラリには自然数の性質に関する補題(?)がいくつか用意されている。そのうち、後で使う2つを紹介する。

一つ目は、和と後続者に関するものである。ペアノ算術にはちょうどこのような公理があったはずだ。plusSuccRightSucc

もう一つは、和の可換性である。ペアノ算術ではこれを証明するためには帰納法を使うが、証明(関数の実装)は標準ライブラリで用意されているので、ここではそれを使うだけにとどめておく。plusCommutative

さて、前回証明したいと言っていた命題は\[x=y+\mathtt{difference}\ p\]であった。ただし、\(p\) の表す命題は \(y\le x\), Idrisの型で言うと LTE y x である。この等式を証明するには、
[sourcecode lang=”plain”]
differenceP : {x:Nat} -> {y:Nat} -> (p : LTE y x) -> (x = y + difference p)
[/sourcecode]
という型の関数を実装すれば良い。

例によって場合分けする。

  • \(y=0\), つまり p=lteZero の場合:
    • この等式は自明に成り立つ。refl を返せば良い。
  • それ以外、つまり \(x=Sx’\), \(y=Sy’\), p=lteSucc p' の場合:
    • 帰納法の仮定として differenceP p' : x' = y' + difference p' が使える。
    • この仮定を使って S x' = (S y') + difference p' を示せば良い。(difference の定義を見れば difference p = difference p' なので、S x' = (S y') + difference p は自動的に導いてくれるようだ)
    • この仮定と、先に紹介した公理、補題を使って、自然演繹っぽい証明図を書くと、prooftree1となる。
    • この証明図にIdrisの項を書き加えるとprooftree1-tとなる。ただし、一部をアンダースコアで省略している。

これらをまとめて、Idrisの式として書けば
[sourcecode lang=”plain”]
differenceP : {x:Nat} -> {y:Nat} -> (p : LTE y x) -> (x = y + difference p)
differenceP lteZero = refl
differenceP {x = S x’} {y = S y’} (lteSucc p’) = trans (trans (cong {f=S} (trans (differenceP p’) (plusCommutative _ _))) (plusSuccRightSucc _ _)) (plusCommutative _ (S y’))
[/sourcecode]
となる。一部の式をアンダースコアで省略しているが、最後の (S y') まで省略してしまうと処理系に怒られたので、そこは明示的に書いている。

今回は、自然演繹の証明図をCurry-Howard対応で翻訳することによって、命題を証明(関数を実装)した。等式の変形のところがやや煩雑になっているが、それはもとの自然演繹の証明図がそういうことになっていたことに起因するもので、Curry-Howard対応のせいでどうこうというものではない。

次回は量化子を含む命題を証明する予定である。Tacticを使った証明は筆者が不勉強なのでまだ扱わない。

Idrisで遊んでみた (1)

前回はIdrisで2つの自然数を比較し、結果とその証明を返す関数 compareNats を作った。今回は、それに追加していくつかの関数を実装してみる。

まず、二つの自然数 \(x\) と \(y\) の差、\(x-y\) を計算する関数 difference を作る。\(y\le x\) の場合は良いが、\(x\lt y\) の場合はどうするか考えどころである。0を返すようにしてもよいが、ここは「\(y\le x\) であることの証拠を要求する」ことにして、\(x\lt y\) の場合は考える必要がない、という仕様にする。型を書けば
[sourcecode lang=”plain”]
difference : (x:Nat) -> (y:Nat) -> LTE y x -> Nat
[/sourcecode]
となる。前回書いたように、LTE y x は命題 \(y\le x\) に対応する型である。

実装には、例によって引数によるパターンマッチをする。

  • \(y=0\) の場合
    • \(x\) を返す。第三引数は lteZero である。
    • コードで書けば difference x 0 lteZero = x となる。
  • \(x=Sx’, \>y=Sy’\) の場合
    • 第3引数は lteSucc p の形をしている。
    • difference x' y' p によって \(x’-y’\) を計算できる。\(x-y=x’-y’\) なので、それをそのまま返せば良い。
    • コードで書けば difference (S x') (S y') (lteSucc p) = difference x' y' p となる。

まとめると
[sourcecode lang=”plain”]
difference : (x:Nat) -> (y:Nat) -> LTE y x -> Nat
difference x 0 lteZero = x
difference (S x’) (S y’) (lteSucc p) = difference x’ y’ p
[/sourcecode]
となる。

さて、よく考えてみると、第1引数 x と第2引数 y の情報は第3引数の型を見れば自動的に得られる。したがって、第1引数と第2引数は省略できるのではないか?

Idrisにはimplicit argumentsという機能があるので、この場合実際に xy を省略できる。Implicit argumentsを使うと、difference は次のように書ける:
[sourcecode lang=”plain”]
difference : {x:Nat} -> {y:Nat} -> LTE y x -> Nat
difference {x} lteZero = x
difference (lteSucc p) = difference p
[/sourcecode]

この difference という関数は、\(y \le x\) を満たす任意の自然数 \(x\), \(y\) および \(y \le x\) の「証拠」p : LTE y x に対して、\[x = y + \mathtt{difference}\>p\]という関係を満たす。この関係式を「証明」するにはどうすればよいだろうか?

その前に、述語論理における等号の公理と、それがIdrisではどうなっているか見ておこう。自然演繹チックに書けば、述語論理の等号の公理と関連する推論規則は

  • 等号公理(反射律)equality
  • 等号規則equality2
  • 等号規則(派生)equality2b
  • 推移律equality3
  • 対称律equality4

などがある。あとの三つは最初の二つから出るので、本質的というわけではない。それぞれにIdrisの式を対応させていくと、

  • 等号公理(反射律)equality-t
  • 等号規則equality2-t
  • 等号規則(派生)equality2b-t
  • 推移律equality3-t
  • 対称律equality4-t

となる。つまり、refl, replace, cong, trans, sym という関数があって、それぞれ適切な型を持っている。

以後、命題の左に項を書いていくと上下で重複が多くなって大変なので、以後

  • 等号公理(反射律)equality-u
  • 等号規則equality2-u
  • 等号規則(派生)equality2b-u
  • 推移律equality3-u
  • 対称律equality4-u

のように、使った規則を横棒の右に書くことにする。

長くなるので、一旦記事を分割することにする。今回はここまで。

2月28日:「等号規則(派生)」を追加。

Idrisで遊んでみた (0)

依存型のあるプログラミング言語に慣れてみたいと思ってIdrisで遊んでみた。

依存型とは何かと言うと、普通の型付きラムダ計算はCurry-Howard対応によって直観主義命題論理に対応するが、依存型のあるラムダ計算は直観主義述語論理に対応する(という認識でいる)。

Curry-Howard対応によれば「型」は「命題」に対応する。そして、対応する論理は述語論理なので、「\(x=y\)」や「\(x\le 3\)」のような、項を含む命題(型)を扱うことができる。Integer型やString型をもつ値はまあ分かるが、\(x=y\) や \(x\le 3\) という型をもつ値は一体何を表しているのかと思われるかもしれない。まあ「その命題が成り立つことの証拠、あるいは証明」だと考えればいいだろう。たぶん。

なんでIdrisをやろうかと思ったかというと、文法がHaskellに似ていて取っつきやすそうだったからである。ただ、まだまだ発展途上で実用するには向かなさそうだ。まとまったリファレンスマニュアルみたいなものは見つけられなかったので、標準ライブラリの細かいところはソースコードを参照した。また、この記事のコードはVersion 0.9.11.2に向けて書いてある。

今回はIdrisで自然数(Nat 型)を扱ってみる。最初の自然数0に対応するものは Z : Nat である。ただし、型クラスがいい感じに定義されているので0と書いても良い。自然数 x の後続者(successor)は S x である。さて、さっき命題の例として \(x\le 3\) と書いたが、Idrisには LTE : Nat -> Nat -> Type というデータ構築子があって、自然数についての \(x\le 3\) という命題は LTE x 3 という型に対応する。LT : Nat -> Nat -> Type というものもあるが、これは LT x y = LTE (S x) y と定義されている。それぞれ、”Less Than or Equal to,” “Less Than” の略だと思われる。携帯電話の通信規格とは関係ない。

簡単な例として、二つの自然数を比較する関数を書いてみよう。普通の言語なら compareNats : Nat -> Nat -> Bool とでもするところだが、さっき書いたように「型は命題に対応」し、「値は証明、あるいは証拠を表す」と考えられるので、どうせなら「大きい方はどちらか、およびその証拠」を返す関数を作ってみよう。この関数の型は
[sourcecode lang=”plain”]
compareNats : (x : Nat) -> (y : Nat) -> Either (LT x y) (LTE y x)
[/sourcecode]
とする。この結果の Either (LT x y) (LTE y x) という型は、数学っぽく書けば \((x<y)\vee(y\le x)\) に相当する。

LTE x y型の値を作る(\(x\le y\) を証明する)には、2通りの関数(公理)がある。つまり、

  • lteZero : LTE Z y、つまり \(0\le y\)
  • lteSucc : LTE x y -> LTE (S x) (S y)、つまり \(x \le y \rightarrow Sx \le Sy\)

だ。compareNats の実装では xy の値によって場合分けして、これらの関数を使って「証拠」を作ってやれば良い。

  • \(y=0\) の場合:
    • \((x<y)\vee(y\le x)\) のうち \(y \le x\) が成り立つ。LTE y x 型の値を返したいが、\(y=0\) なので lteZero が欲しい型を持つ値である。
    • したがって、Eitherのデータ構築子(\(\vee\) の導入則)と組み合わせて、Right lteZeroを返せば良い。
  • \(x=0, y=S y’\) の場合:
    • \((x<y)\vee(y\le x)\) のうち \(x<y\) が成り立つ。したがって、LT x y 型、つまり LTE (S Z) (S y') 型の値を返したい。
    • LTE Z y' 型の値であれば lteZero で得られた。これを lteSucc に食わせて lteSucc lteZero とすれば LTE (S Z) (S y') 型の値が得られる。
    • したがって、Eitherのデータ構築子(\(\vee\) の導入則)と組み合わせて、Left (lteSucc lteZero)を返せば良い。
  • \(x=Sx’,y=Sy’\) の場合:
    • まず \(x’\) と \(y’\) を比較する。比較には compareNats x' y' が使える。
    • compareNats x' y'Left p を返してきた場合:
      • p は \(x’ < y’\) の「証拠」である。これの両辺の後続者をとれば(両辺に1加えれば)、\(Sx'<Sy’\) の「証拠」、つまり \(x<y\) の「証拠」になる。
      • つまり lteSucc p が \(x<y\) の「証拠」である。\((x<y)\vee(y\le x)\) のうち成り立つのは \(x<y\) の方なので、Left (lteSucc p) を返せば良い。
    • compareNats x' y'Right q を返してきた場合:
      • q は \(y’ \le x’\) の「証拠」である。これの両辺の後続者をとれば(両辺に1加えれば)、\(Sy’\le Sx’\) の「証拠」、つまり \(y\le x\) の「証拠」になる。
      • つまり lteSucc q が \(y\le x\) の「証拠」である。\((x<y)\vee(y\le x)\) のうち成り立つのは \(y\le x\) の方なので、Right (lteSucc q) を返せば良い。

以上をまとめると次のようなコードになる:
[sourcecode lang=”plain”]
compareNats : (x : Nat) -> (y : Nat) -> Either (LT x y) (LTE y x)
compareNats x Z = Right lteZero
compareNats Z (S y’) = Left (lteSucc lteZero)
compareNats (S x’) (S y’) = case compareNats x’ y’ of
Left p => Left (lteSucc p)
Right q => Right (lteSucc q)
[/sourcecode]
あるいは、either : (a -> c) -> (b -> c) -> (Either a b) -> c 関数を使えば次のようにも書ける:
[sourcecode lang=”plain”]
compareNats : (x : Nat) -> (y : Nat) -> Either (LT x y) (LTE y x)
compareNats x Z = Right lteZero
compareNats Z (S y’) = Left (lteSucc lteZero)
compareNats (S x’) (S y’) = either (Left . lteSucc) (Right . lteSucc) (compareNats x’ y’)
[/sourcecode]

今回はここまで。

なお、今回の記事では

  • Curry-Howard対応
  • Haskell(またはそれに類した言語)
  • ペアノ算術

あたりの知識を仮定した。初心者向けではなく、自分のための備忘録的な感じで書いた。

自動微分(4) — 数学的(幾何学的)に見る

過去の記事:

これまで数回に渡って自動微分を取り上げてきたが、それでは、自動微分は,数学的に見るとどのような概念と対応しているのだろうか?

簡単のため,ここでは実数 \(\mathbf{R}\) についての関数の、1変数による1階微分の自動微分のみを考えるものとする.1変数関数の自動微分を表す型を \(\mathrm{AutoDiff}(\mathbf{R})\) とする.\(\mathrm{AutoDiff}(\mathbf{R})\) には,\(\mathbf{R}\) から誘導された加法や乗法などの演算が入っている.また,\(\mathbf{R}^n\to\mathbf{R}\) の微分可能な関数から,\(\mathrm{AutoDiff}(\mathbf{R}^n)\to\mathrm{AutoDiff}(\mathbf{R})\) の関数が誘導される.加法や乗法などの演算は,\(\mathbf{R}^2\to\mathbf{R}\) の微分可能な関数と見ることができるので,関数の誘導だけを考えれば良い.

つまり,自動微分を考えるということは,実数 \(\mathbf{R}\) あるいはユークリッド空間 \(\mathbf{R}^n\) について,それぞれ \(\mathrm{AutoDiff}(\mathbf{R})\), \(\mathrm{AutoDiff}(\mathbf{R}^n)\) という空間を考えるということである.そして,(微分可能な)関数 \(\mathbf{R}^n\to\mathbf{R}^m\) に対しては,誘導された関数 \(\mathrm{AutoDiff}(\mathbf{R}^n)\to\mathrm{AutoDiff}(\mathbf{R}^m)\) が存在する.
これは圏論の言葉で言う「関手」(functor)によく似ている.しかし,\(\mathrm{AutoDiff}\) を圏の関手として見る場合,対象とする空間をユークリッド空間 \(\mathbf{R}^n\) に制限しているのは気持ち悪い.1階微分が定義できるようなもっと広い空間に対して,自動微分を一般化して考えられないだろうか?

数学では,ユークリッド空間を一般化したような空間を「多様体」(manifold)と呼ぶ.特に,「滑らか」で「関数の微分が定義できる」ものを「滑らかな多様体」とか「可微分多様体」と呼ぶ.多様体についての詳細はここでは省くので、数学の授業を受けるか、専門書を参照するかして欲しい。何が言いたいかというと,自動微分のようなものを,可微分多様体の圏から可微分多様体の圏への関手として定義できるはずである.

実は,この概念はすでに多様体論ではよく知られていて,「接バンドル」(tangent bundle)と呼ばれている.\(\mathrm{AutoDiff}(\mathbf{R}^n)\) に対応する空間は \(T\mathbf{R}^n\) と書かれる.一般に,多様体 \(M\) の接バンドルは \(TM\) と表される.写像の方は,\(f\colon M\to N\) から誘導される写像は \(df\colon TM\to TN\) と書かれる.

\(\mathbf{R}\) の自動微分、\(\mathrm{AutoDiff}(\mathbf{R})\) は \(\mathbf{R}\times\mathbf{R}\) のような直積で表された。このような場合、 \(T\mathbf{R}\) は自明(trivial)であるというような言い方をする。また、\(\mathrm{AutoDiff}(\mathbf{R}^n)\cong\mathrm{AutoDiff}(\mathbf{R})^n\)(自然な同型)となるが、これは接バンドルについて \(T(M\times N)\cong TM\times TN\)(自然な同型)となることに対応する。

筆者の知識が足らないので今書けるのはここまでだが、要するにプログラミングのテクニックの一つである自動微分は、幾何学の概念と対応していそうだということである。あとそれから、幾何ガチ勢からまさかりが飛んできそうで怖い。

自動微分(3) — 高階導関数の計算

前回まででは、1階の(偏)微分係数を考察した。今度は、高階の導関数の値(微分係数)も計算するように拡張しよう。簡単のため、1変数関数について考えることにする。

高階微分と演算、関数の合成との関係を見る。変数は \(t\) とする。

もとの関数 1階導関数、2階導関数
\(h(t)=\mathrm{const.}\) \(h'(t)=0,\quad h^{\prime\prime}(t)=0\)
\(h(t)=t\) \(h'(t)=1,\quad h^{\prime\prime}(t)=0\)
\(h(t)=f_1(t)\pm f_2(t)\) \(h'(t)=f_1′(t)\pm f_2′(t),\quad h^{\prime\prime}(t)=f_1^{\prime\prime}(t)\pm f_2^{\prime\prime}(t)\)
\(h(t)=f_1(t)f_2(t)\) \(\begin{array}{c}
h'(t)=f_1′(t)f_2(t)+f_1(t)f_2′(t), \\
h^{\prime\prime}(t)=f_1^{\prime\prime}(t)f_2(t)+2f_1′(t)f_2′(t)+f_1(t)f_2^{\prime\prime}(t)
\end{array}\)
\(h(t)=\frac{f_1(t)}{f_2(t)}\) \(\begin{array}{c}
h'(t)=\frac{f_1′(t)}{f_2(t)}-\frac{f_1(t)f_2′(t)}{f_2(t)^2}, \\
h^{\prime\prime}(t)=\frac{f_1^{\prime\prime}(t)}{f_2(t)}-\frac{2f_1′(t)f_2′(t)}{f_2(t)^2}-\frac{f_1′(t)f_2^{\prime\prime}(t)}{f_2(t)^2}+\frac{2f_1(t)f_2′(t)^2}{f_2(t)^3}
\end{array}\)
\(h(t)=g(f(t))\) \(\begin{array}{c}
h'(t)=g'(f(t))f'(t), \\
h^{\prime\prime}(t)=g^{\prime\prime}(f(t))f'(t)^2+g'(f(t))f^{\prime\prime}(t)
\end{array}\)
\(h(t)=g(f_1(t),f_2(t))\) \(\begin{array}{c}
h'(t)=g_x(f_1(t),f_2(t))f_1′(t)+g_y(f_1(t),f_2(t))f_2′(t), \\
h^{\prime\prime}(t)=g_{xx}(f_1(t),f_2(t))f_1′(t)^2+2g_{xy}(f_1(t),f_2(t))f_1′(t)f_2′(t)+g_{yy}(f_1(t),f_2(t))f_2′(t)^2 \\
\quad +g_x(f_1(t),f_2(t))f_1^{\prime\prime}(t)+g_y(f_1(t),f_2(t))f_2^{\prime\prime}(t)
\end{array}\)

2階微分の自動微分では、「値」「1階の微分係数」「2階の微分係数」の3つの値を保持する必要があるので、実装するには \(\mathbf{R}\times\mathbf{R}\times\mathbf{R}\) に適切な演算を入れれば良い。上の表の \(f(t),f'(t),f^{\prime\prime}(t)\) を \(\mathbf{R}\times\mathbf{R}\times\mathbf{R}\) の成分で置き換えると、\(\mathbf{R}\times\mathbf{R}\times\mathbf{R}\) に入れるべき演算が分かる。すなわち、\begin{align*}
(x,x’,x^{\prime\prime})\pm(y,y’,y^{\prime\prime})&=(x\pm y,x’\pm y’,x^{\prime\prime}\pm y^{\prime\prime}) \\
(x,x’,x^{\prime\prime})\cdot(y,y’,y^{\prime\prime})&=(xy,x’y+xy’,x^{\prime\prime}y+2x’y’+xy^{\prime\prime}) \\
\frac{(x,x’,x^{\prime\prime})}{(y,y’,y^{\prime\prime})}&=\left(\frac{x}{y},\frac{x’}{y}-\frac{xy’}{y^2},\frac{x^{\prime\prime}}{y}-\frac{2x’y’}{y^2}-\frac{x’y^{\prime\prime}}{y^2}+\frac{2xy’^2}{y^3}\right),
\end{align*}\(g\colon\mathbf{R}\to\mathbf{R}\)(微分可能) に対し、\[g_*((x,x’,x^{\prime\prime}))=(g(x),g'(x)x’,g^{\prime\prime}(x)x’^2+g'(x)x^{\prime\prime}),\]\(g\colon\mathbf{R}\times\mathbf{R}\to\mathbf{R}\)(微分可能) に対し、\[g_*((x,x’,x^{\prime\prime}),(y,y’,y^{\prime\prime}))=\left(\begin{array}{c}g(x,y), \qquad
g_x(x,y)x’+g_y(x,y)y’, \\
g_{xx}(x,y)x’^2+2g_{xy}(x,y)x’y’+g_{yy}(x,y)y’^2+g_x(x,y)x^{\prime\prime}+g_y(x,y)y^{\prime\prime}\end{array}\right),\]
とすればよい。ただし、変数についている \({}^\prime\) は微分ではなく変数名の一部と考える。関数についている \({}^\prime\) や添字は微分の意である。

これを実際のプログラミング言語で実装するのは容易だろう。

自動微分(2) — 偏微分係数への一般化

前回まで考察したのは,1変数関数とその導関数の値を計算する場合だった.今回は,2変数関数とその偏導関数を計算する場合について考える.

自動微分を考える前に、2変数関数の微分(偏微分)と演算、関数の合成との関係を見ておこう。2つの変数を \(u\), \(v\) とおく.偏導関数は,変数を添字にして表すことにする.

もとの関数 偏導関数
\(h(u,v)=\mathrm{const.}\) \(h_u(u,v)=0,\quad h_v(u,v)=0\)
\(h(u,v)=u\) \(h_u(u,v)=1,\quad h_v(u,v)=0\)
\(h(u,v)=v\) \(h_u(u,v)=0,\quad h_v(u,v)=1\)
\(h(u,v)=f_1(u,v)\pm f_2(u,v)\) \(\begin{aligned}
h_u(u,v)&=f_{1u}(u,v)\pm f_{2u}(u,v), \\
h_v(u,v)&=f_{1v}(u,v)\pm f_{2v}(u,v)
\end{aligned}\)
\(h(u,v)=f_1(u,v)f_2(u,v)\) \(\begin{aligned}
h_u(u,v)&=f_{1u}(u,v)f_2(u,v)+f_1(u,v)f_{2u}(u,v), \\
h_v(u,v)&=f_{1v}(u,v)f_2(u,v)+f_1(u,v)f_{2v}(u,v)
\end{aligned}\)
\(h(u,v)=f_1(u,v)/f_2(u,v)\) \(\begin{aligned}
h_u(u,v)&=(f_{1u}(u,v)f_2(u,v)-f_1(u,v)f_{2u}(u,v))/f_2(u,v)^2, \\
h_v(u,v)&=(f_{1v}(u,v)f_2(u,v)-f_1(u,v)f_{2v}(u,v))/f_2(u,v)^2
\end{aligned}\)
\(h(u,v)=g(f(u,v))\) \(\begin{aligned}
h_u(u,v)&=g'(f(u,v))f_u(u,v), \\
h_v(u,v)&=g'(f(u,v))f_v(u,v)
\end{aligned}\)
\(h(u,v)=g(f_1(u,v),f_2(u,v))\) \(\begin{aligned}
h_u(u,v)&=g_x(f_1(u,v),f_2(u,v))f_{1u}(u,v)+g_y(f_1(u,v),f_2(u,v))f_{2u}(u,v), \\
h_v(u,v)&=g_x(f_1(u,v),f_2(u,v))f_{1v}(u,v)+g_y(f_1(u,v),f_2(u,v))f_{2v}(u,v)
\end{aligned}\)

2変数の自動微分を表すときに必要な情報は

  • 関数の値 \(f(u,v)\)
  • 関数の \(u\) による偏微分係数 \(f_u(u,v)\)
  • 関数の \(v\) による偏微分係数 \(f_v(u,v)\)

の3つである。そこで、\(\mathbf{R}\) の自動微分を表す型を \(\mathrm{AutoDiff2}\mathbf{R}=\mathbf{R}\times\mathbf{R}\times\mathbf{R}\) とおき、演算を定義する。上の関係式の \(f_1\),\(f_2\) とその微分係数を \((x,x_u,x_v),(y,y_u,y_v)\in\mathrm{AutoDiff2}\mathbf{R}\) で置き換えると、\(\mathrm{AutoDiff2}\mathbf{R}\) に入れるべき演算が分かる:\begin{align*}
(x,x_u,x_v)\pm(y,y_u,y_v)&=(x\pm y,x_u\pm y_u,x_v\pm y_v) \\
(x,x_u,x_v)\cdot(y,y_u,y_v)&=(x y,x_u y+x y_u,x_v y+x y_v) \\
(x,x_u,x_v)/(y,y_u,y_v)&=(x/y,(x_u y-x y_u)/y^2,(x_v y-x y_v)/y^2)
\end{align*}\(g\colon\mathbf{R}\to\mathbf{R}\) に対し,\[
g_*((x,x_u,x_v))=(g(x),g'(x)x_u,g'(x)x_v),\]
\(g\colon\mathbf{R}\times\mathbf{R}\to\mathbf{R}\) に対し,\[
g_*((x,x_u,x_v),(y,y_u,y_v))=(g(x,y),g_x(x,y)x_u+g_y(x,y)y_u,g_x(x,y)x_v+g_y(x,y)y_v)
\]

Haskellで実装すると,たとえば次のようになる:
[sourcecode lang=”plain”]
data AutoDiff2 a = AutoDiff2 a a a deriving(Eq)
instance Num a => Num (AutoDiff2 a) where
(AutoDiff2 f f_u f_v) + (AutoDiff2 g g_u g_v) = AutoDiff2 (f+g) (f_u+g_u) (f_v+g_v)
(AutoDiff2 f f_u f_v) – (AutoDiff2 g g_u g_v) = AutoDiff2 (f-g) (f_u-g_u) (f_v-g_v)
(AutoDiff2 f f_u f_v) * (AutoDiff2 g g_u g_v) = AutoDiff2 (f*g) (f_u*g+f*g_u) (f_v*g+f*g_v)
abs _ = undefined
signum _ = undefined
fromInteger n = AutoDiff2 (fromInteger n) 0 0
[/sourcecode]

さて、1変数の場合の自動微分(\(\mathbf{R}\times\mathbf{R}\) に演算を入れる)についてもう一度見てみよう。\begin{align*}
(x_0,x_1)\pm(y_0,y_1)&=(x_0\pm y_0,x_1\pm y_1) \\
(x_0,x_1)\cdot(y_0,y_1)&=(x_0 y_0,x_1 y_0+x_0 y_1) \\
(x_0,x_1)/(y_0,y_1)&=(x_0/y_0,(x_1 y_0-x_0 y_1)/y_0^2),
\end{align*}\(g\colon\mathbf{R}\to\mathbf{R}\)(微分可能) に対し、\[g_*((x_0,x_1))=(g(x_0),g'(x_0)x_1),\]\(g\colon\mathbf{R}\times\mathbf{R}\to\mathbf{R}\)(微分可能) に対し、\[g_*((x_0,x_1),(y_0,y_1))=(g(x_0,y_0),g_x(x_0,y_0)x_1+g_y(x_0,y_0)y_1),\]
であった。
よく見てみると、第2成分の \(x_1\), \(y_1\) は \(\mathbf{R}\) でなくても,\(\mathbf{R}\)-ベクトル空間なら何でもよいことが分かる。2変数の偏微分の自動微分は、\(x_1\), \(y_1\) として \(\mathbf{R}^2\) の元をとったものだと考えることができる。一般に、\(n\) 変数の自動微分は \(x_1\), \(y_1\) として \(\mathbf{R}^n\) の元をとることにより得られる。

自動微分(1) — 1変数の場合

前回の記事が適当すぎたので仕切り直しを。

「値と微分係数を同時に計算する型」AutoDiff型は,\(\mathbf{R}\times\mathbf{R}\)に適当な演算を入れたものと考えられる。微分を考える変数を \(x\),微分を計算したい点を \(t\) とするとき,第一成分を「\(x=t\) における値」,第二成分を「\(x=t\) における微分係数」と考える。AutoDiff型の値を \(y=(y_0,y_1)\) とすると,\(y\) はある関数 \(f\) について \(y=(f(t),f'(t))\) となっている。\(h(x)=g(f(x))\) という関数があったとき,すでに得られている \(f(t),f'(t)\) の値を使って \(h(t),h'(t)\) を計算したいとする。合成関数の微分は\[h'(x)=g'(f(x))f'(x)\]で与えられるので、ここに \(x=t\) を代入すれば\[h'(t)=g'(f(t))f'(t)=g'(y_0)y_1\]となる。\(h(t)\) はもちろん \(h(t)=g(y_0)\) で与えられる。

今度は \(k(x)=g(l(x))\) という関数を考えてみよう。\(x=t\) における \(l\) の値と \(l’\) の値はすでに計算されているとして、これを \(z=(z_0,z_1)=(l(t),l'(t))\) とおく。このとき,先ほどと同じように \(k(t)=g(z_0)\),\(k'(t)=g'(z_0)z_1\) がわかる。

このようにして,\(g\) という関数は、AutoDiff型からAutoDiff型への関数を定める。この関数を \(g_*\) とすると,AutoDiff型の値 \((y_0,y_1)\) を与えた時の \(g_*\) の値は\[g_*(y_0,y_1)=(g(y_0),g'(y_0)y_1)\]となる。AutoDiff型が「これまで計算した関数の値と微分係数」を表す型だとすれば,\(g_*\) という関数は「これまで計算した関数に \(g\) という関数を合成した時の、関数の値と微分係数」を表す型である。

上の説明では \(g\) と \(g_*\) を区別して書いたが、実際のプログラミングでは関数や演算子のオーバーロード、あるいは型クラスを使って、\(g_*\) に相当する関数も同じ演算子・関数名で書けるようにするのが普通である。

いくつかの演算について、AutoDiff型を \(\mathbf{R}\times\mathbf{R}\) とみなしたときの計算方法を与えておく。\begin{align*}
(x_0,x_1)\pm(y_0,y_1)&=(x_0\pm y_0,x_1\pm y_1) \\
(x_0,x_1)\cdot(y_0,y_1)&=(x_0 y_0,x_1 y_0+x_0 y_1) \\
(x_0,x_1)/(y_0,y_1)&=(x_0/y_0,(x_1 y_0-x_0 y_1)/y_0^2),
\end{align*}\(g\colon\mathbf{R}\to\mathbf{R}\)(微分可能) に対し、\[g_*((x_0,x_1))=(g(x_0),g'(x_0)x_1),\]\(g\colon\mathbf{R}\times\mathbf{R}\to\mathbf{R}\)(微分可能) に対し、\[g_*((x_0,x_1),(y_0,y_1))=(g(x_0,y_0),g_x(x_0,y_0)x_1+g_y(x_0,y_0)y_1),\]となる。

文章書くのだるい。気が向いたら続く。

自動微分(0) — 導入

関数の値 \(f(t)\) とその微分係数 \(f'(t)\) を同時に計算したい状況があるとする.ここでは,具体的な2次関数
\[f(x)=x^2+2x\]
で話を進める.この関数は,C++を使うと次のように実装できる.

[sourcecode lang=”cpp”]
template<typename T>
constexpr T f(T x) {
return x*x+2*x;
}
[/sourcecode]

あるいは,Haskellでは次のように書けるだろう.

[sourcecode lang=”plain”]
f :: Num a => a -> a
f x = x*x+2*x
[/sourcecode]

さて,賢明な読者は,\(f\)の導関数が\(2x+2\)であることに気づいたかもしれない.そうすると,\(f(x)\)と\(f'(x)\)を同時に計算する関数\(g\)は次のように実装できる.

[sourcecode lang=”cpp”]
// in C++
template<typename T>
std::pair<T,T> g(T x) {
return {x*x+2*x, 2*x+2};
}
[/sourcecode]

 

[sourcecode lang=”plain”]
— in Haskell
g :: Num a => a -> (a,a)
g x = (x*x+2*x, 2*x+2)
[/sourcecode]

しかし,せっかく\(f(x)\)の定義式をプログラム中に書いているのに,\(f'(x)\)の計算式を改めて書くのは無駄,と言わざるを得ない.また,プログラムを変更して\(f(x)\)の定義式を書き換えたときに\(f'(x)\)の式を直し忘れる,というミスが起こる可能性もある.できることなら,\(f(x)\)の定義式を使って,\(f'(x)\)を自動的に計算させたい.

そこで登場するのが,自動微分(automatic differentiation)と呼ばれるテクニックである。\(f(x)\)の定義の中で使っているのは和と積だけである.和と積の微分は高校で習ったと思うが,次の式で計算できる.

もとの関数 導関数
\(h(x)=f(x)\pm g(x)\) \(h'(x)=f'(x)\pm g'(x)\)
\(h(x)=f(x)g(x)\) \(h'(x)=f'(x)g(x)+f(x)g'(x)\)

では, f の中での和と積の挙動をいじくって,微分もついでに計算させることはできないだろうか?

幸い,C++には演算子オーバーロードという便利な機能がある.先のコードで, f はテンプレートを使って定義したので,引数の型は intdouble でなくてもよい.演算子をオーバーロードしてやった自前の型を与えても良い.そこで, AutoDiff というクラスを次のように定義しよう:

[sourcecode lang=”cpp”]
template<typename T>
class AutoDiff
{
public:
T value; // 値
T diff; // 微分係数の値
constexpr AutoDiff() : value(), diff() {} // 定数0
constexpr AutoDiff(T value) : value(value), diff() {} // 定数で初期化(微分係数は0)
constexpr AutoDiff(T value, T diff) : value(value), diff(diff) {} // 値 value, 微分係数 diff で初期化

friend constexpr AutoDiff<T> operator+(AutoDiff<T> f, AutoDiff<T> g) {
return {f.value+g.value, f.diff+g.diff};
}
friend constexpr AutoDiff<T> operator*(AutoDiff<T> f, AutoDiff<T> g) {
return {f.value*g.value, f.diff*g.value+f.value*g.diff};
}
};
[/sourcecode]

このクラスを使って,先ほどの関数 g を実装しなおしてみると,次のようになる:

[sourcecode lang=”cpp”]
template<typename T>
std::pair<T,T> g(T x) {
AutoDiff<T> y = f(AutoDiff<T>(x, 1));
return {y.value, y.diff};
}
[/sourcecode]

実際にこの定義し直した g が関数の値と微分係数を正しく計算することは,各自確かめてほしい.

Haskellでは,型クラスを使うと同じようなことができる.

[sourcecode lang=”plain”]
data AutoDiff a = AutoDiff a a deriving(Eq)
instance Num a => Num (AutoDiff a) where
(AutoDiff f f’) + (AutoDiff g g’) = AutoDiff (f+g) (f’+g’)
(AutoDiff f f’) – (AutoDiff g g’) = AutoDiff (f-g) (f’-g’)
(AutoDiff f f’) * (AutoDiff g g’) = AutoDiff (f*g) (f’*g+f*g’)
abs _ = undefined
signum _ = undefined
fromInteger = flip AutoDiff 0 . fromInteger

g :: Num a => a -> (a,a)
g x = case f (AutoDiff x 1) of AutoDiff f f’ -> (f,f’)
[/sourcecode]

ここでのコード例は加法と乗法しか定義しなかったが,残りの演算は,次の表を参考にすれば容易に実装できるだろう.

もとの関数 導関数
\(h(x)=f(x)\pm g(x)\) \(h'(x)=f'(x)\pm g'(x)\)
\(h(x)=f(x)g(x)\) \(h'(x)=f'(x)g(x)+f(x)g'(x)\)
\(h(x)=f(x)/g(x)\) \(h'(x)=(f'(x)g(x)-f(x)g'(x))/g(x)^2\)
\(h(x)=g(f(x))\) \(h'(x)=g'(f(x))f'(x)\)
\(h(x)=g(f_1(x),f_2(x))\) \(h'(x)=g_u(f_1(x),f_2(x))f_1′(x)+g_v(f_1(x),f_2(x))f_2′(x)\)
ただし\(g_u,g_v\)は\(g\)の偏導関数

文章を書くのに疲れたので,今日はここまで.(続く)