💥

量子アニーラと数理計画ソルバーでシュレッダーに立ち向かう

2020/12/01に公開

数理最適化

量子コンピュータ

tech

この記事は BrainPad Advent Calendar 2020 初日の記事です。

突然ですが皆さん、 半沢直樹（新シリーズ） はご覧になりましたか？私は見ました。前作も見ていなかったので、あわせて見ました。
さて、半沢直樹の第三話にこんなシーンがありました。

証券取引等監視委員会の黒崎氏がシュレッダーで裁断された書類を見つけて、その場で復元してしまう話です。裁断された書類の山を持ち帰ろうとする部下と「ここでおやんなさいよ」と言い放つ黒崎氏。

河野大臣のブログ危機に直面する霞が関でも指摘されている通り、今若手の国家公務員離れが問題になっています。いや、全然くわしくないのですが。日経新聞でも転職希望の公務員が急増　外資やITへ流れる20代と報道されていました。公務員は我々国民を支えてくれている大事な仕事です。彼らの労働環境を改善しなければ、日本の未来は危ういと言えるでしょう。

というわけで、今回はシュレッダーで裁断された書類の復元を題材に、量子アニーラーや数理計画ソルバーで遊んでみたいと思います。検証コードは以下を御覧ください。

DARPA Shredder Challenge 2011

本取り組みにあたって少し調べたところ、2011年に面白い取り組みが行われていました。米国国防総省の研究機関である DARPA (Defense Advanced Research Projects Agency)] の開催した DARPA Shredder Challenge 2011というコンペティションです。複数のチームが、裁断された書類の復元に取り組みました。

DARPA Shredder Challenge 2011 の問題例

こちらのニュース記事で Shredder Challenge の結果が紹介されていますが、1989年のベルリンの壁の崩壊の際に東ドイツが破棄した書類の復元が未だに完了していない（2012年時点）という話や、クラウドソーシングを利用して人海戦術での復元を試みたり、そこに妨害工作員を忍び込ませたり、という話が載っていて面白かったです。

裁断された書類の復元と巡回セールスマン問題

さて、話を戻しましょう。裁断された書類を自動で復元するにはどうしたらよいでしょう？
人手の場合、とりあえず繋がりそうな紙片を見つけ出してはつないでいく作業を繰り返すと思います。繋がりそうかどうかは、1つの紙片の右端ともう一つの紙片の左端が似ているかどうかで判断します。

似ているかどうかを類似度と呼ぶことにすれば、「隣り合う紙片の類似度が最大となる紙片の並べ方を見つけなさい」という問題と捉えることができます。また、どれくらい似ていないかを非類似度と呼べば、「隣り合う紙片の非類似度が最小となる紙片の並べ方を見つけなさい」という問題とみなせます。これはいわゆる巡回セールスマン問題です。

巡回セールスマン問題は数理最適化の題材としてわかりやすいため、10月に開催されたGDG DevFest 2020でもお話させていただきました。興味のある方はその時の講演資料をご覧いただければと思いますが、DevFest で取り扱った問題は、説明のために色々とシンプルにしており、実際の裁断書類の復元とはいくつか違いがありました。例えば

実際の紙はシュレッダーによって裁断された切り口と紙のもともとの端とで滑らかさが違うため、周囲から埋めていく戦法が使える
紙片同士の類似度の計算に、切り口の形を考慮するなどの工夫が必要
実際には複数の書類をまとめてシュレッダーにかける事が多いので、問題は更に難しくなる
最近のシュレッダーは縦方向の祭壇だけでなく、横方向の裁断も同時に行う

などです。特に 4. についてはわかっていたのですが、裁断のされ方が違うと巡回セールスマン問題でなくなってしまうので、 DevFest では知らないふりをしていました（すいません）。そこで、今回は横方向も裁断された場合を取り上げます。

問題設定

少し長くなりましたが、こういうことです。

素朴な定式化

まずは素朴に定式化してみましょう。シュレッダーの復元問題は

$N$ 枚の紙片を縦 $R$ 、横 $C$ のマスに割当てる
各紙片との隣接する紙片との類似度の合計を最大化する

と捉えれば良さそうです。

この問題の目的関数は、紙片 $a$ の横（下）に紙片 $b$ が割り当てられた場合の類似度を表す定数 $C_{ab}$ （ $R_{ab}$ ）と、紙片 $a$ を位置 $(i, j)$ に割当てる場合のみ $1$ となる0-1変数 $x_{a,i,j}$ を使って、

\mathrm{Maximize}\quad \sum_{abij}\left\{C_{ab}x_{a,i,j}x_{b,i+1,j} + R_{ab}x_{a,i,j}x_{a,i,j+1}\right\}\tag{1}

と書けます。 $x_{a,i,j}x_{b,i+1,j}$ は、紙片 $a$ が位置 $(i,j)$ にあって、かつ、紙片 $b$ が位置 $(i + 1, j)$ にあった場合のみ $1$ となるので、第1項目が横方向の類似度の合計となります。同様に同様に第2項目が縦方向の類似度の合計を表します。

制約はどうなるでしょうか。位置 $(i, j)$ に複数の紙片を置くことはできませんし、紙片 $a$ を複数箇所に置くこともできません。これは

\begin{aligned}\tag{2} s.t.\quad \sum_{a} x_{a,i,j} &= 1, \\ \sum_{i,j} x_{a,i,j} &= 1. \end{aligned}

とすればよさそうです。では、この問題をいくつかの方法で解いてみましょう。

（混合）整数線形計画法による解法

上記の目的関数は $x_{a,i,j}x_{b,i+1,j}$ という2次の項があり、一見非線形のように見えます。しかし、 $x\in\{0, 1\}$ であることを使うと、実は線形な形で書き直す事ができます。２つの変数 $x_{a,i,j}, x_{b,k,l}$ と $x_{a,i,j}x_{b,k,l}$ の関係を見てみましょう

$x_{a,i,j}$	$x_{b,k,l}$	$x_{a,i,j}x_{b,k,l}$
0	0	0
0	1	0
1	0	0
1	1	1

上表から、

x_{a,i,j}x_{b,k,l} = \min(x_{a,i,j}, x_{b,k,l})

であることがわかります。そのため、前述の目的関数は新たな変数 $y_{a,b,i,j,k,l}\in\{0,1\}$ を使って、

\mathrm{Maximize}\quad \sum_{abij}\left\{C_{ab}y_{a,b,i,j,i+1,j} + R_{ab}y_{a,b, i,j,i,j+1}\right\}

\begin{aligned} s.t.\quad y_{a,b,i,j,k,l} &\leq x_{a,i,j},\\ y_{a,b,i,j,k,l} &\leq x_{b,k,l}. \end{aligned}

と書き直せます^[1]。もともとの $x_{a,b,i,j}$ に対する制約も合わせると、結局以下の問題を解けば良いことがわかります。

\mathrm{Maximize}\quad \sum_{abij}\left\{C_{ab}y_{a,b,i,j,i+1,j} + R_{ab}y_{a,b, i,j,i,j+1}\right\}

\begin{aligned} s.t.\quad \sum_{a} x_{a,i,j} &= 1, \\ \sum_{i,j} x_{a,i,j} &= 1, \\ \quad y_{a,b,i,j,k,l} &\leq x_{a,i,j},\\ y_{a,b,i,j,k,l} &\leq x_{b,k,l}. \end{aligned}

PuLP による実装と結果

実際に PuLP を用いて実装してみました。 PuLPは混合整数線形計画問題 (MIP; Mixed Integer Programming) のモデリング用のPythonライブラリーで、目的関数や制約式を直感的に記述することができます。また、PuLPと一緒に Cbc (Coin-or branch and cut) というソルバーもインストールされるので、お手軽に混合整数線形計画問題を解くことができます。

MILP定式化のPuLPによる実装コード

後で出てくる量子アニーリングと比較するため、目的関数の符号を反転させて最小化問題にしています。

def build_model(data):
    model = pulp.LpProblem(sense=pulp.LpMinimize)
    x = np.array(pulp.LpVariable.matrix(
        'x',
        (range(data.size), range(data.rows), range(data.cols)),
        cat=pulp.LpBinary
    ))
    y = {}

    costs = []
    for a in range(data.size):
        img_a = data.images[a]
        for b in range(data.size):
            if a == b:
                continue
            img_b = data.images[b]
            sim = (sim_y(img_a, img_b), sim_x(img_a, img_b))
            for row in range(data.rows):
                for col in range(data.cols):
                    if row < data.rows - 1:
                        y[a, b, row, col, row + 1, col] = pulp.LpVariable(
                            f'y_{a}_{b}_{row}_{col}_{row + 1}_{col}',
                            cat=pulp.LpBinary
                        )
                        costs.append(-sim[0]*y[a, b, row, col, row + 1, col])
                    if col < data.cols - 1:
                        y[a, b, row, col, row, col + 1] = pulp.LpVariable(
                            f'y_{a}_{b}_{row}_{col}_{row}_{col + 1}',
                            cat=pulp.LpBinary
                        )
                        costs.append(-sim[1]*y[a, b, row, col, row, col + 1])
    model.setObjective(pulp.lpSum(costs))


    for a in range(data.size):
        for b in range(data.size):
            if a == b:
                continue
            for row in range(data.rows):
                for col in range(data.cols):
                    if row < data.rows - 1:
                        model.addConstraint(
                            x[a, row, col] >= y[a, b, row, col, row + 1, col]
                        )
                        model.addConstraint(
                            x[b, row + 1, col] >= y[a, b, row, col, row + 1, col]
                        )
                    if col < data.cols - 1:
                        model.addConstraint(
                            x[a, row, col] >= y[a, b, row, col, row, col + 1]
                        )
                        model.addConstraint(
                            x[b, row, col + 1] >= y[a, b, row, col, row, col + 1]
                        )

    for a in range(data.size):
        model.addConstraint(
            x[a].sum() == 1
        )

    for row in range(data.rows):
        for col in range(data.cols):
            model.addConstraint(
                x[:, row, col].sum() == 1
            )

    return model, x

結果は以下のとおりでした。2x10くらいまでは割とすぐに答えが出せるようです。一方で 4x10 ではおかしな結果が得られてしまいました^[2]。

MIPによる定式化の結果

行数	列数	求解時間(sec)	最適値
2	3	0.047	-2080
2	5	26.9	-3070
4	5	270	-4126
2	10	354	-6340
4	10	2074	-3218

量子アニーラ（量子コンピューター）による解法

さて、話を少し戻して素朴な定式化 (1)-(2) をもう一度見てみると、これは量子アニーリングで解けそうだ、とすぐに分かるでしょう。

D-Wave に代表される量子アニーラは、量子効果を利用して QUBO (Quadratic Unconstrainted Binary Optimization) と呼ばれる組合せ最適化問題を解く専用機です。「Unconstrainted」とある通り、制約（今回の場合は (2)式）があると解けないので、少し問題を変形します。具体的には、適当に大きな定数 $M$ を用意して、以下のように書き換えます。

\begin{aligned} \mathrm{Minimize}\quad &-\sum_{abij}\left\{C_{ab}x_{a,i,j}x_{b,i+1,j} + R_{ab}x_{a,i,j}x_{a,i,j+1}\right\} \\ &+ M\left(\sum_{a} x_{a,i,j} - 1\right)^2 \\ &+ M\left(\sum_{i,j} x_{a,i,j} - 1\right)^2 \tag{3} \end{aligned}

まず、目的関数 (1) の符号を反転して、最大化ではなく最小化問題にします。次に制約式 (2) を満たすと自動的に 0 となるような項を追加します。 $M$ が大きな定数なので、目的関数 (3) を最小化するには、まず第2項と第3項が 0 である必要があり、結果として制約 (2) が満たされます。このあたりの式変形や量子アニーリングの基礎については昔ブログを書いているので、興味のある方は御覧ください。

上のブログを書いた2017年は、まだまだ量子アニーリングが普及しておらず量子アニーリングを利用するには高額な契約が必要でしたが、時代は変わりました。今年の8月に AWS がリリースした Amazon Braket を利用すると、誰でも手軽に初期費用無しで量子アニーラーを使うことができます。

PyQUBO を用いた実装と結果

QUBOの実装には、 PyQUBO を使うのがとても便利です。

PyQUBOを用いたモデル構築

def build_model(data, c=10):
    x = pyqubo.Array.create('x', shape=(data.size, data.rows, data.cols), vartype='BINARY')

    costs = []
    sim_max = {'x': 0, 'y': 0}
    for a in range(data.size):
        for b in range(data.size):
            if a == b:
                continue
            sim = {
                'x': sim_x(data.images[a], data.images[b]),
                'y': sim_y(data.images[a], data.images[b])
            }
            sim_max['x'] = max(sim_max['x'], sim['x'])
            sim_max['y'] = max(sim_max['y'], sim['y'])
            for row in range(data.rows):
                for col in range(data.cols):
                    if col < data.cols - 1:
                        costs.append(-sim['x']*x[a, row, col]*x[b, row, col + 1])
                    if row < data.rows - 1:
                        costs.append(-sim['y']*x[a, row, col]*x[b, row + 1, col])


    coeff = c*(sim_max['x'] + sim_max['y'])
    for a in range(data.size):
        costs.append(coeff*pyqubo.Constraint((np.sum(x[a]) - 1)**2, f'img_{a}'))
    for row in range(data.rows):
        for col in range(data.cols):
            costs.append(coeff*pyqubo.Constraint((np.sum(x[:, row, col]) - 1)**2, f'pos_{row}_{col}'))

    return sum(costs).compile(), x

大きな定数 $M$ を決める必要はあるものの、余分な変数 $y$ を導入する必要がないため、PuLP よりも少しだけシンプルなコードになっています。

では、D-Waveの実機を使う前にシミュレーション（シミュレーテッドアニーリング）で求解してみましょう。アニーリングは必ずしも最適解を見つけられるアルゴリズムではないので、何回も繰り返して実行して最も良かったものを採用するのが通常の使い方です。今回は1000回実行することとしました。以下のコードでは、 sampler.sample で結果を取得したあとで、最適なものを選びだしています。

シミュレーテッドアニーリングによる求解

import neal

sampler = neal.SimulatedAnnealingSampler()
sampleset = sampler.sample(bqm, num_reads=1000)
decoded_samples = model.decode_sampleset(sampleset)
best_sample = min(decoded_samples, key=lambda x: x.energy)

結果は以下のとおりです。MIPの結果と比較すると、求解時間は短いものの2x3以外では最適値にたどり着いていない事がわかります。シミュレーテッドアニーリングをそのまま利用するのは難しそうで、問題ごとにパラメータの調整の必要がありそうです^[3]。

シミュレーテッドアニーリングの結果

行数	列数	求解時間(sec)	目的関数	制約充足
2	3	1.12	-2080	true
2	5	2.7	-1489	true
4	5	12.6	-1224	true
2	10	11.1	-2256	true
4	10	67.1	-1611	true

D-Wave で解くには少し準備が必要です。まず Amazon Braket では2つデバイスが用意されているので、好きな方を選びます。次に、そのデバイスの各量子ビットに解きたい問題を対応させます（埋め込み）。
あとはシミュレーテッドアニーリングと同じです。Amazon Braket では、量子アニーリング 1-shot 毎にお金がかかるので注意しましょう。僕は10,000円使いました。

D-Waveによる求解

import boto3
from braket.aws import AwsDevice
from braket.ocean_plugin import BraketSampler, BraketDWaveSampler
from dwave.system.composites import EmbeddingComposite

aws_account_id = boto3.client("sts").get_caller_identity()["Account"]
bucket = '<your bucket name>'
prefix = "d-wave/output"
s3_location = (bucket, prefix)

# 2種類のデバイスから選択
# device = AwsDevice("arn:aws:braket:::device/qpu/d-wave/DW_2000Q_6")
device = AwsDevice("arn:aws:braket:::device/qpu/d-wave/Advantage_system1")

# トポロジーに埋め込む
sampler = BraketDWaveSampler(s3_location, device.arn)
sampler = EmbeddingComposite(sampler)

# 1000サンプル取得
%time sampleset = sampler.sample(bqm, num_reads=1000, auto_scale=True)
decoded_samples = model.decode_sampleset(sampleset)
best_sample = min(decoded_samples, key=lambda x: x.energy)

結果は以下のとおりです。最適解の前に実行可能解（制約を満たす解）を見つけられていないのでなかなか厳しいですね。
もちろんパラメーター調整でどうにかなる可能性はありますが、そもそも論として問題規模が大きいと問題を埋め込めない事がわかります。

D-Waveの結果

D-Waveをそのまま使っただけでは、 2x3 の小さな問題も解けないようです。また、4x5 以上では ValueError: no embedding found というエラーを吐いて解くことができません。これは、量子アニーラーのハードウェア的な制約で、解ける問題のサイズが限られているために起こるエラーです。現在の D-Wave では非常に小さな問題しか解けないようです。

行数	列数	求解時間(sec)	目的関数	制約充足	結果
2	3	10.6	34156	false	-
2	5	144.8	61046	false	-
4	5	-	-	-	`ValueError: no embedding found`

シミュレーテッド分岐マシンによる解法

さて、D-Wave をそのまま利用するだけでは、あまり大きな問題を解けないことがわかりました。もちろん、量子アニーラーや量子コンピューターの発展は目覚ましいものがあるので、今後に期待したいところですが、現状のままでは我らが公務員の待遇改善には寄与できません。

そこで、東芝が開発したシミュレーテッド分岐マシン(SBM; Simulated Bifurcation Machine) を試してみることにします。SBMは2019年に発表されたアルゴリズムで、当時結構話題になりました。現在その実装が試験的に無料^[4]で公開されていて、10000ビットまで等の制約はあるものの、AWSのアカウントさえあればだれでも使えます。HTTP インターフェースが用意されているため、cURL を使って問題を解くことができるのもの面白い点かと思います。

cURLをつかっても良いのですが、今回はせっかく PyQUBO を使ってモデルを記述したので、 PyQUBO で使える Sampler を書いてみました。

PyQUBO用に書いたSBMSampler

from io import BytesIO
import json
from urllib import request

import dimod
import scipy.sparse as sp
from scipy.io import mmwrite

class SBMSampler(dimod.Sampler, dimod.Initialized):
    parameters = None
    properties = None

    SCHEME = 'http'
    PATH = 'solver/ising'
    HTTP_HEADERS = {'Content-Type': 'application/octet-stream'}

    def __init__(self, host, port=8000):
        self.host = host
        self.port = port
        self.parameters = {
            'steps': [],
            'loops': [],
            'timeout': [],
            'maxwait': [],
            'target': [],
            'prefer': ['speed', 'auto'],
            'stats': ['none', 'summary', 'full'],
            'dt': [],
            'C': []
        }
        self.properties = {}

    def sample(self, bqm, steps=0, loops=1, timeout=None, maxwait=None, target=None, prefer='auto', stats=None, dt=1.0, C=0):
        biases, offset = bqm.to_qubo()
        keys = list(bqm.variables)
        key2idx = {key:idx for idx, key in enumerate(keys)}
        n_vars = len(keys)

        mat = sp.dok_matrix((n_vars, n_vars))
        for (key1, key2), bias in biases.items():
            row, col = key2idx[key1], key2idx[key2]
            mat[row, col] = bias
            if row != col:
                mat[col, row] = bias

        buff = BytesIO()
        mmwrite(buff, mat, symmetry='symmetric')
        buff.seek(0)
        req_body = buff.read()

        query = urlencode(
            self._delete_if_none(
                dict(steps=steps, loops=loops, timeout=timeout, maxwait=maxwait, target=target, prefer=prefer, stats=stats, dt=dt, C=C)
            )
        )
        url = f'{self.SCHEME}://{self.host}:{self.port}/{self.PATH}?{query}'
        req = request.Request(url, req_body, self.HTTP_HEADERS)

        with request.urlopen(req) as res_:
            res = json.loads(res_.read())

        samples = np.array([res.pop('result')])
        energy = res.pop('value') + bqm.offset

        return dimod.SampleSet.from_samples(
            (samples, keys),
            energy=energy,
            info=res,
            vartype=dimod.SPIN
        )

    def _delete_if_none(self, d):
        return {k:v for k, v in d.items() if v != None}

使い方は SimulatedAnnealingSampler とほぼ同じです。

# SBMの動いているEC2インスタンスを指定する
sampler = SBMSampler(SBM_HOST)
# num_reads ではなく loops。 SimulatedAnealingSampler と違い、最良の結果1つだけが返される
samples = sampler.sample(bqm, loops=1000)

SBMの結果

結果は以下のとおりです。確かに大きな問題を入力できましたが、だからと言って良い結果が得られるわけではないようです。もちろん、他の手法と同じく今回はパラメーターの調整をほとんど行っていませんので、これだけで結果の良し悪しが決まるわけではない点に注意が必要です。

Simulated Bifurcation Machine の結果

行数	列数	求解時間(sec)	目的関数	制約充足	結果
2	3	2.0	-2080	true
2	5	2.0	2389	false
4	5	2.7	150952	false
2	10	2.7	209275	false
4	10	-	-	-	HTTPError: HTTP Error 413: Payload Too Large^[5]

D-Wave Hybrid による解法

話は D-Wave に戻ります。先程、大きな問題は D-Wave では解けないと書きましたが、問題を小さな部分問題の組合せに分割できれば、ある程度大きな問題でも D-Wave で解くことができます。以前は qbsolv というライブラリがその役を担っていたのですが、現在は D-Wave Hybrid というライブラリの利用が推奨されているようです。D-Wave Hybrid は問題を分割するだけでなく、古典的なアルゴリズムと量子アルゴリズムを並行して実行することで、より良い解を高速に求めることができます。

D-Wave Hybrid の利用イメージ https://docs.ocean.dwavesys.com/projects/hybrid/en/latest/intro/overview.html

せっかくなので使ってみましょう。定式化の部分はこれまでのものと同じものを利用できます。

D-Wave Hybrid を用いた求解

def sample(sampler, bqm):
    iteration = hybrid.RacingBranches(
        hybrid.InterruptableTabuSampler(),
        hybrid.EnergyImpactDecomposer(size=2)
        | hybrid.QPUSubproblemAutoEmbeddingSampler(qpu_sampler=sampler)
        | hybrid.SplatComposer()
    ) | hybrid.ArgMin()
    workflow = hybrid.LoopUntilNoImprovement(iteration, convergence=3)

    init_state = hybrid.State.from_problem(bqm)
    final_state = workflow.run(init_state).result()

    return final_state.samples

sampler = BraketDWaveSampler(s3_location, device.arn)
sampler = EmbeddingComposite(sampler)
samples = sample(sampler, bqm)

大規模な問題を部分問題に分割して解こうとすると、どうしても処理が複雑になってしまいます。D-Wave Hybrid を使うと、各処理をパイプ | でつないでワークフローを構築することで、複雑な処理の見通しが立ちやすくなります。上のコードの例では、以下の流れのワークフローを構築しています。

hybrid.InterruptableTabuSampler() でサンプルを取得する（古典的なアルゴリズムのタブーサーチを使った求解）
1.と同時に問題を分割して D-Wave で求解
- hybrid.EnergyImpactDecomposer(size=2) で問題を分割して
- hybrid.QPUSubproblemAutoEmbeddingSampler(qpu_sampler=sampler) で部分問題を D-Wave で解いて
- hybrid.SplatComposer() で D-Wave の出力結果をまとめ上げる
hybrid.ArgMin() で1. と 2. で良い方を採用
hybrid.LoopUntilNoImprovement で解が更新されなくなりまで繰り返す

D-Wave Hybrid の結果

結果は以下のとおりです。素の D-Wave やシミュレーテッドアニーリングを利用した場合と比べて、随分良くなっています。2x30 でも最適解でこそありませんが、制約は満たしていますし、画像を見る限りそこそこ良い結果のように思います。

D-Wave Hybrid による求解結果

行数	列数	求解時間(sec)	目的関数	制約充足
2	3	20.6	-2080	true
2	5	15.3	-3070	true
4	5	19.6	-4126	true
2	10	19.3	-6340	true
4	10	41.6	-6340	true
2	30	302	-25937	true

本当にDWaveが必要なのか？

D-Wave Hybrid、なかなか良さそうですね。しかし、ここで一つ疑問が出てきます。D-Wave Hybrid では、古典アルゴリズムと量子アルゴリズムの両方を実行し、良いほうを採用しています。本当に量子アルゴリズムが必要なのでしょうか？前述の通り、シミュレーテッドアニーリングとD-Waveではシミュレーテッドアニーリングのほうがむしろ良い結果を出していました。

そこで、ワークフローを少し修正し、 D-Wave を利用していた部分を古典的なアルゴリズムのに置き換えた結果を見てみたいと思います。

タブーサーチとシミュレーテッドアニーリングの Hybrid 実装

def sample(bqm):
    iteration = hybrid.RacingBranches(
        hybrid.InterruptableTabuSampler(),
        hybrid.EnergyImpactDecomposer(size=2)
        | hybrid.SimulatedAnnealingSubproblemSampler()
        | hybrid.SplatComposer()
    ) | hybrid.ArgMin()
    workflow = hybrid.LoopUntilNoImprovement(iteration, convergence=3)

    init_state = hybrid.State.from_problem(bqm)
    final_state = workflow.run(init_state).result()
    
    return final_state.samples

結果は以下のとおりです。結果はだいたい D-Wave を使った時と同じで、速度はシミュレーテッドアニーリングのほうが早いようですが、問題サイズが大きくなるにつれて、速度の差はなくなてくるようです。より大規模な問題になると、D-Wave のほうが有利なのかもしれません。4x30になると ValueErrorが発生しているので、原因を追求する必要があります。

D-Wave Hybrid （シミュレーテッドアニーリング）による求解結果

行数	列数	求解時間(sec)	目的関数	制約充足	結果
2	3	0.44	-2080	true
2	5	0.56	-3070	true
4	5	2.71	-4126	true
2	10	3.01	-6340	true
4	10	25.5	-7396	true
2	30	125.4	-25293	true
4	30	-	-	-	`ValueError: mismatch between variables in 'initial_states' and 'bqm'`

問題の構造を使った定式化

ここまでは、素直な定式化 (1)-(2) をもとに最適解を求めていましたが、問題の構造を利用することで解きやすくなることがあります。
今回の問題についても、もう少し効率的な定式化方法がありそうです。いろいろな方法があると思いますが、ここでは、行数=1のときに巡回セールスマン問題に帰着できたことに注目し、それを拡張した形で定式化してみます。これまでの定式化と違い、各紙片同士が直接つながっているかどうかだけを考えます。

インデックスと定数

定数は今まで出てきたものと変わりません。定式化をシンプルにするため上下左右端を表すダミーの紙片を考え、追加しています。

$R$ : 縦方向の分割数
$C$ : 横方向の分割数
$N = R \times C$ : 紙片の総数
$S^{x}_{ab}$ : 紙片 $a$ の右に紙片 $b$ が置かれている場合の類似度
$S^{y}_{ab}$ : 紙片 $a$ の下に紙片 $b$ が置かれている場合の類似度
$a,b\in[0,1,\dots,N]$ : 紙片を表すインデックス（ $N$ は上下左右端を表すダミーの紙片を表す）

変数

巡回セールスマン問題にも様々な定式化方法がありますが、ここではMTZ形式と呼ばれる方法を拡張することを考え、ポテンシャル変数を導入します。

$x_{ab} \in \{0, 1\}$ : 紙片 $a$ の右に紙片 $b$ が置かれている場合は $1$ となるフラグ
$y_{ab} \in \{0, 1\}$ : 紙片 $a$ の下に紙片 $b$ が置かれている場合は $1$ となるフラグ
$u^{x}_{a}$ : 横方向のポテンシャル（横方向の位置を表している）
$u^{y}_{a}$ : 縦方向のポテンシャル（縦方向の位置を表している）

目的関数

目的関数はこれまでと似た形をしていますが、紙片の絶対的な位置に依存していない点が異なります。

\max \sum_{ab}S^{x}_{ab}x_{ab} + \sum_{ab}S^{y}_{ab}x_{ab}

制約条件

制約条件は以下のとおりです。はじめの4つは、各紙片の左右上下に1つづつ別の紙片がおかれていること、のこりの4つは、１つ右（下）にずれると $u_x$ （ $u_y$ ）が1つ増えること、同じ列（行）であれば $u_x$ （ $u_y$ ）は等しいことを表しています。(12/1 コードに合わせて数式を修正)

\begin{aligned} \sum_{b} x_{ab} &= 1, \quad (a \neq N) \\ \sum_{b} x_{ba} &= 1, \quad (a \neq N) \\ \sum_{b} y_{ab} &= 1, \quad (a \neq N) \\ \sum_{b} y_{ba} &= 1, \quad (a \neq N) \\ \sum_{b} x_{Nb} &= R, \\ \sum_{b} x_{bN} &= R, \\ \sum_{b} y_{Nb} &= C, \\ \sum_{b} y_{bN} &= C, \\ u^{x}_{a} + 1 - N(1 - x_{ab}) &\leq u^{x}_{b}, \\ u^{y}_{a} + 1 - N(1 - y_{ab}) &\leq u^{y}_{b}, \\ u^{x}_{a} - N(1 - y_{ab} - y_{ba}) &\leq u^{x}_{b}, \\ u^{y}_{a} - N(1 - x_{ab} - x_{ba}) &\leq u^{y}_{b}, \\ x_{ab} + x_{ba} &\leq 1, \\ y_{ab} + y_{ba} &\leq 1, \\ x_{ab} + y_{ab} &\leq 1, \\ x_{ab} + y_{ba} &\leq 1. \end{aligned}

PuLPによる実装

PuLPでの実装は以下のようになります。これまでのものと比べて若干複雑になっています。

PuLPによる実装

def build_model(data):
    model = pulp.LpProblem(sense=pulp.LpMaximize)
    x = np.array(pulp.LpVariable.matrix(
        'x',
        (range(data.size + 1), range(data.size + 1)),
        cat=pulp.LpBinary
    ))
    y = np.array(pulp.LpVariable.matrix(
        'y',
        (range(data.size + 1), range(data.size + 1)),
        cat=pulp.LpBinary
    ))
    ux = np.array(pulp.LpVariable.matrix(
        'ux',
        (range(data.size),),
        lowBound=1.0,
        upBound=data.cols,
        cat=pulp.LpContinuous
    ))
    uy = np.array(pulp.LpVariable.matrix(
        'uy',
        (range(data.size),),
        lowBound=1.0,
        upBound=data.rows,
        cat=pulp.LpContinuous
    ))
    for a in range(data.size + 1):
        x[a, a].upBound = 0
        x[a, a].lowBound = 0
        y[a, a].upBound = 0
        y[a, a].lowBound = 0

    model.setObjective(pulp.lpSum(
            sim_x(data.images[a], data.images[b])*x[a, b]
            for a in range(data.size)
            for b in range(data.size)
            if a != b
        ) + pulp.lpSum(
            sim_y(data.images[a], data.images[b])*y[a, b]
            for a in range(data.size)
            for b in range(data.size)
            if a != b
        )
    )

    for a in range(data.size):
        model.addConstraint(pulp.lpSum(x[a]) == 1)
        model.addConstraint(pulp.lpSum(x[:, a]) == 1)
        model.addConstraint(pulp.lpSum(y[a]) == 1)
        model.addConstraint(pulp.lpSum(y[:, a]) == 1)
    model.addConstraint(pulp.lpSum(x[-1]) == data.rows)
    model.addConstraint(pulp.lpSum(x[:-1, -1]) == data.rows)
    model.addConstraint(pulp.lpSum(y[-1]) == data.cols)
    model.addConstraint(pulp.lpSum(y[:, -1]) == data.cols)


    for a in range(data.size):
        for b in range(data.size):
            if a != b:
                model.addConstraint(
                    ux[a] + 1 - data.cols*(1 - x[a, b]) <= ux[b]
                )
                model.addConstraint(
                    ux[a] - data.cols*(1 - y[a, b] - y[b, a]) <= ux[b]
                )
                model.addConstraint(
                    uy[a] + 1 - data.rows*(1 - y[a, b]) <= uy[b]
                )
                model.addConstraint(
                    uy[a] - data.rows*(1 - x[a, b] - x[b, a]) <= uy[b]
                )
                model.addConstraint(
                    x[a, b] + x[b, a] <= 1
                )
                model.addConstraint(
                    y[a, b] + y[b, a] <= 1
                )
                model.addConstraint(
                    x[a, b] + y[a, b] <= 1
                )
                model.addConstraint(
                    x[a, b] + y[b, a] <= 1
                )
    return model, x, y, ux, uy

問題の構造を利用した定式化の結果

結果は以下の様になりました。2x30についても最適解を見つけることができました。4x30についてもタイムアウトで最適解は見つけられなかったものの、なかなか良い解が得られているようです。問題の構造を利用するのは大事だというのがわかる一方で、複数の書類を同時にシュレッダーにかけたような現実的な状況で使えるかといわれると、少なくともこの方法では難しそうだと言えるかと思います。

MIPによる定式化の結果

行数	列数	求解時間(sec)	最適値
2	3	0.035	-2080
2	5	0.11	-3070
4	5	0.57	-4126
2	10	0.79	-6340
4	10	6.25	-7396
2	30	429	-27213
4	30	1860	-28358(Not Solved)

商用ソルバで解いてみる

MIPソルバーには商用のものも非商用のものもありますが、商用ソルバーのほうが圧倒的に性能が高く、大規模な問題を解くには商用ソルバーが必須であることが知られています。そこで、本記事の締めくくりとして、今手元にあって自由に使える Gurobi 7.0 の結果を示そうと思います。数年前のモデルで、かつインストールされている環境（PC）も古いので、Gurobi に不利な比較となっていますが、非商用のソルバーと比較してどれくらい差があるかはわかるかと思います。

結果は以下のとおりです。想像通り、Cbcでは厳しかった 4x30 や 2x50 も最適解を見つけられています。一方で 4x50 については30分以内に結果を得ることはできませんでした。問題サイズが大きくなると急激に時間がかかるのは、Cbcと同じようです。

Gurobiによる求解

行数	列数	求解時間(sec)	最適値
2	3	0.015	-2080
2	5	0.035	-3070
4	5	0.11	-4126
2	10	0.11	-6340
4	10	0.50	-7396
2	30	2.43	-27213
4	30	34.7	-28269
2	50	36.8	-33692

まとめ

シュレッダーで裁断された書類の復元問題に色々な方法で取り組んでみました。まだまだ量子アニーリングは難しいなと思った一方、数年前に量子アニーリングに取り組んだ時と比べると、格段にできることが増えているな、という印象をもちました。残念ながら国家公務員の労働環境改善にはつながらなそうですが、個人的には楽しかったので満足です。

ちなみにシュレッダーの復元問題は色々研究がされているようです。特徴量をどうするか、とか深層学習を使ってみると、とかいろいろやられていて面白そうです。

脚注

$x$ の値によらず常に $y=0$ としても制約を満たしてしまいますが、目的関数によって $y$ が可能な限り大きな値を取ろうとするので、この制約を課すだけで、最適解では $y_{a,b,i,j,k,l} = \min(x_{a,i,j},x_{b,k,l})$ となります。 ↩︎
明らかに最適解ではなさそうですが、PuLPの status は Optimal となっています。計算誤差の可能性もあるかもしれません ↩︎
今回はあくまでいろいろな手法を試してみることを目的としているので、デフォルトパラメーターをそのまま利用しています ↩︎
ソフトウェア利用料は無料ですが、EC2のインスタンス利用料はかかります。利用できるインスタンスタイプは p3.2xlarge のみなので、$3.06/hr かかります（2020年12月現在）。 ↩︎
payloadの最大値は設定ファイルを書き換えることで変更可能です。 ↩︎

DARPA Shredder Challenge 2011

裁断された書類の復元と巡回セールスマン問題

問題設定

素朴な定式化

（混合）整数線形計画法による解法

PuLP による実装と結果

量子アニーラ（量子コンピューター）による解法

PyQUBO を用いた実装と結果

シミュレーテッド分岐マシンによる解法

SBMの結果

D-Wave Hybrid による解法

D-Wave Hybrid の結果

本当にDWaveが必要なのか？

問題の構造を使った定式化

インデックスと定数

変数

目的関数

制約条件

PuLPによる実装

問題の構造を利用した定式化の結果

商用ソルバで解いてみる

まとめ

Discussion