SICP 2.3.2 Example: Symbolic Differentiation - プログラミング再入門

数値計算的な微分ではなくて代数的に数式そのものを扱った微分。
Lispにとっては昔から重要なテーマで、シンボルを扱えるプログラミング言語を開発する動機になっていたとの事。まぁ、確かにFORTRANでは難しい。

ノート

The differentiation program with abstract data

抽象化されたデータを使った微分プログラム。
$\frac{d(u+v)}{dx}=\frac{du}{dx}+\frac{dv}{dx}$
$\frac{d(uv)}{dx}=u\frac{dv}{dx}+v\frac{du}{dx}$
uとかvとか何かなと思ったら、基本的にはxを含んだ式（含まないかも知れないけど）。

(define (deriv exp var)
  (cond ((number? exp) 0)
        ((variable? exp)
         (if (same-variable? exp var) 1 0))
        ((sum? exp)
         (make-sum (deriv (addend exp) var)
                   (deriv (augend exp) var)))
        ((product? exp)
         (make-sum
           (make-product (multiplier exp)
                         (deriv (multiplicand exp) var))
           (make-product (deriv (multiplier exp) var)
                         (multiplicand exp))))
        (else
         (error "unknown expression type -- DERIV" exp))))

プログラムを解釈すると、導関数とは：

expがただの数字だったら０
expが変数だったら、
1. 微分しようとしている変数だったら1、そうでなければ0
expが足し算の式だったら、最初の項をの導関数と２番目の項の導関数の和の式
expが掛け算の式だったら、掛ける数（×の後ろ）と掛けられる数（×の前）の導関数の積と、掛ける数の導関数と掛けられる数の積の式の和のしき
取り敢えずそれ以外はエラー

DrRacketのR5RSモードでerrorを使うのは色々とややこしい話が有る様なのでSchemeモード（昔のDrSchemeの時との互換モード？）を使う。
このプログラムは数式がどの様な方法で実装されていても抽象化された（コンストラクタとアクセサで実装方法が隠されている）数式を扱っているので動作する。

Representing algebraic expressions

代数式を表現（実装）する。
ここでは中置記法ではなく普通のLispの記法で表現する。
ここまで出て来た事はない気がするが、(cadr s)は(car (cdr s))つまりリストの２番目の要素、(caddr s)は(car (cdr (cdr s)))つまりリストの３場面目の要素。
これまでの定義を動作させてみる。Racketモードなので結果はクオートされた形式で表示される。

> (deriv '(+ x 3) 'x)
(+ 1 0)
> (deriv '(* x y) 'x)
(+ (* x 0) (* 1 y))
> (deriv '(* (* x y) (+ x 3)) 'x)
(+ (* (* x y) (+ 1 0)) (* (+ (* x 0) (* 1 y)) (+ x 3)))
>

make-sum、make-productで式の無駄な部分を省略する実装に換えて、動作確認。

> (deriv '(+ x 3) 'x)
1
> (deriv '(* x y) 'x)
y
> (deriv '(* (* x y) (+ x 3)) 'x)
(+ (* x y) (* y (+ x 3)))
>

$\frac{d}{dx}xy(x+3)$
$=xy+y(x+3)$
これを
$=2xy+3y$
とするか
$=y(2x+3)$
まで持って行くかは目的によるので、式の整理はなかなか難しい問題。

Exercise 2.56

引き算の式が作れないので、とりあえず(make-sum a -1)と言う形にする。

(define (base s) (cadr s))
(define (exponent s) (caddr s))

(define (make-exponentiation b e)
  (cond ((=number? b 0) 0)
        ((or (=number? b 1) (=number? e 0)) 1)
        ((=number? e 1) b)
        (else (list '** b e))))

(define (deriv exp var)
  (cond ((number? exp) 0)
        ((variable? exp)
         (if (same-variable? exp var) 1 0))
        ((sum? exp)
         (make-sum (deriv (addend exp) var)
                   (deriv (augend exp) var)))
        ((product? exp)
         (make-sum
           (make-product (multiplier exp)
                         (deriv (multiplicand exp) var))
           (make-product (deriv (multiplier exp) var)
                         (multiplicand exp))))
        ((exponentiation? exp)
         (make-product
          (make-product (exponent exp)
                        (make-exponentiation (base exp) (make-sum (exponent exp) -1)))
          (deriv (base exp) var)))
        (else
         (error "unknown expression type -- DERIV" exp))))

動作確認

> (exponentiation? '(** a b))
#t
> (exponentiation? '(* a b))
#f
> (base '(** a b))
a
> (exponent '(** a b))
b
> (make-exponentiation 'a 1)
a
> (make-exponentiation 'a 0)
1
> (make-exponentiation 2 16)
(** 2 16)
> (make-exponentiation 2 'b)
(** 2 b)
> (make-exponentiation 'a 3)
(** a 3)
> (make-exponentiation 1 'b)
1
> (make-exponentiation 0 'b)
0
> (deriv '(** 2 3) 'x)
0
> (deriv '(** x 3) 'x)
(* 3 (** x 2))
> (deriv '(** x y) 'x)
(* y (** x (+ y -1)))
> (deriv '(* (+ 3 y) (+ y (** x 4))) 'x)
(* (+ 3 y) (* 4 (** x 3)))
> (deriv '(+ (** x 2) (* 3 x)) 'x)
(+ (* 2 x) 3)
>

$(2^3)'=0$
$(x^y)'=yx^{y-1}$
$\left{(3+y)(y+x^4)\right}'=(3+y)4x^3$
$(x^2+3)'=2x+3$

Exercise 2.57

要するに

> (augend '(+ x (+ y 2)))
(+ y 2)
>

と言う動作を

'(x y 2)

に対しても同じ結果になる様にしたい。このインターフェースを使っている限り他の部分の変更は不要。

(define (the-rest exp op)
  (if (pair? exp)
      (if (null? (cdr exp)) (car exp)
          (cons op exp))
      (exp)))
(define (augend s)
  (the-rest (cddr s) '+))
(define (multiplicand p)
  (the-rest (cddr p) '*))

動作確認

> (deriv '(+ (** x 2) (* 3 x) 4) 'x)
(+ (* 2 x) 3)
>

（掛け算の方は使ってないけど）

Exercise 2.58

中置記法を実現したい。
a)
ここでは引数の数は必ず２個。つまり(1+2+3)の様に連結されずに必ず( (1+2)+3)の様に演算子の引数が必ず２個になる様な式になっている前提。
まずはコンストラクタ。

(define (make-sum a1 a2)
  (cond ((=number? a1 0) a2)
        ((=number? a2 0) a1)
        ((and (number? a1) (number? a2)) (+ a1 a2))
        (else (list a1 '+ a2))))
(define (make-product m1 m2)
  (cond ((or (=number? m1 0) (=number? m2 0)) 0)
        ((=number? m1 1) m2)
        ((=number? m2 1) m1)
        ((and (number? m1) (number? m2)) (* m1 m2))
        (else (list m1 '* m2))))
(define (make-exponentiation b e)
  (cond ((=number? b 0) 0)
        ((or (=number? b 1) (=number? e 0)) 1)
        ((=number? e 1) b)
        (else (list b '** e))))

最後にリストを形成する部分だけが違う。
動作確認

> (make-sum 1 2)
3
> (make-sum 'x 2)
(x + 2)
> (make-product 'x 2)
(x * 2)
> (make-product 'x 1)
x
> (make-exponentiation 'x 2)
(x ** 2)
> (make-exponentiation 'x 0)
1
>

次に演算子の判定

(define (infix? expr op)
  (and (pair? expr)
       (pair? (cdr expr))
       (eq? (cadr expr) op)))
(define (sum? x)
  (infix? x '+))
(define (product? x)
  (infix? x '*))
(define (exponentiation? x)
  (infix? x '**))

項の間に置いてある演算子が指定の演算子と一致するか判定するinfix?を定義。
アクセサはそれぞれ以下の通り。

(define (addend s)
  (car s))
(define (augend s)
  (caddr s))
(define (multiplier p)
  (car p))
(define (multiplicand p)
  (caddr p))
(define (base e)
  (car e))
(define (exponent s)
  (caddr s))

引数は２個限定とするのでaugend、multiplicand、exponentは元の定義のまま。
動作確認

> (addend '(x + 1))
x
> (augend '(x + 1))
1
> (addend '((x + 1) + (y + 2)))
(x + 1)
> (augend '((x + 1) + (y + 2)))
(y + 2)
>

これでderivは変更無しで動作する筈。動作確認

> (deriv '(x * (y * (x + 3))) 'x)
((x * y) + (y * (x + 3)))
> (deriv '((x ** 2) + ((3 * x) + 4)) 'x)
((2 * x) + 3)
>

数式で追うと：
$\left[x{y(x+3)}\right]'=x\left{y(x+3)\right}'+x'\left{y(x+3)\right}=x\left{y(x+3)'+y'(x+3)\right}+y(x+3)=xy+y(x+3)$
$\left{x^2+(3x+4)\right}'=(x^2)'+(3x+4)'=2x+(3x)'+4'=2x+3$
きちんと再帰して各項を微分していそう。

b)
引数は必ず２個の縛りを外す。コンストラクタ、セレクタ、述語を変更してderivをそのまま動作させるにはどうするか。
元々式の最初の演算子に対しその第１項と第２項を取り出して導関数を導き出していた。これを式の最初に現れる最も優先度の低い演算子に対しその第１項と第２項を取り出して導関数を導き出す方法に変更する。これは式を分解して細かくして行く際に優先度＝結合度が低い所から分解される為。最も優先度の低い演算子で式を分割する事で優先度の高い演算に括弧が付いているのと同じ評価の仕方になる。分解の必要のない単項の処理は従来通り。
つまりsum?は最初に現れる+を探して、それよりも前をaddend、後ろをaugendとする。product?も同様。
この方法であればderivの定義はそのままで行けそう。

まず第２項を取り出すにはmemqを使えば良さそうな事は簡単に思いつく。演算子を取り除くのにcdrを使う。

> (memq '+ '(1 + 2))
(+ 2)
> (cdr (memq '+ '(1 + 2)))
(2)
> (cdr (memq '+ '(1 + 2 + 3)))
(2 + 3)
> (cdr (memq '+ '(1 + (2 + 3))))
((2 + 3))
>

これだとただの数字なのにリストになっていたり、一段余計なリストに入っていたりでderivでは扱えない状態になる。
そこで以下の手続きthe-restを定義する。

(define (the-rest op s)
  (let ((r (cdr (memq op s))))
    (if (null? (cdr r))
        (car r)
        r)))

(cdr (memq op s))の結果が単項のリストであれば項のみを返し、複数の項が有る場合にはリストのままとする。

> (the-rest '+ '(1 + 2))
2
> (the-rest '+ '(1 + 2 + 3))
(2 + 3)
> (the-rest '+ '(1 + (2 + 3)))
(2 + 3)
>

これを使うとaugend、multiplicand、exponentが定義出来る。

(define (augend s)
  (the-rest '+ s))
(define (multiplicand p)
  (the-rest '* s))
(define (exponent e)
  (the-rest '** e))

これと同様な動きで第１項を取り出す手続きが必要。memqは後ろしか取り出せないので自前で作るしかなさそう。まずはmemqとはある意味逆の動きをするrmemqを定義。これは先頭から指定の項の前までのリストを返す。

(define (rmemq i s)
  (if (or (null? s) (eq? (car s) i))
      null
      (cons (car s) (rmemq i (cdr s)))))

動かしてみると

> (rmemq '+ '(1 + 2))
(1)
> (rmemq '+ '(1 + 2 + 3))
(1)
> (rmemq '+ '((1 + 2) + 3))
((1 + 2))
>

やはりリストに入っているのが邪魔。同様に単項のリストであれば項のみを取り出す様にしてfirst-termを定義。the-restと同じ様にも定義出来るが、少し凝って式を引数に取る関数を返してみる。それを使ってaddend、multiplier、baseを定義する。

(define (first-term op)
  (lambda (s)
    (let ((f (rmemq op s)))
      (if (null? (cdr f))
          (car f)
          f))))

(define (addend s)
  ((first-term '+) s))
(define (multiplier p)
  ((first-term '*) p))
(define (base e)
  ((first-term '**) e))

動作確認

> (addend '(1 + 2 + 3))
1
> (addend '((1 + 2) + 3))
(1 + 2)
> (addend '(3 * x ** 2 + 4 * x))
(3 * x ** 2)
> (augend '(3 * x ** 2 + 4 * x))
(4 * x)
> (multiplier '(3 * x ** 2))
3
> (multiplicand '(3 * x ** 2))
(x ** 2)
> (base '(x ** 2))
x
> (exponent '(x ** 2))
2
>

述語達はそれぞれの演算子が式に含まれているかをチェックするだけなのでmemqが使える。memqは入れ子になったリストの中まではチェックしないので都合が良い。入れ子の場合はその部分だけを取り出して処理するまでは単一の項として評価される必要があるので。

(define (sum? x)
  (memq '+ x))
(define (product? x)
  (memq '* x))
(define (exponentiation? x)
  (memq '** x))

memqは指定の項がリストに含まれないとfalseで、含まれている場合にはリストが返るが、if等の条件分岐はfalseでなければtrueとして扱ってくれる。

取り敢えず教科書の例

> (deriv '(x + 3 * (x + y + 2)) 'x)
4
> (deriv '(x + 3 * (x + y + 2)) 'y)
3
>

もう少しややこしい奴

> (deriv '(4 * x ** 3 + x ** 2 * 3 + 5 * x + 6) 'x)
((4 * (3 * (x ** 2))) + (((2 * x) * 3) + 5))
>

make-sum、make-product辺りで上手く簡約は出来ていないが、取り敢えずは出来ていそう。