Simple Boltzmann Machine Simulation

bm.lisp

(in-package :cl-user)
(eval-when (:compile-toplevel :load-toplevel :execute)
  (ql:quickload '(:cl-num-utils :lla :array-operations) :silent t))
(defpackage :bm
  (:use :cl)
  (:import-from :clnu
                :e*
                :transpose)
  (:import-from :clnu.mx
                :mx)
  (:import-from :lla
                :mm)
  (:import-from :aops
                :flatten
                :reshape
                :reshape-col
                :nrow
                :ncol
                :sub)
  (:export :bm))
(in-package :bm)

(defparameter *node-count* 3)
(defparameter *ws* (clnu.mx:mx t
                     ( 0 -2 -2)
                     (-2  0 -2)
                     (-2 -2  0)))
(defparameter *bs* (clnu.mx:mx t
                     (-5)
                     (-5)
                     (-5)))
(defparameter *temperature* 1)
(defparameter *update-count* 100)

(defun sigmoid (s temperature)
  (/ 1 (+ 1 (exp (/ (- s) temperature)))))

(defun prob-sigmoid (s temperature)
  "シグモイド関数を使った確率的2値モデル。"
  (if (< (random 1d0) (sigmoid s temperature)) 1 0))

(defun append1 (lst x)
    (append lst (list x)))

(defun power-set (n)
  "n bitの状態変数の集合をリストとして生成する。"
  (labels ((rec (n acc)
             (if (> n 0)
                 (rec (1- n)
                      (append (mapcar #'(lambda (x) (append1 x 0)) acc)
                              (mapcar #'(lambda (x) (append1 x 1)) acc)))
                 acc)))
    (rec n (list nil))))

(defun power-set-array (n)
  "n bitの状態変数の集合をベクタとして生成する。"
  (mapcar #'(lambda (x) (aops:reshape-col (coerce x 'vector))) (power-set n)))

(defun energy (ws bs xs)
  "エネルギー関数。状態変数の存在確率の理論値を求めるために使用。"
  (+ (reduce #'+ (aops:flatten (clnu:e* -1/2 ws (lla:mm xs (clnu:transpose xs)))))
     (aref (lla:mm (clnu:transpose bs) xs) 0 0)))

(defun softmax (ws bs xs temperature)
  "重みws, 閾値bs, 温度temperatureにおける状態変数xsの存在確率を求める。"
  (loop for s in (power-set-array (aops:nrow xs))
     sum (exp (- (/ (energy ws bs s) temperature))) into total
     finally (return (/ (exp (- (/ (energy ws bs xs) temperature))) total))))

(defun weighted-sum (ws xs bs row)
  (- (lla:mm (aops:sub ws row) (aops:flatten xs)) (aref bs row 0)))

(defun update-nodes! (ws bs xs temperature)
  "状態変数xsを重みws, 閾値bsにより1要素ずつ順番に更新する。xsは変更されることに注意。"
  (loop for i below (aops:nrow xs)
     for s = (weighted-sum ws xs bs i)
     for x = (prob-sigmoid s temperature)
     do (setf (aref xs i 0) x)
     finally (return xs)))

(defun make-nodes (n)
  "0/1にランダムに初期化した状態変数を作成する。"
  (let ((xs (make-array (list n 1))))
    (dotimes (i n)
      (setf (aref xs i 0) (random 2)))
    xs))

(defun make-hist (n)
  (make-array (expt 2 n)))

(defun hist-count (hist bin)
  "ヒストグラムのbinをインクリメントする。"
  (incf (aref hist bin)))

(defun vector-to-number (xs &key (radix 10))
  "vectorをradix進数の数値として解釈する。状態変数をヒストグラムのbinにエンコードするために使用する。"
  (parse-integer (format nil "~{~a~}"
                         (let ((*print-base* radix))
                           (map 'list #'princ-to-string
                                (reverse (aops:flatten xs)))))
                 :radix radix))

(defun bm (ws bs &key (update-count *update-count*) (temperature *temperature*))
  "ボルツマンマシンをシミュレートして状態の存在確率の理論値と実測値を表示する。
エネルギー関数から求めた重みws, 閾値bsを与える。"
  (assert (= (aops:ncol ws) (aops:nrow bs)))
  (let* ((node-count (aops:nrow ws))
         (xs (make-nodes node-count))
         (hist (make-hist node-count)))
    (loop repeat update-count
       for bin = (vector-to-number (update-nodes! ws bs xs temperature) :radix 2)
       do (hist-count hist bin))
    (let ((expected (mapcar #'(lambda (xs) (softmax ws bs xs temperature)) (power-set-array 3)))
          (actual (map 'list #'(lambda (x) (float (/ x update-count))) hist)))
      (format t "理論値: ~{~1,3f~^ ~}~%" expected)
      (format t "実測値: ~{~1,3f~^ ~}~%" actual))))

#+ (or) (bm:bm (clnu.mx:mx t (0 -2 -2) (-2 0 -2) (-2 -2 0))
               (clnu.mx:mx t (-5) (-5) (-5))
               :update-count 1000
               :temperature 1)

#+ (or) (bm:bm (clnu.mx:mx t (0 -2 -2) (-2 0 -2) (-2 -2 0))
               (clnu.mx:mx t (-1) (-1) (-1))
               :update-count 1000
               :temperature 1)