Think Stitch
PRINCIPIA  最近の更新


非決定的計算オペレータ amb の並列化

Structure and Interpretation of Computer Programs [1] (計算機プログラムの構造と解釈,以下 SICP)の第4章に,非決定的計算のためのオペレータ amb というものが出てきます.On Lisp [2] の第22章には choose と fail という同じく非決定的計算のためのオペレータが出てきます.元は Lisp の創案者である John McCarthy さんによって論文 [3] で紹介されたものだそうです.このオペレータを並列化してみようというのがこの文書の主題です.

言語はもちろん Scheme, 処理系は Gauche を使います.

非決定的選択オペレータ amb による天使の選択

非決定的計算も amb も初めて聞いたという人のためにかんたんに説明します.より詳しいことは SICP を見てください.パズルなど面白い例が出ています.

非決定的選択オペレータ amb を使って (amb e1 e2 ...) という式を書くと,e1, e2, ... の中からどれか1つが非決定的に選ばれて評価され,その値が返ってきます.「非決定的」というくらいですから,どれが選ばれるかわからない,任意である,ということなのですけど,1つ条件があるのです.これを例で説明します.

いま,大きさ 10 の配列 v があるとして,こんなコードを書いてみます.やりたいことは,配列の要素が 0 であるようなインデックスを求めることです.必ずどこかに 0 があるとしておきます.

(let ((k (amb 0 1 2 3 4 5 6 7 8 9)))
  (let ((x (vector-ref v k)))
    (if (= x 0)
        k
        (amb))))

はじめに amb を使って配列のインデックス k を1つ選びます.そうしたら対応する要素を x に取り出します.もしそれが 0 であれば見つけたことになりますから,インデックス k を返します.そうでない場合は (amb)に出会います.

この (amb) は何かというと,引数が1つもないので「選択するものがない」という意味です.選択するものがないということは「ここに来てはいけない」ということを意味するのです.インデックスを選ぶ最初の amb が与えられた式の中からどれを選ぶのかは基本任意なのですが,1つ条件があって,それは選んだあとの計算で (amb) に出会わないようなものを選ぶのです.したがってこのコードは要素が 0 であるインデックスを返します.

なんとも不思議な話です.非決定的選択オペレータ amb は,まるで未来が見通せるかのように振る舞います.コード上では最初の amb がインデックスを選んだ後で条件判定がなされているのに,その結果を予知してるかのように (amb) に出会わない選択を行います.この不思議な選択を「天使の選択」と呼ぶ人もいます.

未来を予知する選択は1度しかできないわけではありません.次のようなコードを考えます.これはピタゴラスの定理を満たす3つの数,いわゆるピタゴラス数を求めることを意図したものです.

(let ((x (amb 2 3 4 5 6 7 8 9))
      (y (amb 2 3 4 5 6 7 8 9))
      (z (amb 2 3 4 5 6 7 8 9)))
  (if (= (+ (* x x) (* y y)) (* z z))
      (list x y z)
      (amb)))

選択を行う amb が3つあります.コード上ではそれぞれ独立して選択を行うように見えるのに,なぜか3つがからんだ条件 \(x^2 + y^2 = z^2\) が満たされるように選択が行われるのです.

amb が持つこの不思議な能力のおかげで,何かを見つける探索のコードがとても短く直感的に書けることがわかります.単に「これこれの中からこの条件を満たすものを選べ」と書くだけで済むわけですから.

非決定的選択オペレータ amb に初めて出会ったとき,変分原理のことを思い出しました.光が進む経路に関するフェルマーの原理というのをご存知でしょうか.フェルマーの最終定理で有名なフェルマーさんが見つけた原理で,光は進むのに必要な時間が最小になるような経路を進むというものです.まるで光は未来を予知できるといっているように聞こえます.これと同じような神秘さを amb に感じたものでした.しかし実際には魔法はなく,光にも amb も,そう見えるような「しくみ」が内部にあるわけです.

条件を満たす値がないときはどうなるのだろう,とか,複数あるときはどうなるのだろう,といったことが気になっている人もいるでしょう.これらについても「しくみ」を見るとわかります.

探索による非決定的選択オペレータ amb の実装

amb を実装する方法の1つはバックトラックによる探索です.処理系は amb に出会ったら,1つ選択を行うと同時に,後で戻ってこられるように場所と残りの選択肢を記録しておきます.そうしておいてから計算を続けて,もし (amb) に出会ったら,最後に行った選択ポイントまで戻って,別の選択肢を選び,再度試行をします.もし選択肢がなくなってしまった場合は,さらにその前の選択ポイントまで戻ります.もし選択ポイントがなくなってしまったら,それは解がないということなので,諦めて失敗します.

このしくみを基に考えてみると,まず解がない場合にどうなるかはわかりました.解が複数ある場合は,選択の順序に応じて最初に見つかったものが返ってくるということになります.したがって実は選択は決定的であるということになります.また,最初の説明のとおりに未来が見通せるのであれば,コードの長さに応じた時間で答えを見つけられそうでしたが,実際には試行錯誤をしているので残念ながらそれなりの時間がかかるということになります.なんだかちょっと恋が冷めた感じでしょうか.

それでも実装の詳細はとても魅力的なのです.魔法の種は魔法以上に面白いので,ぜひ SICP を見ていただきたいと思います.実際に動かせるコードを丁寧に説明してくれます.On Lisp では複数の実装方法が解説されています.

並列化による非決定的選択オペレータ amb の実装

バックトラックによる実装とは別の話として,SICP には次のようなパラグラフがあります.

抽象的にはamb式の評価は時を多くの分岐に分け, 各分岐では式の一つの可能な値を使って計算が続くと考えられる. ambは 非決定的な選択点(nondeterministic choice point)を表現するという. 動的に割り振られる十分な数の処理装置を持つ計算機があれば, 探索は直截的な方法で実装出来る. 実行はamb式に出会うまで逐次的な計算機のように進行する. amb式に出会うとより多くの処理装置が割り振られ, 選択を示す並列実行のすべてが継続するよう初期化される. 各処理装置は, それが唯一の選択であるかのように, 失敗に出会って終了するか, 更に分岐するか, 終了するかまで逐次的に進行する

つまり amb に出会ったら選択肢の数だけ処理装置を用意して,並列に計算すればいいわけです.これならば高速な探索が期待できます.そこでここではこの並列化を実際にやってみたいと思います.

スレッドプール

処理装置としてスレッドを使うことをすぐに思いつきます.しかし OS が提供するネイティブなスレッドは amb が対象とするような小さな計算を実行するためにはあまり向きません.計算そのもののコストに対して,スレッドの生成等々のコストが大きすぎるからです.加えて amb による分岐は莫大な数になる可能性がありますが,ネイティブスレッドはそのような使い方を想定していないからです.

そこで代わりにスレッドプールを使うことにします.スレッドプールという名前を初めて聞いたという人のためにかんたんに説明すると,膨大な数の小さな仕事があるとき,それぞれにスレッドを割り当てるのではなくて,比較的少数のスレッドを用意しておいて,それらで仕事を分け合いながら実行する手法のことをいいます.たいていの場合は,スレッドの数を実際に物理的に存在する処理装置(CPU あるいはコア)の数にしておくと効率よく計算できます.

終了処理

アイデアは単純なのですが,考えるべきことがあと3つほどあります.まず1つ目は,どんどん並列化していった結果,誰かが解にたどり着いたときにどうするのかということです.この場合は見つけた解を即座に返しつつ,残りの仕事をすべてキャンセルする必要があります.2つ目は,誰も解にたどり着かなかった場合はどうするかです.その場合はそれを検出して失敗で返すようにします.いずれの場合もすべてのスレッドを安全に終了させるためにひと工夫要ります.

並列化する仕事の切り出しと CPS 変換

3つ目は並列に実行する仕事をスレッドに分配できるように切り出すことです.これはちょっと説明を要します.たとえば

(amb e1 e2 e3)

という式に出会ったら,各式を thunk に変換して

(amb-par
  (lambda () e1)
  (lambda () e2)
  (lambda () e3))

とすれば,thunk を単位としてスレッドに分配できます.amb-par は分配を行うシステムコールのようなものと考えることができます.

しかしこれだけでは十分ではありません.最初の例をもう1度見てみます.

(let ((k (amb 0 1 2 3 4 5 6 7 8 9)))
  (let ((x (vector-ref v k)))
    (if (= x 0)
        k
        (amb))))

これを見ると,amb の引数だけでなく amb で選択が行われたあとに実行する処理も並列化する必要があることがわかります.たとえばこんなふうになります.

(amb-par
  (lambda ()
    (let ((k 0))
      (let ((x (vector-ref v k)))
        (if (= x 0)
            k
            (amb)))))
  (lambda ()
    (let ((k 1))
      (let ((x (vector-ref v k)))
        (if (= x 0)
            k
            (amb)))))
  (lambda ()
    (let ((k 2))
      (let ((x (vector-ref v k)))
        (if (= x 0)
            k
            (amb)))))
  ...)

これは冗長なので,共通部分を切り出すと次のようになるでしょう.

(let ((cont
        (lambda (k)
          (let ((x (vector-ref v k)))
            (if (= x 0)
                k
                (amb))))))
  (amb-par
    (lambda () (cont 0))
    (lambda () (cont 1))
    (lambda () (cont 2))
    ...))

この cont は何かといえば,amb 式の結果を待っている継続ということになります.amb に出会ったら,その選択肢である引数だけでなく,amb 式の継続を含めて thunk にする必要があるということです.そこで,継続が取り出せるようにするためにプログラムを継続渡しスタイル(Continuation-Passing Style, CPS)に変換します.CPS 変換の結果,式 (amb e1 e2 e3) の継続が k であったとすると

(amb-par
  (lambda () (k e1))
  (lambda () (k e2))
  (lambda () (k e3)))

とすればいいわけです.実際には e1, e2, e3 自身も CPS に変換する必要があります.なぜなら,その中にも amb が出現する可能性があるからです.

SICP の実装を読んだ人は,あれがちょっとかわった CPS 変換を含んでいることを知っているでしょう.バックトラックによる実装でも,選択ポイントに戻るためと,異なる選択について計算を再試行するために継続が必要なのでした.On Lisp ではこれを call/cc を使って実装した例が出ています.

実装

では,具体的な実装の話に入っていきたいと思います.

Scheme のコア言語

CPS 変換をかんたんに行うために,対象とする式を最低限のものに限定します.これをコア言語と呼ぶことにします. コア言語の要素は以下のとおりです.

コードを処理するプログラムを書くときは,この小さなコア言語だけを考えればよいところが Scheme の好きな点の1つです.

Scheme 本来の"フル"言語からコア言語への変換はかんたんなので省略します.勢いあまってマクロ展開子を作ってしまってもいいかもしれません :P

CPS 変換

CPS 変換のアルゴリズムはいろいろありますが,ここではわかりやすさを優先して,いちばんかんたんなものを使うことにします.冗長なコードが生成されますけど,それは別に簡約することにします.

(define (cps-convert expr cont)
  (match expr
    ((? symbol?)
     (if (eq? expr 'call/cc)
         `(,cont
           (lambda (cc f)
             (f cc (lambda (k v) (cc v)))))
         `(,cont ,expr)))
    (('quote liteal) `(,cont ,expr))
    (('set! var e)
     (let ((r (gensym)))
       (cps-convert
        e
        `(lambda (,r) (,cont (set! ,var ,r))))))
    (('if e1 e2 e3)
     (let ((r (gensym)))
       (cps-convert
        e1
        `(lambda (,r)
           (if ,r
               ,(cps-convert e2 cont)
               ,(cps-convert e3 cont))))))
    (('begin . es)
     (if (null? (cdr es))
         (cps-convert (car es) cont)
         (let ((r (gensym)))
           (cps-convert
            (car es)
            `(lambda (,r)
               ,(cps-convert `(begin ,@(cdr es)) cont))))))
    (('lambda formals e)
     (let ((k (gensym)))
       `(,cont (lambda ,(cons k formals)
                 ,(cps-convert e k)))))
    (('amb . es)
     (let ((k (gensym "k")))
       `((lambda (,k)
           (amb-par
            ,@(map (lambda (e)
                     `(lambda () ,(cps-convert e k)))
                   es)))
         ,cont)))
    ((f . args)
     (if (primitive-operator? f)  ; eq? car cdr cons ...
         (let ((vs (map (lambda (a) (gensym)) args)))
           (fold-left
             (lambda (a v e)
               (cps-convert e `(lambda (,v) ,a)))
             (list cont (cons f vs))
             vs args))
         (let ((vs (map (lambda (a) (gensym)) expr)))
           (fold-left
             (lambda (a v e)
               (cps-convert e `(lambda (,v) ,a)))
             (cons (car vs) (cons cont (cdr vs)))
             vs expr))))
    (_  ; constant
     `(,cont ,expr))))

ポイントは amb のところです.先に説明したように継続を含めて thunk に変換し,amb-par を呼び出す形にします.

余談ですが,call/cc が次のようなかんたんな関数になってしまうところが CPS 変換の面白いところです.call/cc で取得した"継続"(内側の lambda 式)を呼び出すと,そこでの継続 k を捨ててしまっていることが見て取れます.

(lambda (cc f)
  (f cc (lambda (k v) (cc v))))

簡約

上記の素朴な CPS 変換を適用すると ((lambda (x ...) e) a ...) という式がたくさんできてしまいます.これを簡約することにします.ちょっと手抜きですけど許してください.もちろん簡約しなくても動きます.

(define (beta-reduce expr)
  (match expr
    ((('lambda formals body) . args)
     (if (every atomic? args)  ; symbol? or constant? other than lamnda-expr
         (let ((body-u (subst (map cons formals args) body)))
           (if (equal? body body-u)
               body
               (beta-reduce body-u)))
         `((lambda ,formals ,(beta-reduce body))
           ,@(map beta-reduce args))))
    ((? symbol?) expr)
    (('lambda formals body)
     `(lambda ,formals ,(beta-reduce body)))
    ((f . args)
     (map beta-reduce expr))
    (_ expr)))
(define (subst s expr)
  (match expr
    ((? symbol?)
     (let ((p (assq expr s)))
       (if p (cdr p) expr)))
    (('lambda formals body)
     (let ((s-u (filter (lambda (p) (not (memq (car p) formals))) s)))
       `(lambda ,formals ,(subst s-u body))))
    (('set! var e)
     `(set! ,var ,(subst s e)))
    ((f . args)
     (map (lambda (e) (subst s e)) expr))
    (_ expr)))

スレッドプール

続いてスレッドプールを作ります.Gauche にはスレッドプールモジュールがありますけど,ここでは自分で作ってみることにします.並行プログラミングのよい練習になります.

スレッドプールの設計でポイントになるのは,どうやって仕事を分配するかというところです.ここでは work stealing という方式を採用することにします.

まず各スレッドには固有のキューを持たせます.各スレッドは自分のキューから仕事を取り出して順番に処理します.仕事は thunk で表すので呼び出すだけです.

呼び出しの結果起こりうることは2つあります.1つは解が見つかる場合で,そのときには関数 succeed を呼び出すと約束しておきます.つまり成功時に呼び出す継続です.succeed が呼び出されたら他のスレッドを安全に終了させたのち,発見した値を返します.

もう1つは amb-par が呼び出される場合です.これはさらに2つに場合分けできて,選択肢がない場合とある場合があります.選択肢がない場合というのは (amb) に出会ったということです.したがってこの仕事は終わりで,次の仕事をキューから取り出して実行します.選択肢がある場合はそれをすべて自分のキューに入れて,さらにどれかを選んで実行します.

問題は自分のキューが空になってしまった場合です.自分のキューが空でも,他のスレッドはまだ仕事を抱えているかもしれません.そこで他のスレッドから仕事を「盗む」ことにします.これが work stealing と呼ばれる所以です.乱数を生成してスレッドを1つ選び,そこから仕事を盗みます.もし選んだ相手のキューも空だった場合は,再び乱数を生成するところからやり直します.

もしすべてのスレッドのキューが空になったとしたら,それは解がなかったということですから全体を失敗で終了させる必要があります.これを検出するために,仕事を抱えているスレッドの数を数えておくことにします.これが 0 になったら全体を終了させます.

グローバル変数

では具体的なコードを見ていきます.まずはグローバル変数です.

(define mutex (make-mutex))
(define cv (make-condition-variable))
(define notfound (gensym "notfound"))
(define shutdown-flag #f)
(define result notfound)
(define qv #f)
(define rgv #f)
(define num-running 0)

排他制御のための mutex と,それからメインスレッドがスレッドプールの終了を待つための条件変数 cv を用意します.notfound は解が見つからなかったことを表すトークン(マーカー)です.shutdown-flag は解が見つかった場合と仕事がなくなった場合にスレッドプールを終了させるためのフラグです.result は見つけた解を入れる変数です.

qv は各スレッド用のキューを並べたベクタです.キューとしては Gauche の mtqueue を使います.複数のスレッドがキューにアクセスしますからスレッドセーフな版を使う必要があります.ふつう work stealing というと lock-free の deque を使うのですが,ここでは mtqueue で代用します.一般にはスレッドが自分のキューから仕事を取り出したり追加したりするときは前から,他のスレッドのキューから仕事を盗むときは後ろからするのですが,今回はすべて前から操作します.ですのでキューといいながら実際はスタックです.この構造は探索の順序,あるいは戦略に関係します.探索の戦略というと深さ優先探索と幅優先探索がよく使われます.これから調べる予定のリスト,いわゆるフロントをキューにすると幅優先,スタックにすると深さ優先になります.探索を並列化すると,やや乱暴にいえば幅優先傾向が出ると思われます.問題と戦略には相性があるので,一概にどちらがいいとはいえませんが,work stealing の立場では,スレッドが操作する端点が異なっている方が衝突が少ないのでよいといわれています.

rgv は乱数ジェネレータのベクタです.排他制御を避けるためにスレッドの数だけ用意しておきます.num-running は仕事を抱えているスレッドの数です.

起動関数

つづいて起動の関数です.まず各グローバル変数を初期化して,id=0 のスレッドのキューに最初の仕事を入れます. それからスレッドを起動します.スレッド固有領域には id を入れておきます.あとは仕事が終わるまで条件変数を待って,起こされたらスレッドの終了を待ち,結果を引き取って返ります.

(define (init-amb-thread-pool n thunk)
  (mutex-lock! mutex)
  (set! result notfound)
  (set! shutdown-flag #f)
  (set! num-running n)
  (set! qv (list->vector (map (lambda (_) (make-mtqueue)) (iota n))))
  (set! rgv (list->vector (map (lambda (_) (integers$ n)) (iota n))))
  (queue-push! (vector-ref qv 0) thunk)
  (let ((thread-list
         (map (lambda (id)
                (let ((th (make-thread amb-thread)))
                  (thread-specific-set! th id)
                  (thread-start! th)
                  th))
              (iota n))))
    (mutex-unlock! mutex cv)
    (for-each thread-join! thread-list)
    result))

スレッドの処理

スレッドのエントリポイントは amb-thread です.プログラムが CPS ですからスレッドの処理はすべて tail call でつないでいきます.まず shutdown-flag を検査して,#t の場合は終了します.そうでなければ自分のキューから thunk を取り出して実行します.キューが空の場合は work-steal-begin へ行きます.

(define (amb-thread)
  (unless shutdown-flag
    (let ((id (thread-specific (current-thread))))
      (let ((q (vector-ref qv id)))
        (let ((thunk (queue-pop! q notfound)))
          (if (eq? thunk notfound)
              (work-steal-begin (vector-ref rgv id))
              (thunk)))))))

キューが空になったら他のスレッドから仕事を盗むのですが,その前に空ではないスレッドの数を更新しておきます.これは終了判定のためです.数の更新には排他制御が必要なので mutex をロックして行います.もしすべてのスレッドが無職になったら shutdown-flag#t にし,条件変数を signal してスレッドプールを終了します.そうでなければ work-steal に行きます.

(define (work-steal-begin rg)
  (mutex-lock! mutex)
  (dec! num-running)
  (if (= num-running 0)
      (begin
        (set! shutdown-flag #t)
        (condition-variable-signal! cv)
        (mutex-unlock! mutex))          ; exit
      (begin
        (mutex-unlock! mutex)
        (work-steal rg))))

乱数を生成してスレッドを選び,仕事を盗みます.もし成功したらスレッドの数を更新してから実行します.盗めなかった場合は乱数生成から再試行します.まぬけにも自分のキューを覗くことがありますけど,気にしないことにします.

(define (work-steal rg)
  (thread-yield!)
  (unless shutdown-flag
    (let ((k (rg)))
      (let ((q (vector-ref qv k)))
        (let ((thunk (queue-pop! q notfound)))
          (if (eq? thunk notfound)
              (work-steal rg)
              (begin
                (mutex-lock! mutex)
                (inc! num-running)
                (mutex-unlock! mutex)
                (thunk))))))))

つづいては amb-par です.まず選択肢,つまり thunk のリストが空であるかどうかを調べます.空の場合,この仕事は終わりですから,就職活動のためエントリポイント amb-thread へ行きます.

選択肢がある場合は自分のキューに格納して,それから1つ選んで実行します.選択肢が1つしかない場合はちょっとだけ効率化します.

(define (amb-par . thunk-list)
  (unless shutdown-flag
    (cond ((null? thunk-list)
           (amb-thread))
          ((null? (cdr thunk-list))
           ((car thunk-list)))
          (else
           (let ((id (thread-specific (current-thread))))
             (let ((q (vector-ref qv id)))
               (dolist (thunk (cdr thunk-list))
                 (queue-push! q thunk))
               ((car thunk-list))))))))

最後は成功の継続 succeed です.結果を result に格納し,shutdown-flag を立て,条件変数を signal してメインスレッドを起こしてから終了します.

(define (succeed v)
  (mutex-lock! mutex)
  (set! result v)
  (set! shutdown-flag #t)
  (condition-variable-signal! cv)
  (mutex-unlock! mutex))

以上で実装は終わりです.

実行例

いくつかの例を実行してみます.

並列化されていることを確認するための実行例

並列化されていることを確認するために,次のようなコードを動かしてみました.各選択が行われるときに,sys-sleep で待ち時間が発生するようにしてあります.スレッドの数を1にすると深さ優先探索になるので11秒かかります.しかしスレッド数を4にすると4秒で終わります.

(time
 (init-amb-thread-pool
  4
  (lambda ()
    (amb-par
     (lambda ()
       (sys-sleep 1)
       (amb-par
        (lambda () (sys-sleep 3) (amb-par))
        (lambda () (sys-sleep 2) (amb-par))))
     (lambda ()
       (sys-sleep 2)
       (amb-par
        (lambda () (sys-sleep 1) (amb-par))
        (lambda () (sys-sleep 2) (succeed #t))))))))

N-Queens 問題

Lisp で探索問題の定番といえば N-Queens 問題かなと思いました.N-Queens は解がたくさんあるタイプの問題で,探索する空間の中にそれらが散らばっています.

一般に探索問題では探索の順序によって解にたどり着くまでのステップ数に違いが出ます.プログラムの実行に非決定性があると,実行のたびに早く見つかったり時間がかかったりするわけです.

故意に失敗させた場合

そこでまず,並列化の性能を見るために,かならず探索に失敗するようにしてみました.具体的には,各列に順番に Queen を配置していくのですが,最後の列だけどこにおいても必ず conflict するという設定にしてみます.解なしにしているだけでなく,探索空間を広げている感じです.プログラムはすべての可能性を調べたのちに失敗で終了します.

(define (req p)
  (if (not p) (amb)))

(define (an-element-of items)
  (req (not (null? items)))
  (amb (car items) (an-element-of (cdr items))))

(define (conflict? m k qs) (= (length qs) (- m 1)))

(define (q m n qs)
  (if (= n m)
      qs
      (let ((k (an-element-of (iota m))))
        (if (conflict? m k qs)
            (amb)
            (q m (+ n 1) (cons k qs))))))

このコードで全探索にかかる実行時間を計測してみました.測定環境は Gauche 0.9.4, Debian 8 amd64, Core i7 4770 3.7GHz (4core/8HT) です.

N=7 の場合

まず Queen の数 N=7 の場合です.濃い青い線が計測値です.結構早い段階で飽和しています.workers=8 で少し落ちています.

Amdahl の法則

黄色とオレンジ色の線は Amdahl の法則による曲線です.Amdahl の法則というのは,処理装置の数 \(N\) を増やしていったときにどれだけ速度が向上するのかを表した法則です.

大まかに考えて,プログラムの全実行時間を \(1\) としたとき,比率 \(p\ \ (0 \le p \le 1)\) の部分が並列に実行できて,残りの部分は排他処理等で逐次的にしか実行できないとします.このプログラムを \(N\) 個の処理装置で実行すると \(p\) の部分の実行時間は理論的には \(p/N\) になります.逐次的にしか実行できない部分 \(1-p\) はそのままです.したがって全体の実行時間は \((1-p) + p/N\) となり,その逆数を求めればどれだけ速くなったかがわかります.

\[ Speedup = \frac{1}{(1-p)+p/N} \]

この法則にしたがっていくつかの \(p\) についてグラフを描いてみると次のようになります.

見てのとおり,この法則はけっこう恐ろしい法則で,プログラムの 90% もの部分が並列に実行可能でも,10 倍以上には速くならないのですね.

やや乱暴ですが,この法則を N-Queens のグラフに当てはめてみると,\(p\) は 0.8 前後という感じでしょうか.したがって 20% くらいの部分がロック等で逐次的にしか実行できていない部分という感じです.あくまでもおおざっぱな話です.実際,Amdahl の法則による曲線と実測値とはカーブが異なるので,何か別の要因があるのだと思われます.

N=8 の場合

続いて N=8 の場合です.6スレッドのところに何か落ち込みがありますが,原因はわかりません.

N-Queens: 通常の探索

では本来の問題で探索をやってみます.衝突判定関数を以下のものと入れ替えます.

(define (conflict? m k qs)
  (if (memv k qs)
      #t
      (let loop ((a (+ k 1)) (b (- k 1)) (qs qs))
        (if (null? qs)
            #f
            (let ((q (car qs)))
              (if (or (eqv? a q)
                      (eqv? b q))
                  #t
                  (loop (+ a 1) (- b 1) (cdr qs))))))))

Queen の数 N=16 の場合です.各場合とも1,000回測定した平均です.6スレッドからかなり落ちてしまいます.

いちばん速かった5スレッドの場合の度数分布です.かなり分散が大きくなっています.この非決定性は並列化から来ているわけです.スレッドが実行される順序により探索の順序が変わるということです.暴れまくりですね.かなり早く解を見つける場合がある一方で,遅い場合は5倍くらいかかっています.

つづいて N=20 の場合です.えらく暴れています.もう法則がどうこういうレベルではありません.Gauche を1度起動したインスタンスの中で連続して測定したので,メモリ管理回りの影響があるのかもしれないと後で思いました.1回ごとにフレッシュな実行をすべきでした.たまたまこうなったわけではなく何回かの計測で同じ傾向が見られました.

4 スレッドで6倍はおかしいと思われるかもしれません.実は問題として N=20 の場合は深さ優先だと発見までのパスが長い方なのです.ですから探索の順序が変わるとそういうことは起こり得るとは思います.

各スレッド数毎に分布を見ていきます.

worker = 1.これはシングルスレッドで深さ優先ですから安定した分布です.

worker = 2.なぜか4か所に山があります.

3.4sec あたりを細かく分けてみました.

worker = 3.山が2つあります.

前の部分を細かく.

worker = 4.やはり山が2つあります.実行中に GC 関係で何かあって,山が2つになるのかと推測しましたが,確認はできていません.

worker = 5

worker = 6

worker = 7

worker = 8

まとめ

非決定的計算オペレータ amb を並列化してみました.並列化による高速化の効果が少しだけありました.もっと高速化するには mtqueue を deque で置き換えるとか,ロックを減らすなどのくふうが必要だと思われます.また並行性の方からくる非決定性によるばらつきがかなり大きいことがわかりました.他にも説明のできていない現象も課題として残りました.

メモ

雑多に思いつくことを書いておきたいと思います.

この並列化実装には SICP の実装と大きく振る舞いが異なる部分が1つあります.それは代入の扱いです.SICP の実装では,バックトラックの際に代入の効果もロールバックするという実にかっこいいことをやっていますが,ここでの並列実装にはありません.局所的な束縛は同じ効果がありますが,グローバル変数の変更は SICP でいうところの permanent-set!になり,他のルートに影響します.

代入と関連して,探索空間に循環や重複がある場合は訪問済みノードを記録する必要があります.これは排他制御を必要とするので結構な足かせになります.

たまに amb が深さ優先であることを前提としたコードがありますが,当然のことながら動きません.たとえばすべての解を列挙する場合などです.そのための並列化というのはまた面白そうではあります.On Lisp では列挙を効率よく行うためのカットというしくみが紹介されています.これを並列化された計算でも実現することはできるでしょうか.

N-Queens 問題は深さ優先が適した問題なので,幅優先が適した問題でやってみると何か別の面白いことがあるかもしれません.加えてスレッドプールの分配戦略を変えると楽しそうです.

プログラムの意味を状態空間 \(\Sigma\) から状態空間 \(\Sigma\) への関係 \(R \subseteq \Sigma \times \Sigma\) で表したとすると,部分的な関係というものを考えることができます.関係の定義域が状態空間のすべてを網羅していない,つまり真部分集合になっている場合です.もし関係 \(R\) で表されるあるコードが,それに続くコードの意味 \(S\) の定義域に入っていないような状態を生成した場合,計算を続けることができません.これは (amb) に相当します.2つのコードの連接を単純に関係の合成

\[ R; S = \{(x, y) | \exists z. (x, z) \in R \land (z, y) \in S\} \]

で表したとすると,この存在限量のおかげで \(S\) の定義域に入ってこない \(R\) の計算は取り除かれてしまいます.つまり関係の合成は未来を見通す「天使の選択」なのですね.事後条件に存在限量が入っているときはふつう探索になりますから amb との関連が見えます.プログラマの仕事の1つの見方は,非決定性と天使の選択を使って書かれた仕様を,より決定的で予測可能で効率の良いアルゴリズムに変換することですが,amb は非決定的で予測不可能で非効率でもよければ,仕様をそのまま動かしてくれるという感じです.Prolog にも似ています.どちらも探索ですからあたりまえかもしれませんが.

「こうしたらこんなに速くなった!」ということがあったら教えてください.

次は,きっと誰かが「amb を GPU で実装してみた」とか「FPGA で実装してみた」という記事を書いてくれると期待しています :P

参考文献

[1] Structure and Interpretation of Computer Programs,
   https://mitpress.mit.edu/sicp/
  (邦訳)計算機プログラムの構造と解釈,
   http://sicp.iijlab.net/

[2] On Lisp,
   http://www.paulgraham.com/onlisp.html

[3] A Basis For a Mathematical Theory of Computation, John McCarthy, 1961-1963,
   http://ropas.snu.ac.kr/~kwang/4190.310/mccarthy63basis.pdf

2015/12/05
© 2013,2014,2015 PRINCIPIA Limited