お絵描きロジック（ピクロス）を最適化で解く [MATLAB]

この記事は 数理最適化 Advent Calendar 2024 の4日目の記事です。

はじめに

Switchで昔懐かしいテレビゲームが遊べるということで、子供と『マリオのスーパーピクロス』をやっていたところ、ふと「・・・ん？これ最適化問題じゃね？（俺ソルバやらされてる…）」と気づき、数理最適化で解いてみようと思いました。

先行の取り組みはいくつか見つかった（§ 参考リンク を参照）のですが、全く同じアプローチというのも面白くないので、ここでは整数表現・バイナリ表現を併用する定式化で攻めてみます。また 宗教 職業上の理由で MATLAB の『Optimization Toolbox』を使います。

ちなみに『ピクロス』という言葉は任天堂の登録商標とのことで、一般的には『お絵描きロジック』などと呼ばれるようです。ここでは後者で呼ぶことにします。

ルール

『お絵描きロジック』は、各行／各列の隣に並んだ数字のヒントを基に、黒マス／白マスで表現されるドット絵を復元するゲームです。以下の任天堂の公式動画を観るとなんとなくイメージがつかめるでしょうか。

各ヒント数字は「黒マスブロックの長さ」を表現しています。例えば 幅10マス の行においてヒントが【 1 2 3 】であった場合、左から長さ1、長さ2、長さ3の黒ブロックが並んでいることを示しています。つまり、その行は以下の 9パターン のどれかとなります。ここで黒マスは ■ で、白マスは ○ で表しています。

1 2 3 : [ ■○■■○■■■○○ ]
1 2 3 : [ ■○■■○○■■■○ ]
1 2 3 : [ ■○■■○○○■■■ ]
1 2 3 : [ ■○○■■○■■■○ ]
1 2 3 : [ ■○○■■○○■■■ ]
1 2 3 : [ ■○○○■■○■■■ ]
1 2 3 : [ ○■○■■○○■■■ ]
1 2 3 : [ ○■○○■■○■■■ ]
1 2 3 : [ ○○■○■■○■■■ ]

このように各黒ブロックについて、長さは事前にわかっているのですが、どこに位置するのかは通常わかりません。しかしヒントによってはパターンによらず黒／白が確定するマス目もあります。実際上記の例のどのパターンにおいても、右から３つ目 は常に ■ です。こういったところを突破口に順次確定するマスを埋めていく…というのが人間がとる基本戦略となります。特にヒント数字がより大きな値であるほど、より多くの確定マスが生じます。このような思考を各行／各列について繰り返し考え、全体として矛盾のないドットパターンを見つけるというのが『お絵描きロジック』というゲームです。

実を言うと、与えられたヒント数字によってはドットパターンが一意に決まらない場合もありえます。例えば 2x2 の盤面での市松模様（オセロの初期配置）はそれにあたります。これでは絵を復元するゲームとしては成立しないため、ここでは解は一意であると仮定して話を進めます。…ちなみに『マリオのスーパーピクロス』では解が一意でない問題が実は含まれているらしいです。私は全クリしましたが全く気づきませんでした。

問題の定式化

ここからが話の本題です。この『お絵描きロジック』を最適化問題として定式化しましょう。各マスの状態はバイナリで表現され、ヒント数字は制約条件を構成し、最適化によって具体的に求めたい値（MATLAB流に「最適化変数」と呼びます）は各黒マスブロックの開始位置となります。ゲームの目的は制約条件を満たすパターン（＝実行可能解）を見つけることに他なりません。これを MATLAB (Optimization Toolbox) を使って解きます。

MATLAB では最適化問題を次のように定式化します。今回は正当なクイズであれば実行可能解は一意なはずですので、複数の解の良し悪しを比べる指標である目的関数は設けません（適当な定数として考えます）。ここでは最小化問題として構成しますが、最大化問題としても結果は同じです。

% 最小化問題として定式化をはじめる
optprob = optimproblem(ObjectiveSense="minimize");

MATLAB では「問題ベースの最適化」と「ソルバーベースの最適化」という２種類のアプローチが提供されています。今回は前者を使います。問題を抽象的に記述できる点や、問題の構成から適切なソルバーが自動的に選択される点が便利です。

問題の記述

ヒント数字を次のような配列に格納しておきます。左詰めで数字を並べ、空欄は 0 とします。ここでは前述の任天堂公式動画の問題をセットしました。

hintH = [  % (M, Kh)
	5;
	2;
	2;
	2;
	2;
];
hintV = [  % (N, Kv)
	1, 0;
	1, 1;
	1, 3;
	3, 1;
	2, 0;
];

盤面サイズを $(M,N)$、行／列のヒント数字の個数を $K_h$, $K_v$ とそれぞれ表します。またそれらの添字をそれぞれ $n (=1,\ldots,N)$, $m (=1,\ldots,M)$, $k_h (=1,\ldots,K_h)$, $k_v (=1,\ldots,K_v)$ と書きます。

% 盤面のサイズ
M = height(hintH);
N = height(hintV);
% 行・列の最大ブロック数
Kh = width(hintH);
Kv = width(hintV);

各行／各列のヒント数字をあらためて数式としても表します。

• $H_h(m, k_h)$： 第 $m$ 行における $k_h$ 番目の黒マスブロックの長さです。
• $H_v(n, k_v)$： 第 $n$ 列における $k_v$ 番目の黒マスブロックの長さです。

最適化変数（整数表現）の定義

最適化によって求めたい具体的な値は「各行／各列の各黒マスブロックの開始位置」です。これを整数の最適化変数として定義します。

• $P_h(m, k_h)$： 第 $m$ 行における第 $k_h$ 黒ブロックの開始位置（整数）です。
• $P_v(n, k_v)$： 第 $n$ 列における第 $k_v$ 黒ブロックの開始位置（整数）です。

これをコードに落とすと以下となります。

% ブロックの開始位置（整数表現）
Ph = optimvar("Ph", M, Kh, Type="integer", LowerBound=1, UpperBound=N);
Pv = optimvar("Pv", N, Kv, Type="integer", LowerBound=1, UpperBound=M);

最適化変数（バイナリ表現）の定義

上記の整数表現に加えて バイナリ表現 についても考えます。バイナリ変数を導入すると論理計算などを線形演算として記述でき、線形の枠組みでありながら比較的幅広い問題をカバーできます。実際にはそれぞれで表現の得手・不得手はありますが、今回はどちらも使いたいので上記の整数表現に同期するバイナリ表現も一緒に用意します。

• $Q_h(m, k_h, n)$： 第 $m$ 行における第 $k_h$ 黒ブロックの開始位置が $n$ マス目であるか？（バイナリ）
• $Q_v(n, k_v, m)$： 第 $n$ 列における第 $k_v$ 黒ブロックの開始位置が $m$ マス目であるか？（バイナリ）

% ブロックの開始位置（バイナリ表現）
Qh = optimvar("Qh", M, Kh, N, Type="integer", LowerBound=0, UpperBound=1);
Qv = optimvar("Qv", N, Kv, M, Type="integer", LowerBound=0, UpperBound=1);

各黒ブロックの開始位置はどこか一か所のマスに限られるため、それを制約条件として記述します。これはいわゆる整数値の One-hot表現 となっています。

• $\sum_{n=1}^{N} Q_h(m, k_h, n) = 1$, $(\forall m, k_h)$： 各黒ブロックの開始位置はどこか一箇所だけ
• $\sum_{m=1}^{M} Q_v(n, k_v, m) = 1$, $(\forall n, k_v)$： 各黒ブロックの開始位置はどこか一箇所だけ

% One-hot表現とする制約
optprob.Constraints.QhOnehot = (sum(Qh, 3) == 1);  % One-hot表現
optprob.Constraints.QvOnehot = (sum(Qv, 3) == 1);  % One-hot表現

また、上記の整数表現とバイナリ表現が互いに整合する（矛盾しない）ように制約条件を設けます。

• $(\forall m, k_h)$, $P_h(m, k_h) = \sum_{n=1}^{N} n Q_h(m, k_h, n)$： どちらの表現でも同じ開始位置を表す
• $(\forall n, k_v)$, $P_v(n, k_v) = \sum_{m=1}^{M} m Q_v(n, k_v, m)$： どちらの表現でも同じ開始位置を表す

% 整数値とOne-hot表現を整合させる
optprob.Constraints.PhQh = (Ph(:) == reshape(Qh, [], N) * (1:N)');
optprob.Constraints.PvQv = (Pv(:) == reshape(Qv, [], M) * (1:M)');

整数表現とバイナリ表現の整合の取り方はこのほかにも様々あります。気になる方は参考リンクの『整数計画法による定式化入門（著：藤江哲也）』をご覧ください。

制約：黒ブロックは必ず１マス以上空けて並べること

続いてヒント数字に基づく制約条件を設けます。これは整数表現を用いると簡潔に記述できます。

ある黒ブロックの一つ右の黒ブロックは必ず１マス以上空けて並べる必要があります。

• $P_h(m, k_h) + H_h(m, k_h) + 1 \leq P_h(m, k_h+1)$, $(\forall m, k_h = 1,2,\ldots,K_h-1)$
• $P_v(n, k_v) + H_v(n, k_v) + 1 \leq P_v(n, k_v+1)$, $(\forall n, k_v = 1,2,\ldots,K_v-1)$

また最後の黒ブロックは盤外に出ないように収める必要があります。

• $P_h(m, K_h) + H_h(m, K_h) - 1 \leq N$, $(\forall m)$
• $P_v(n, K_v) + H_v(n, K_v) - 1 \leq M$, $(\forall n)$

これらをコードに落とし込みます。

% 隣り合うブロック間での開始位置に関する制約
optconsPositionH = optimconstr(M, Kh);
optconsPositionV = optimconstr(N, Kv);
for y = 1 : M  % 各行について
	K = nnz(hintH(y, :));  % その行でのブロック数
	if K == 0
		continue
	end
	optconsPositionH(y, 1:K-1) = (Ph(y, 1:K-1) + hintH(y, 1:K-1) + 1 <= Ph(y, 2:K));
	optconsPositionH(y, K) = (Ph(y, K) + hintH(y, K) - 1 <= N);  % その行での最後のブロック
end
for x = 1 : M  % 各列について
	K = nnz(hintV(x, :));  % その列でのブロック数
	if K == 0
		continue
	end
	optconsPositionV(x, 1:K-1) = (Pv(x, 1:K-1) + hintV(x, 1:K-1) + 1 <= Pv(x, 2:K));
	optconsPositionV(x, K) = (Pv(x, K) + hintV(x, K) - 1 <= M);  % その列での最後のブロック
end
optprob.Constraints.PositionH = optconsPositionH;
optprob.Constraints.PositionV = optconsPositionV;

制約：各行／各列から求めた２つのパターンが厳密に一致すること

各マスが黒／白かを判定するための前準備として、以下のバイナリ変数を導入します。

• $B_h(m, n, K_h)$： マス $(m,n)$ がその行の $k_h$ 番目の黒マスブロックと重なっているか？
• $B_v(m, n, K_v)$： マス $(m,n)$ がその列の $k_v$ 番目の黒マスブロックと重なっているか？

% マス(m,n)が、その行/列の第kh/kvブロック上か？
Bh = optimvar("Bh", M, N, Kh, Type="integer", LowerBound=0, UpperBound=1);
Bv = optimvar("Bv", M, N, Kv, Type="integer", LowerBound=0, UpperBound=1);

マス $(m,n)$ がその行／列におけるいずれかの黒ブロックと重なっていれば True となります。この判定を制約条件で記述します。この際、バイナリ表現を導入したことで論理和を加算で表現できています。

• $B_h(m, n, k_h) = \sum_{b = \max(1, n - H_h(m, k_h))}^{n} Q_h(m, k_h, b)$, $(\forall m, n, k_h)$
• $B_v(m, n, k_v) = \sum_{b = \max(1, m - H_v(n, k_v))}^{m} Q_v(n, k_v, b)$, $(\forall m, n, k_v)$

constrH = optimconstr(M, N, Kh);
constrV = optimconstr(M, N, Kv);
for y = 1 : M
	for x = 1 : N
		for kh = 1 : Kh
			onblock = (max(1, x-hintH(y, kh)+1) : x);
			constrH(y, x, kh) = (Bh(y, x, kh) == sum(Qh(y, kh, onblock), 3));
		end
	end
end
for x = 1 : N
	for y = 1 : M
		for kv = 1 : Kv
			onblock = (max(1, y-hintV(x, kv)+1) : y);
			constrV(y, x, kv) = (Bv(y, x, kv) == sum(Qv(x, kv, onblock), 3));
		end
	end
end

最後に、列から求めるドットパターンと行から求めるドットパターンが互いに厳密に一致しなければならないという制約を設けます。

• $\sum_{k_h = 1}^{K_h} B_h(m, n, k_h) = \sum_{k_v = 1}^{K_v} B_v(m, n, k_v)$, $(\forall m, n)$

optprob.Constraints.BoardConsistency = (sum(Bh, 3) == sum(Bv, 3));

ソルバで解く

さて、これで一連の定式化が完了しました。まずは定式化の結果を確認します。複雑な数式がズラズラと表示されます。

show(optprob);

次のこれをソルバに解いてもらいましょう。ここでの定式化の結果は 《混合整数線形計画問題》 となっており、効率よく解けるタイプの問題クラスであることが知られています。MATLAB では intlinprog というソルバが適用されるはずです。

sol = solve(optprob);

ソルバ実行後、返値の sol に最適化変数の求解結果が格納されます。この結果から答えのドットパターンを復元できます。実際に解いたところ以下のパターンが得られました。カタカナの《マ》です。正しく答えが出せてますね。

[ Answer ]
-- 1 1 1 3 2
-- 0 1 3 1 0
5 ■■■■■
2 ○○○■■
2 ○○■■○
2 ○■■○○
2 ○○■■○

ちなみに 30x30 程度のサイズになりますと、最適化問題の内部構成に時間がかかるようですが、ソルバの計算自体はほぼ一瞬で終わりました。

おわりに

『お絵描きロジック』を数理最適化で解いてみました。最適化の特に難しい点（と同時に醍醐味）は「対象となる問題をどう定式化するか？」にあると思いますが、どういった形で変数定義・制約記述すればよいのか、資料が少なく数学も高度であるために、なかなか慣れないものです。本記事がそういった学習の参考になれば幸いです。

コード全体

% 盤面のサイズ
M = height(hintH);
N = height(hintV);
fprintf("盤面のサイズ:  (M, N) = (%d, %d)\n", M, N);

% 行・列の最大ブロック数
Kh = width(hintH);
Kv = width(hintV);
fprintf("最大ブロック数:  (Kh, Kv) = (%d, %d)\n", Kh, Kv);

%% 最適化問題の定式化

optprob = optimproblem(ObjectiveSense="minimize");

% ブロックの開始位置（整数表現）
Ph = optimvar("Ph", M, Kh, Type="integer", LowerBound=1, UpperBound=N);
Pv = optimvar("Pv", N, Kv, Type="integer", LowerBound=1, UpperBound=M);

% 隣り合うブロック間での開始位置に関する制約
optconsPositionH = optimconstr(M, Kh);
optconsPositionV = optimconstr(N, Kv);
for y = 1 : M  % 各行について
	K = nnz(hintH(y, :));  % その行でのブロック数
	if K == 0
		continue
	end
	optconsPositionH(y, 1:K-1) = (Ph(y, 1:K-1) + hintH(y, 1:K-1) + 1 <= Ph(y, 2:K));
	optconsPositionH(y, K) = (Ph(y, K) + hintH(y, K) - 1 <= N);  % その行での最後のブロック
end
for x = 1 : M  % 各列について
	K = nnz(hintV(x, :));  % その列でのブロック数
	if K == 0
		continue
	end
	optconsPositionV(x, 1:K-1) = (Pv(x, 1:K-1) + hintV(x, 1:K-1) + 1 <= Pv(x, 2:K));
	optconsPositionV(x, K) = (Pv(x, K) + hintV(x, K) - 1 <= M);  % その列での最後のブロック
end
optprob.Constraints.PositionH = optconsPositionH;
optprob.Constraints.PositionV = optconsPositionV;

% ブロックの開始位置（バイナリ表現）
Qh = optimvar("Qh", M, Kh, N, Type="integer", LowerBound=0, UpperBound=1);
Qv = optimvar("Qv", N, Kv, M, Type="integer", LowerBound=0, UpperBound=1);

% 整数値をOne-hot表現するための制約
optprob.Constraints.QhOnehot = (sum(Qh, 3) == 1);  % One-hot表現
optprob.Constraints.QvOnehot = (sum(Qv, 3) == 1);  % One-hot表現
optprob.Constraints.PhQh = (Ph(:) == reshape(Qh, [], N) * (1:N)');  % One-hot表現を整数値と紐づけ
optprob.Constraints.PvQv = (Pv(:) == reshape(Qv, [], M) * (1:M)');  % One-hot表現を整数値と紐づけ

% (y,x)マスが、その行/列の第kh/kvブロック上か？
Bh = optimvar("Bh", M, N, Kh, Type="integer", LowerBound=0, UpperBound=1);
Bv = optimvar("Bv", M, N, Kv, Type="integer", LowerBound=0, UpperBound=1);
constrH = optimconstr(M, N, Kh);
constrV = optimconstr(M, N, Kv);
% Bh = optimvar("Bh", M, N, Type="integer", LowerBound=0, UpperBound=1);
% Bv = optimvar("Bv", M, N, Type="integer", LowerBound=0, UpperBound=1);
% constrH = optimconstr(M, N);
% constrV = optimconstr(M, N);
for y = 1 : M
	for x = 1 : N
		for kh = 1 : Kh
			onblock = (max(1, x-hintH(y, kh)+1) : x);
			constrH(y, x, kh) = (Bh(y, x, kh) == sum(Qh(y, kh, onblock), 3));
		end
	end
end
for x = 1 : N
	for y = 1 : M
		for kv = 1 : Kv
			onblock = (max(1, y-hintV(x, kv)+1) : y);
			constrV(y, x, kv) = (Bv(y, x, kv) == sum(Qv(x, kv, onblock), 3));
		end
	end
end
optprob.Constraints.constrH = constrH;
optprob.Constraints.constrV = constrV;

% 各行・各列のブロック配置による盤面が両者一致すること
optprob.Constraints.BoardConsistency = (sum(Bh, 3) == sum(Bv, 3));

% 最適化問題を確認
show(optprob);

% 最適化問題の求解
sol = solve(optprob);
pv = sol.Pv;
ph = sol.Ph;

if isempty(ph) || isempty(pv)
	warning("解が見つかりませんでした！");
	return
end

% 盤面を復元
boardH = false(M, N);
for y = 1 : M
	for kh = 1 : Kh
		boardH(y, ph(y, kh) + (0:hintH(y, kh)-1)) = true;
	end
end
boardV = false(M, N);
for x = 1 : N
	for kv = 1 : Kv
		boardV(pv(x, kv) + (0:hintV(x, kv)-1), x) = true;
	end
end

% 両者の整合性の確認
if ~all(boardH == boardV, "all")
	error("両者の盤面が一致していません。");
end

fprintf("[ Answer ]\n");
for y = 1 : M
	for x = 1 : N
		if boardV(y, x)
			fprintf("■");
		else
			fprintf("・");
		end
	end
	fprintf("\n");
end

参考リンク