２次割り当て問題の解を評価する(1)

２次割り当て問題(QAP)は巡回セールスマン問題と似たようなNP困難な問題である。QAPで検索すれば出てくるけど、大雑把にいうと、物資の移動量(流量)×距離がシステム全体で最小になる配置を探す問題だ。答えが分かっている問題例はQAPLIB Home Pageで得られるので、ここのデータを使えば解を探すアルゴリズムを試すことが出来る。

配布データは地点間の距離と、要素間の流量がそれぞれ行列で表されている。解（位置と要素の組合せ）に対するコスト値を求めるには、要素数×要素数の計算をする必要があるので、それなりに処理が重いのだ。lispで書いたのがこんな感じである*1。

(defvar *size*)
(defvar *distance*)
(defvar *flow*)

(defun qap-read (file)
  (with-open-file (fp file :direction :input
                      :if-does-not-exist :error)
    (let ((size (read fp))
          ll ll2)
      (dotimes (a size (setq ll (reverse ll)))
        (let (l)
          (dotimes (b size (push (reverse l) ll))
            (push (read fp) l))))
      (dotimes (a size (setq ll2 (reverse ll2)))
        (let (l)
          (dotimes (b size (push (reverse l) ll2))
            (push (read fp) l))))
      (values size ll ll2))))

(defun qap-load (file)
  (multiple-value-bind (size distance flow)
      (qap-read file)
    (setq *size* size
          *distance* distance
          *flow* flow)))

(defun %qap-eval (x)
  (macrolet ((distance (i j)
               `(elt (elt *distance* ,i) ,j))
             (flow (i j)
               `(elt (elt *flow* ,i) ,j)))
    (let ((sum 0))
      (dotimes (i *size* sum)
        (dotimes (j *size*)
          (incf sum (* (distance i j) (flow (elt x i) (elt x j)))))))))

;読み込み
(qap-load "tai30b.dat")

;解を作る
(let (l)
  (dotimes (a 30 (setq x (shuffle l)))
    (push a l)))

(compile '%qap-eval)
(time (dotimes (a 10000) (%qap-eval x)))
|
12390

10000回で12390msだから、1回に1ms以上かかっている。この問題の場合は解が30!通りあるので、もちろん10000個の解程度では探しきれない。xyzzy上では難しそうだ。といっても、他の処理系でも高速化の仕方がわからない。

(defun %qap-eval (x)
  (declare (type (list fixnum) x))
  (macrolet ((distance (i j)
               `(the fixnum (elt (elt *distance* ,i) ,j)))
             (flow (i j)
               `(the fixnum (elt (elt *flow* ,i) ,j))))
    (let ((sum 0))
      (declare (type fixnum sum))
      (dotimes (i *size* sum)
        (dotimes (j *size*)
          (incf sum (the fixnum
                         (* (distance i j)
                            (flow (elt x i)
                                  (elt x j))))))))))

一応これでGCL上で試したところ2550msにはなったが、まだ遅い。もっと速くする方法を教えて欲しい・・・。

*1:shuffleは過去の記事を参照。