ニューロンに数独を解かせてみよう
電通大生による電通大生のためのAdvent Calendar 2024 3日目
UECアドカレ3日目です。2日目は我らがSHINNさんの「ε論法を論理学からみつめてみよう！」でした。
離散数学で苦しんでいる身としてはとてもありがたい記事でした。全人類読もう！
はじめに
UEC24のbubuzukeです。気づいたら入学してから8ヶ月経過していて戦慄しています。年をとると時間の流れが早く感じますね...。ノリと勢いでAdvent Clandarに登録してしまったので(2敗)泣きながら記事を書いています。せっかくなので技術系の話をしようと思います。(間違ってることを書いてたら教えてください修正して切腹します...)
内容
今回の記事では今流行りのニューロモルフィック計算について少し触れたいと思います。このテーマを選んだ理由は締切前日から書き始めて間に合いそうだったから 提案者であるホップフィールドさんが今年のノーベル物理学賞を受賞して話題的にアツかったからです。
ニューロンとは
まず今回の記事でもっとも重要になってくるニューロンの性質について軽く触れようと思います。我々の脳はニューロンという神経細胞が相互に結合することで情報を処理、記憶しています。ニューロンは他のニューロンから入力をもらうと、自分の中で電気を加算し、それがあるレベルを超えるとスパイクと呼ばれる約100ｍｖの短い電気パルスを出力します。一つ一つのニューロンの挙動自体は簡単でも数百億スケールで複雑なネットワークを構築することで知覚や運動など非常に複雑な活動が可能になっています。
ノイマン型コンピュータと脳の共通点
ついでに同じ情報処理装置であるコンピュータと脳の共通点についても触れたいと思います。コンピュータと脳の基本素子は半導体素子とニューロンですが、そのどちらも電気パルスを用いて情報を伝達しています。意外(？)なことに脳はアナログ値ではなく電気パルスによる1か0かで情報を伝達しています。
組み合わせ最適化問題の近似解法
そろそろ本題に入りたいと思います。組み合わせ最適化問題とは多くの解の中から最もよいものを見つける問題のことを言います。有名なものだと巡回セールスマン問題があります。組合せ最適化問題には解を見つけるのは困難でも、与えられた候補が解であるかを判別するのは簡単という性質を持つものがあります。このような場合、解に非常に近い候補を見つけるだけでも意味がある場合が多く、そのような解法は近似解法と呼ばれます。この近似解法の汎用アルゴリズムとして先ほど説明した複数ニューロンのネットワーク(以下スパイキングネットワーク)が用いられることがあります。今回はスパイキングネットワークを用いて数独を解いていこうと思います。
やってみよう
では実際に数独を解くためのスパイキングネットワークを組んでみましょう。今回は簡単のために4×4のマスに0~3の数字を当てはめるケースを考えます。まず、各マスに0~3に対応するニューロンを割り当てます。今回の場合、ニューロンの数は4×4×4=64個となります。次に数独のルールをニューロンの結合パターンで表現します。まず最初に「各マスには0~3のどれか一つの数字しか入らない」という条件を考えます。これは各マスごとに4個のニューロンが互いに活動を抑制するように結合することで表現できます(下図)。
同様に各行、列、ブロックについても互いに抑制性結合をすることで数独のルールを表現することができます。
最初から入る数字がわかっている場合は、外部電流を与えることにします。これで数独のルールをスパーキングネットワークで表現することができました。すべての結合が抑制性というところが重要です。また、電流にランダムなノイズを加えることでネットワークの状態が数独の解に収束するようになります。
実装してみる
以下に実装例を示します。
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <unistd.h>
#include <stdint.h>
#include <stdbool.h>
#include <SFMT.h>
#define NN (2)
#define N ((NN)*(NN)) // size of a board
int32_t puzzle[N][N] = {
{0, -1, -1, -1},
{2, -1, 0, 1},
{3, 2, -1, 0},
{1, -1, -1, 2}
};
#define N_NEURON ((N)*(N)*(N))
#define TAU (10.)
#define W_INH (-10.)
#define E_REST (-65.)
#define THETA (-55.)
#define I_EXT (20.)
#define DECAY (0.5)
#define NT (10000)
#define DT (1.0)
static inline int32_t id(const int32_t i, const int32_t j, const int32_t a) {
return (a + N * (j + N * i));
}
static inline int32_t wid(const int32_t post, const int32_t pre) {
return (pre + N_NEURON * post);
}
void createConnection(float *w) {
for(int32_t i = 0; i < N; i++) {
for(int32_t j = 0; j < N; j++) {
for(int32_t a = 0; a < N; a++) {
int32_t post = id(i, j, a);
// grid
for(int32_t b = 0; b < N; b++) {
int32_t pre = id(i, j, b);
int32_t x = wid(post, pre);
w[x] = (a != b) ? W_INH : 0.;
}
// row
for(int32_t l = 0; l < N; l++) {
for(int32_t b = 0; b < N; b++) {
int32_t pre = id(l, j, b);
int32_t x = wid(post, pre);
w[x] = (i != l && a == b) ? W_INH : 0.;
}
}
// column
for(int32_t l = 0; l < N; l++) {
for(int32_t b = 0; b < N; b++) {
int32_t pre = id(i, l, b);
int32_t x = wid(post, pre);
w[x] = (j != l && a == b) ? W_INH : 0.;
}
}
// box
for(int32_t ll = 0; ll < NN; ll++) {
for(int32_t mm = 0; mm < NN; mm++) {
int32_t ib = i / NN, jb = j / NN;
int32_t l = ll + NN * ib, m = mm + NN * jb;
for(int32_t b = 0; b < N; b++) {
int32_t pre = id(l, m, b);
int32_t x = wid(post, pre);
w[x] = (!(i == l && j == m) && a == b) ? W_INH : 0.;
}
}
}
}
}
}
}
void setInput(float *i_ext) {
for(int32_t i = 0; i < N; i++) {
for(int32_t j = 0; j < N; j++) {
int32_t b = puzzle[i][j];
if(b >= 0) {
for(int32_t a = 0; a < N; a++) {
int32_t x = id(i, j, a);
i_ext[x] = (a == b) ? 1 : 0;
}
}
}
}
}
void printAnswer(const int32_t *ns) {
for(int32_t i = 0; i < N; i++) {
for(int32_t j = 0; j < N; j++) {
int32_t maxa = 0, maxs = ns[id(i, j, maxa)];
for(int32_t a = 0; a < N; a++) {
if(maxs < ns[id(i, j, a)]) {
maxa = a;
maxs = ns[id(i, j, maxa)];
}
}
fprintf(stderr, "%d%s", maxa, (j == N - 1) ? "\n" : " ");
}
}
}
int main(void) {
float *w = (float *)calloc(N_NEURON * N_NEURON, sizeof(float));
createConnection(w);
float *i_ext = (float *)calloc(N_NEURON, sizeof(float));
setInput(i_ext);
float *v = (float *)calloc(N_NEURON, sizeof(float));
for(int32_t i = 0; i < N_NEURON; i++) v[i] = E_REST;
float *g_syn = (float *)calloc(N_NEURON, sizeof(float));
bool *s = (bool *)calloc(N_NEURON, sizeof(bool));
int32_t *ns = (int32_t *)calloc(N_NEURON, sizeof(int32_t));
sfmt_t rng;
sfmt_init_gen_rand(&rng, getpid());
for(int32_t nt = 0; nt < NT; nt++) {
for(int32_t i = 0; i < N_NEURON; i++) {
float r = 0;
for(int32_t j = 0; j < N_NEURON; j++) {
int32_t x = wid(i, j);
r += w[x] * s[j];
}
g_syn[i] = DECAY * g_syn[i] + r;
}
for(int32_t i = 0; i < N_NEURON; i++) {
float dv = DT * (-(v[i] - E_REST) + g_syn[i] + I_EXT * (i_ext[i] + sfmt_genrand_real2(&rng))) / TAU;
s[i] = (v[i] > THETA);
v[i] = s[i] ? E_REST : (v[i] + dv);
if(nt >= NT / 2) ns[i] += s[i];
}
}
printAnswer(ns);
free(v);
free(s);
free(i_ext);
free(w);
}
実行例も示します。
bubuzuke sudoku % make
gcc -O3 -std=gnu11 -Wall -I../misc/SFMT-src-1.5.1 -DSFMT_MEXP=19937 -c sudoku.c
gcc -O3 -std=gnu11 -Wall -I../misc/SFMT-src-1.5.1 -DSFMT_MEXP=19937 -o sudoku sudoku.o SFMT.o -lm
bubuzuke sudoku % ./sudoku
0 1 2 3
2 3 0 1
3 2 1 0
1 0 3 2
いい感じに解けてますね(雑)
一応コードの内容も少し解説します。
定数
#define NN (2)
#define N ((NN)*(NN))
int32_t puzzle[N][N] = {
{0, -1, -1, -1},
{2, -1, 0, 1},
{3, 2, -1, 0},
{1, -1, -1, 2}
};
#define N_NEURON ((N)*(N)*(N))
#define TAU (10.)
#define W_INH (-10.)
#define E_REST (-65.)
#define THETA (-55.)
#define I_EXT (20.)
#define DECAY (0.5)
#define NT (10000)
#define DT (1.0)
ここで数独のサイズや計算に必要な定数を定義しています。
DTは時間の刻み幅で、ニューロンのシミュレーションを行う際、微分方程式を数値的に解くために使用します。DTを小さくすることでより詳しくニューロン単体の活動を調べられますが今回はその必要がないので大きめの値になっています。
ネットワークの定義
void createConnection(float *w) {
for(int32_t i = 0; i < N; i++) {
for(int32_t j = 0; j < N; j++) {
for(int32_t a = 0; a < N; a++) {
int32_t post = id(i, j, a);
// grid
for(int32_t b = 0; b < N; b++) {
int32_t pre = id(i, j, b);
int32_t x = wid(post, pre);
w[x] = (a != b) ? W_INH : 0.;
}
// row
for(int32_t l = 0; l < N; l++) {
for(int32_t b = 0; b < N; b++) {
int32_t pre = id(l, j, b);
int32_t x = wid(post, pre);
w[x] = (i != l && a == b) ? W_INH : 0.;
}
}
// column
for(int32_t l = 0; l < N; l++) {
for(int32_t b = 0; b < N; b++) {
int32_t pre = id(i, l, b);
int32_t x = wid(post, pre);
w[x] = (j != l && a == b) ? W_INH : 0.;
}
}
// box
for(int32_t ll = 0; ll < NN; ll++) {
for(int32_t mm = 0; mm < NN; mm++) {
int32_t ib = i / NN, jb = j / NN;
int32_t l = ll + NN * ib, m = mm + NN * jb;
for(int32_t b = 0; b < N; b++) {
int32_t pre = id(l, m, b);
int32_t x = wid(post, pre);
w[x] = (!(i == l && j == m) && a == b) ? W_INH : 0.;
}
}
}
}
}
}
}
ここでニューロン同士の結合を定義しています。
外部入力
void setInput(float *i_ext) {
for(int32_t i = 0; i < N; i++) {
for(int32_t j = 0; j < N; j++) {
int32_t b = puzzle[i][j];
if(b >= 0) {
for(int32_t a = 0; a < N; a++) {
int32_t x = id(i, j, a);
i_ext[x] = (a == b) ? 1 : 0;
}
}
}
}
}
数独の初期状態に基づいて外部入力を設定しています。
時間ステップごとにニューロンの活動を更新
for(int32_t nt = 0; nt < NT; nt++) {
for(int32_t i = 0; i < N_NEURON; i++) {
float r = 0;
for(int32_t j = 0; j < N_NEURON; j++) {
int32_t x = wid(i, j);
r += w[x] * s[j];
}
g_syn[i] = DECAY * g_syn[i] + r;
}
for(int32_t i = 0; i < N_NEURON; i++) {
float dv = DT * (-(v[i] - E_REST) + g_syn[i] + I_EXT * (i_ext[i] + sfmt_genrand_real2(&rng))) / TAU;
s[i] = (v[i] > THETA);
v[i] = s[i] ? E_REST : (v[i] + dv);
if(nt >= NT / 2) ns[i] += s[i];
}
}
ここで他ニューロンからの入力を加算し、膜電位の値を更新しています。
膜電位が一定値を超えた場合スパイクを発射し、結合している他のニューロンを抑制します。
おまけ
一年生はRubyのほうが馴染みがあると思うのでRuby版も載せておきます(ちゃんと動くか確認できてない)。
$NN = 2
$N = $NN * $NN
$PUZZLE = [
[0, -1, -1, -1],
[2, -1, 0, 1],
[3, 2, -1, 0],
[1, -1, -1, 2]
]
$N_NEURON = $N * $N * $N
$TAU = 10.0
$W_INH = -1.0
$E_REST = -65.0
$THETA = -55.0
$I_EXT = 20.0
$DECAY = 0.5
$NT = 2000
$DT = 1.0
def id(i, j, a)
return (a + $N * (j + $N * i));
end
def wid(post, pre)
return pre + $N_NEURON * post
end
def createCoNNection(w)
$N.times do |i|
$N.times do |j|
$N.times do |a|
post = id(i, j, a)
#grid
$N.times do |b|
pre = id(i, j, b)
w[wid(post, pre)] = (a != b) ? $W_INH : 0.0
end
#row
$N.times do |l|
$N.times do |b|
pre = id(l, j, b)
w[wid(post, pre)] = (i != l && a == b) ? $W_INH : 0.0
end
end
#column
$N.times do |l|
$N.times do |b|
pre = id(i, l, b)
w[x = wid(post, pre)] = (j != l && a == b) ? $W_INH : 0.0
end
end
#box
$NN.times do |ll|
$NN.times do |mm|
ib = i / $NN
jb = j / $NN
l = ll + $NN * ib
m = mm + $NN * jb
$N.times do |b|
pre = id(l, m, b)
w[x = wid(post, pre)] = (!(i == l && j == m) && a == b) ? $W_INH : 0.0
end
end
end
end
end
end
end
def setInput(i_ext)
$N.times do |i|
$N.times do |j|
b = $PUZZLE[i][j]
if(b >= 0)
$N.times do |a|
i_ext[id(i, j, a)] = (a == b) ? 1: 0
end
end
end
end
end
def printAnswer(ns)
$N.times do |i|
$N.times do |j|
maxa = 0
maxs = ns[id(i, j, maxa)]
$N.times do |a|
if(maxs < ns[id(i, j, a)])
maxa = a
maxs = ns[id(i, j, maxa)]
end
end
print "#{maxa}#{j == $N - 1 ? "\n" : " "}"
end
end
end
def main
w = Array.new($N_NEURON * $N_NEURON, 0.0)
createCoNNection(w)
i_ext = Array.new($N_NEURON, 0.0)
setInput(i_ext)
v = Array.new($N_NEURON, $E_REST)
g_syn = Array.new($N_NEURON, 0.0)
s = Array.new($N_NEURON, false)
ns = Array.new($N_NEURON, 0)
srand(0)
$NT.times do |nt|
$N_NEURON.times do |i|
r = 0.0
$N_NEURON.times do |j|
r += w[wid(i, j)] * (s[j] ? 1 : 0)
end
g_syn[i] = $DECAY * g_syn[i] + r
end
$N_NEURON.times do |i|
dv = $DT * (-(v[i] - $E_REST) + g_syn[i] + $I_EXT * (i_ext[i] + rand(0.0..1.0))) / $TAU
s[i] = (v[i] > $THETA)
v[i] = s[i] ? $E_REST : (v[i] + dv)
ns[i] += 1 if s[i] && nt >= $NT / 2
#printf("%d %d\n", nt, i) if s[i]
end
end
printAnswer(ns)
end
end
まとめ
今回は4×4の簡単なケースについて試しました。ただこれはあくまで近似解法であるため、より大きな9×9の数独や初期状態での空白が多い場合では正しい解を出力できない場合が多いです。しかし、先にもいったように組み合わせ最適化問題においては近似解を見つけること自体に意味がある場合が多いため全くの無駄というわけではありません。今回は数独を解きましたが実は巡回セールスマン問題もこの近似解法で解くことができます。暇な人はぜひやってみてください。
最後に
最後まで読んでいただきありがとうございます。後半に行くについて説明が雑になってしまって読みづらかったと思います(反省)。
今回軽く触れた微分方程式の数値的解法についてはteam某11のアドカレで書けたらいいなと思っています。
参考文献
Discussion