D-WaveでMaxCutをN番煎じ

概略

色々めんどくさくてやる気が起きないが、D-Waveの新しいQPUが結構よさげだったので情報共有。
簡単にまとめると、

• D-Waveに新しいQPU(Advantage2 Prototype2.2がリリースされていた)
• 何の気なしに試してみたら、結構性能が出た

パラメータ調整はせずに昔のコードを実行した結果Gurobiに一部勝結果をだした。（Size45ではGurobiよりも早い処理時間をたっせいし、40%の確率で最適解をだしている）

線形計画では全くかなわないが、一部の領域では使える目がでてきたのではないか。

D-Wave Advantage 2.2

悲しいことに話題になっていないが、D-Waveの新しいQPUがPrototypeではあるが2024年2月にリリースされている。
これ、無償のIDでもお試しできる。

3年前

3年前に実施した記事は以下。MaxCutは2次のため、当時はCBCなどでは解けず、全探索 vs D-Waveを実施した。
https://qiita.com/___monta___/items/bef8ebe88cc6a5d51944

＠2024年2月

当時は全探索しかなかったが、Gurobipyも小さなサイズであればお試しで使用できるようになった。
コードも動作しなくなっていたため、微修正して動くようにして性能比較をした。

実行前のお約束

ぐぐるなり、この辺ご参考に設定するとよいでしょう。

大体わかる人向けに3行でまとめると

• Leapに登録
• pip install dwave-ocean-sdk
• dwave config create 実行でTokenなど登録

修正したコード

全量以下でサイズを変えたベンチマークが走るようにコードをいじっているので、興味がある方は直接コード修正お願いします。

Maxcut

import itertools
from pyqubo import Array
from dwave.system import DWaveSampler,EmbeddingComposite
import networkx as nx
from matplotlib import pyplot as plt
import time
import argparse
import random
import numpy as np
import gurobipy as gp

class MyTimer():
    '''
    時刻計測用 のクラス
    '''
    def __init__(self):
        self.set_timer()

    def set_timer(self):
        self.MY_TIMER = time.time()

    def get_timer(self):
        current_time = time.time()
        ret_time = current_time - self.MY_TIMER
        self.MY_TIMER = current_time
        return(ret_time)

def get_args():
    parser = argparse.ArgumentParser()
    parser.add_argument('-n', '--node', default=7, type=int, help='This is number of nodes')
    parser.add_argument('-p', '--percentage', default=0.6, type=float, help='This is Probability for edge creation')
    parser.add_argument('-s', '--seed', default=0,type=int, help='This is seeds of random')

    return parser.parse_args()

class Problem():
    def __init__(self,n,p,seed):
        self.p = p
        self.n = n
        self.seed = seed
        self.G = nx.fast_gnp_random_graph(n,p,seed=seed)

        for v in self.G.nodes:
            if len(self.G.edges(v)) == 0:
                self.G.add_edge(v,random.randint(0,n))

        #重みづけ
        for (u,v) in self.G.edges:
            self.G.edges[u,v]['w'] = random.randrange(1,30,1)


        self.pos = nx.layout.spring_layout(self.G)

    def solve_dwave(self):
        timer = MyTimer()
        timer.set_timer()

        X = Array.create('X',shape=(self.n),vartype='SPIN')

        H0 = 0
        for (u,v) in self.G.edges:
            H0 -= (1-X[u]*X[v])*self.G.edges[u,v]['w']/2.0

        model = H0.compile()
        bqm = model.to_bqm()
        bqm.normalize()

        response = EmbeddingComposite(DWaveSampler(solver='Advantage2_prototype2.2')).sample(bqm,num_reads=3000)

        top_response = response.first

        self.cutted = []
        for (u,v) in self.G.edges:
            u_spin = int(top_response.sample[f"X[{u}]"])
            v_spin = int(top_response.sample[f"X[{v}]"])
            if u_spin != v_spin:
                self.cutted.append([u,v])
                
        self.cal_energy()   
        
        print(f"RES:D-Wave size = {self.n} time = {timer.get_timer():.2f} , Energy = {self.energy} ")

    def solve_gurobi(self):
        timer = MyTimer()
        timer.set_timer()

        m = gp.Model("MaxCut")

        X = {}
        for i in range(self.n):
            X[i] = m.addVar(vtype=gp.GRB.BINARY)

        m.setObjective(sum((2*X[u]*X[v]-X[u]-X[v])*int(self.G.edges[u,v]['w']) for u,v in self.G.edges))

        m.optimize()

        self.cutted = []
        for (u,v) in self.G.edges:
            u_spin = int(X[u].X)
            v_spin = int(X[v].X)
            if u_spin != v_spin:
                self.cutted.append([u,v])
                
        self.cal_energy()   

        print(f"RES:Gurobi size = {self.n} time = {timer.get_timer():.2f} , Energy = {self.energy} ")
        
    def cal_energy(self):
        
        self.energy = 0
        for  (u,v) in self.cutted:
            self.energy -= self.G.edges[u,v]['w']
            
    def show(self):
        plt.figure(figsize=(12,12),dpi=100)
        nx.draw_networkx(self.G,pos=self.pos)
        if len(self.cutted)>0:
            nx.draw_networkx_edges(self.G,self.pos,self.cutted,width=2,edge_color='r')
        edge_labels={}
        for (u,v) in self.G.edges:
            edge_labels[(u,v)] = self.G.edges[u,v]['w'] 
        nx.draw_networkx_edge_labels(self.G,self.pos,edge_labels=edge_labels)
        nx.draw_networkx_nodes(self.G,self.pos,nodelist=[idx for idx,v in enumerate(self.v_group) if v == 0],node_color='yellow')

        plt.show()

    def solve_bf(self):
        timer = MyTimer()
        timer.set_timer()
        min_value = 0
        for node_status in itertools.product((0,1),repeat=self.n):
            tmp_min = 0
            for (u,v) in self.G.edges:
                if node_status[u] != node_status[v]:
                    tmp_min -= self.G.edges[u,v]['w'] 

            if tmp_min < min_value:
                min_value = tmp_min

        print(f"RES:Brute-Force size = {self.n} ,time = {timer.get_timer():.2f}, Energy = {min_value}")


if __name__ == '__main__':
    args = get_args()
    #p = Problem(n=args.node,p=args.percentage,seed=args.seed)

    #p.solve_dwave()
    #p.show()
    for j in range(5):
        for i in range(20,65,5):
            p = Problem(n=i,p=args.percentage,seed=args.seed)
            p.solve_dwave()
            p.solve_gurobi()
            #solve_bf()

結果

一部数字が重なっていて申し訳ないですが、

• 問題サイズが非常に小さい（15未満）はどの手法でも大丈夫
• 問題数がまぁまぁ小さい（15以上20程度）は全探索しんどい。GurobiもD-Waveも余裕

D-Waveは前処理、ネットワーク越しなこともあり問題サイズにかかわらず、3-4秒程度

問題サイズが20を超えたもの。

• 45をこえるとGurobiが急に厳しくなる。
• D-Waveは前処理時間などあるが処理時間ののびは緩やかMax10秒程度
• 最適解だけ見るとサイズ45を超えるとD-Waveでは算出できていないが、そこそこの近似解は出ている。