Go to Qiita Advent Calendar 2024 Top
LoginSignup
20
10

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@ynntech(Yushi Nakaya)

頂点被覆問題のコスト関数をPyQUBOで実装しOpenJijで評価する方法

Last updated at Posted at 2019-04-16

本稿では，頂点被覆問題のコスト関数から，QUBO(Quadratic unconstrained binary optimization)を，従来の自分で重み付けを行う方法とPyQUBOによって自動生成する方法の2種類それぞれによって作成し，PyQUBO付属のSA (Simulated Annealing)とOpenJijによるSQA（Simulated Quantum Annealing)によって評価する方法を紹介する．

#定式化
最小頂点被覆問題では，部分集合$V'$の要素数が最小となるものを求めるを考える．
彩色する頂点数($V'$に属する頂点数)を最小化するため，この最小化項がコスト関数の主な項となる．

H_{A} = \sum_{i \in V} q_v

ちなみに，各辺のいずれかの端点はどこかが彩色されるという条件を以下のペナルティ項として与える．

H_{B}=\sum_{(i,j)\in E} (1-q_i)(1-q_j)

以上の最小化項，ペナルティ項より，グラフ$G=(V,E)$の頂点被覆問題のコスト関数は次のように導かれる．

E = H_{A}＋H_{B} = \sum_{i \in V} q_v + \sum_{(i,j)\in E} (1-q_i)(1-q_j)

#自分でQUBOを打ち込む従来の方法

頂点被覆問題をPythonで実装すると以下のようになる．
評価はOpenJijによるSimulated Quantum Annealingを用いている．

今回は頂点被覆問題の初期条件を以下のように与える．


#グラフを定義．ここでは0~7と名前をつけた8個の頂点 を持つグラフを定義．
#(0,1)の場合，頂点0と頂点1がつながっていることを示す．
graph = [(0,1), (0,2), (1,3), (5,6), (2,3), (2,5), (3,4), (5,7)]

##実装

#openJijを使用．
import openjij as oj
graph = [(0,1), (0,2), (1,3), (5,6), (2,3), (2,5), (3,4), (5,7)]

#ここで，QUBOを作成．
indices = set([]) 
for index in graph:
    indices = indices | set(index)   
indices = list(indices)    
    
#ここから，この問題のQUBOを手動で作成．
Q = {}
#まず，最小化項を作成する．すべての頂点に対して1を代入．
for i in indices:
    Q[(i,i)] = 1
    
#ペナルティ項の作成．
#ここで，各辺に対して，各頂点が何回使われるかに注目し，
#より多くの辺の端点になっている点のコストを小さくする．
for i,j in graph:
    Q[(i,j)] = 1
    Q[(i, i)] += -1
    Q[(j, j)] += -1

#結果を出力．SASamplerではSimulated Annealing，SASamplerではQuantum Annealing方式で最適化項を探す．
#iterationで計算回数，step_numで一回のアニーリングあたりの計算時間を表す．
results = oj.SQASampler(iteration=30, step_num=300).sample_qubo(Q)

計算結果を出力すると以下のようになる．

results.states

[[1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [0, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 1, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [0, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0],
 [1, 0, 0, 1, 0, 1, 0, 0]]

また，それぞれの計算に要したエネルギーを以下のように出力できる．

results.energies

[-1.75,
 -1.75,
 -1.75,
 -1.75,
 -1.25,
 -1.75,
 -1.25,
 -1.75,
 -1.75,
 -1.75,
 -1.75,
 -1.75,
 -1.75,
 -1.75,
 -1.75,
 -1.25,
 -1.25,
 -1.75,
 -1.75,
 -1.75]

次に，PyQUBOによってQUBOを作成する方法を紹介する．ここでは，OpenJijによる評価とPyQUBOに搭載されているSA(Simulated Annealing)によって評価する方法の2種類を紹介する．
#PyQUBOによってQUBOを作成する
##使用した環境
※使用した環境は以下の通り．

pyqubo==0.1.3

numpy==1.15.4
dimod==0.7.9
six==1.11.0

dwave-cloud-client==0.4.16
dwave-neal==0.4.2
dwave-networkx==0.6.9
dwave-ocean-sdk==1.0.3
dwave-qbsolv==0.2.7
dwave-system==0.5.4
dwavebinarycsp==0.0.10

順番を間違えるとバージョンが異なるものに書き換えられる可能性もあるので，最後に
pip freeze
で必ず確認する．

##サンプル

配列を作成する際の例

x = Array.create('x', shape=4, vartype='BINARY')

shape=4で4つの1次元配列をつくる．

ちなみに，２次元配列を作るときは以下のとおり．

from pyqubo import Array, Binary
 Array.create('x', shape=(2, 2), vartype='BINARY')
 Array([[Binary(x[0][0]), Binary(x[0][1])],
       [Binary(x[1][0]), Binary(x[1][1])]])

shape=(2, 2)で2×2の上記のような2次元配列を作っている．

##PyQUBOによるSimulated Annealing
ここから，実際にPyQUBOを用いて最適化問題をSAで解く方法を説明していく．
まずは実装例を以下に示す．

###頂点被覆問題 実装例

import openjij as oj
from pyqubo import Array, Constraint, Placeholder, solve_qubo

graph = [(0,1), (0,2), (1,3), (5,6), (2,3), (2,5), (3,4), (5,7)]

# BINARY変数
x = Array.create('x', shape=8, vartype='BINARY')
    
param_cover = Placeholder("cover")
H_cover = Constraint(sum((1-x[u])*(1-x[v]) 
    for(u,v) in graph), "cover")
    
#最小化項
H_A = sum(x) 
#ペナルティ項
H_B = H_cover 

#コスト関数
Q = H_A+H_B 
    
# モデルをコンパイル
model = Q.compile()
feed_dict = {'cover': param_cover}
qubo, offset = model.to_qubo(feed_dict=feed_dict)
    
#PyQUBOによるSA
raw_solution = solve_qubo(qubo)
decoded_solution, broken, energy = model.decode_solution(
    raw_solution, 
    vartype="BINARY", 
    feed_dict=feed_dict
)

#PyQUBOによるSAの結果を出力
print(decoded_solution)

ここで，上記コードの詳細について説明していく．
まず，PyQUBOをimport．

#PyQUBO
from pyqubo import Array, Constraint, Placeholder, solve_qubo

コスト関数に用いるxを以下のように定義．

# BINARY変数 PyQUBOでは，以下のようにArray.createという関数で配列を作ることができる．
x = Array.create('x', shape=(8), vartype='BINARY')

次に，最小化項とペナルティ項を以下のように打ち込み，コスト関数を定義する．
通常のQUBOを作成する方法とはことなり，式をそのまま直感的に打つこむことができる．

# プレースホルダー
param_cover = Placeholder("cover")
# QUBO形式で定式化
H_cover = Constraint(sum((1-x[u])*(1-x[v]) 
    for(u,v) in graph), "cover")
    
#最小化項
H_A = sum(x)
#ペナルティ項
H_B = H_cover*param_cover 

#コスト関数
Q = H_A+H_B

作成したコスト関数は，以下のようにQ.compile().to_qubo() で簡単にQUBOに変換することができる．

# モデルをコンパイル
model = Q.compile()
# プレースホルダーと合わせてQUBOを作成
feed_dict = {'cover': param_cover}
qubo, offset = model.to_qubo(feed_dict=feed_dict)

PyQUBOによってSimulated Annealingするには，
solve_qubo(qubo)を用いる．

#PyQUBOによるSA
raw_solution = solve_qubo(qubo)

# 得られた結果をデコードする
decoded_solution, broken, energy = model.decode_solution(
    raw_solution, 
    vartype="BINARY", 
    feed_dict=feed_dict
)

そして最後に，PyQUBOによるSAの結果を出力する．

print(decoded_solution)

##PyQUBOによってQUBOに変換したモデルのOpenJijによる評価
次に，上のコードに追記する形で，OpenJijにQUBOを受け渡し，SQAで解く方法を以下で紹介する．

自分でQUBOを打ち込む従来の方法でQを代入していたsampler.sample_qubo()にquboを受け渡す．
あとは通常のQUBOを作成してOpenJijで評価する場合と同じ．

sampler = oj.SQASampler(iteration=10, step_num=300)
response = sampler.sample_qubo(qubo)
#結果を出力．
print(response.states)

まとめると，以下のようになる．

import openjij as oj
from pyqubo import Array, Constraint, Placeholder, solve_qubo

graph = [(0,1), (0,2), (1,3), (5,6), (2,3), (2,5), (3,4), (5,7)]

#ここで、自動でQUBOを作成。
# BINARY変数
x = Array.create('x', shape=8, vartype='BINARY')
    
# プレースホルダー    
param_cover = Placeholder("cover")    
# QUBO形式で定式化
H_cover = Constraint(sum((1-x[u])*(1-x[v])        
for(u,v) in graph),"cover")
        
#最小化項   
H_A = sum(x)     
#ペナルティ項   
H_B = H_cover 

#コスト関数
Q = H_A+H_B 
    
# モデルをコンパイル
model = Q.compile()
feed_dict = {'cover': param_cover}
qubo, offset = model.to_qubo(feed_dict=feed_dict)
    
#PyQUBOによるSA
raw_solution = solve_qubo(qubo)

# 得られた結果をデコードする
decoded_solution, broken, energy = model.decode_solution(
    raw_solution, 
    vartype="BINARY", 
    feed_dict=feed_dict
)

#あとは通常のQUBOを作成してOpenJijで評価する場合と同じ．
sampler = oj.SQASampler(iteration=10, step_num=300)
response = sampler.sample_qubo(qubo)

#結果を出力．
print(response.states)

参考：PyQUBO公式ドキュメント
https://pyqubo.readthedocs.io/en/latest/reference/array.html?highlight=arry%20create

20
10
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
20
10

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?