A coherent Ising machine for MAX-CUT problems : Performance evaluation against semidefinite programming relaxation and simulated annealing

0.87 MB·English

by  Yoshitaka Haribara

A coherent Ising machine for maximum cut problems : Performance evaluation
against semidefinite programming relaxation and simulated annealing
Yoshitaka Haribara,1,2,∗ Shoko Utsunomiya,2,† Ken-ichi Kawarabayashi,2,‡ and Yoshihisa Yamamoto3,4,§
1Department of Mathematical Informatics, The University of Tokyo,
Hongo 7-3-1, Bunkyo-ku, Tokyo 113-8656, Japan
2National Institute of Informatics, Hitotsubashi 2-1-2, Chiyoda-ku, Tokyo 101-8403, Japan
3E. L. Ginzton Laboratory, Stanford University, Stanford, CA94305, USA
4ImPACT program, The Japan Science and Technology Agency,
Gobancho 7, Chiyoda-ku, Tokyo 102-0076, Japan
(Dated: March 1, 2016)
6
Combinatorialoptimizationproblemsarecomputationallyhardingeneral,buttheyareubiquitous
1
in our modern life. A coherent Ising machine (CIM) based on a multiple-pulse degenerate optical
0
parametric oscillator (DOPO) is an alternative approach to solve these problems by a specialized
2
physical computing system. To evaluate its potential performance, computational experiments are

b performed on maximum cut (MAX-CUT)problems against traditional algorithms such as semidef-
e inite programming relaxation of Goemans-Williamson and simulated annealing by Kirkpatrick, et
F al. Thenumericalresultsempiricallysuggestthatthealmostconstantcomputationtimeisrequired
to obtain the reasonably accurate solutions of MAX-CUT problems on a CIM with the number of
9
vertices up to2×104 and thenumberof edges up to108.
2

]
h I. INTRODUCTION where is an Ising Hamiltonian defined in Eq. (1) with
p J = Hw . It indicates that the MAX-CUT problem is
ij ij
- −
t Combinatorial optimization appears in many impor- equivalent to the Ising problem except for the constant
n tantfieldssuchascomputerscience[1,2],drugdiscovery factor.
a
u and life-science [3], and information processing technol- TheMAX-CUT problembelongstothe NP-hardclass
q ogy [4]. One of the example of such problems is an Ising ingeneral,eventhoughthere aregraphtopologieswhich
[ problem to minimize the Ising Hamiltonian, which is a can be solved in polynomial time [5, 7–11]. Many at-

function of a spin configuration σ =(σ ) defined as tempts have been made to approximatelysolve NP-hard
5 i
MAX-CUT problems,buttheprobabilisticallycheckable
v
0 H(σ)=− Jijσiσj (1≤i,j ≤N), (1) proof (PCP) theorem states that no polynomial time al-
3 Xi<j gorithms can approximate MAX-CUT problems better
0 than 0.94118 [12, 13]. Currently, the approximation ra-
where each spin takes binary values σ 1 , a real
7 i ∈ {± } tioof0.87856achievedbytheGoemans-Williamsonalgo-
number symmetric matrix J denotes a coupling con-
0 ij rithm (GW) based on semidefinite programming (SDP)
stant, and N is the total number of spins. Despite
.
1 isthebestvalueforperformanceguarantee[14]. Thisal-
its simple statement, it belongs to the non-deterministic
0 gorithm is a well-established benchmark to evaluate any
polynomial-time(NP)-hardclasstofindthegroundstate
5 new algorithms or computing methods.
of the Ising model on the three-dimensional lattice [5].
1
Besides, there exist several heuristic algorithms to
: Similarly,amaximumcut(MAX-CUT)probleminthe
v tackle these NP-hard MAX-CUT problems. The sim-
graph theory is to find the size of the largest cut in a
i ulated annealing (SA) was designed by mimicking the
X given undirected graph. Here, a cut is a partition of the
thermalannealingprocedureinmetallurgy[15]. Aquan-
verticesV into two disjointsubsets S ,S and the size
r 1 2
a ofthecutisthetotalweightofedges{w w}ithonevertex tum annealing technique was also formulated and was
ij
shown to have competitive performance against SA [16–
i in S and the other j in S . The size of the cut can
1 2
20]. Independently, novel algorithms which are superior
becountedbyassigningthebinaryspinvaluestoexpress
either in its speed or its accuracy are proposed [21–23].
which subset the vertex i belongs to σ 1 [6]:
i
∈{± } Werecentlyproposedanovelcomputingsystemtoim-
(1 σ σ )
i j plementthe NP-hardIsing problemsusingthe criticality
C(σ)= w = w −
ij ij
2 of laser [24–27] and degenerate optical parametric oscil-
i∈SX1,j∈S2 Xi<j
lator(DOPO)phasetransition[28,29]. Thearchitecture
1 1
= w (σ), (2) ofthismachineismotivatedbytheprincipleoflaserand
ij
2 − 2H
Xi<j DOPOin whichthe mode with the minimum loss rateis
most likely to be excited first. The energy of the Ising
Hamiltoniancanbemappedontothetotallossrateofthe
∗ [e-mail] laser or DOPO network. The selected oscillation mode
† [e-mail] inthelaserorDOPOnetworkcorrespondstotheground
‡ [e-mail] state of a given Ising Hamiltonian, while the gain acces-
§ [e-mail] sible to all other possible modes is depleted due to the
2
cross-gain saturation. This means that a mode with the In Sections IIID and IIIE, we assume a CIM with a
lowest loss rate reaches a threshold condition first and fiber length of 2 km (or cavity round trip time of 10 µs)
clumps the gain at its loss rate, so that all the other andpulsespacingof10cm(orpulserepetitionfrequency
modes with higher loss rates stay at sub-threshold con- of2GHz),thus2 104independentDOPOpulsescanbe
×
ditions. Moreover, the DOPO is in the linear superposi- prepared for computation. The system clock frequency
tion of 0-phase state and π-phase state at its oscillation for the CIM should be defined by the cavity circulation
threshold [30]. The coupled DOPOs form quantum en- frequency (inverse of cavity round trip time). One clock
tanglement in spite of their inherent dissipative natures cycle (round trip) includes every elements of computa-
[31], so that some form of quantum parallel searchcould tion,suchasparametricamplification,out-couplingport,
be embedded in a DOPO network. and coherent feedback. Thus the clock frequency of the
In this article, the validity of the CIM for MAX-CUT CIM assumed for the present benchmark study is 100
problemsis testedagainstthe representativeapproxima- kHz since the roundtriptime of2km fiber ringis 10µs.
tion algorithms. The DOPO signal pulse amplitudes in We fixed this system clock frequency, just like any dig-
CIM,whichareinterpretedasthesolution,aredescribed ital computer has a fixed clock frequency and chose the
bythec-numberstochasticdifferentialequations(CSDE) appropriate pulse interval to pack the desired number of
as presented in Section II. Then we conduct numerical pulses in the 2km fiber.
simulations for MAX-CUT problems in Section III with
the number of vertices up to N =20000. It is, of course,
difficulttocomparetheperformanceoftheproposedsys- B. c-number stochastic differential equations for
tem as a MAX-CUT solver with the representative ap- multiple-pulse DOPO with mutual coupling
proximationalgorithmswhichcanberunoncurrentdig-
ital computers mainly because the unit of “clock” can- The in-phase and quadrature-phase amplitudes of a
not be uniquely defined. Thus we defined the feasible single isolatedDOPOpulse obey the following c-number
system clock which dominates the computing process in stochastic differential equations (CSDE) [33]:
CIM as mentioned later. Moreover, here we evaluated
the computational ability under either time or accuracy 1 1
was fixed, while a preliminary benchmark study done in dc=(−1+p−c2−s2)cdt+ A rc2+s2+2dW1 (3)
s
the previous paper is focused on the performance after
1 1
physical convergence [29]. ds=( 1 p c2 s2)sdt+ c2+s2+ dW .(4)
− − − − A r 2 2
s
The above CSDE are derived by expanding the DOPO
II. MULTIPLE-PULSE DOPO WITH MUTUAL
field density operator with the truncated Wigner dis-
COUPLING
tribution functions. An alternative approach is to use
two coherent states in the generalized (off-diagonal)
A. A proposed machine P(α ,β )-representationfor the fielddensity matrix [34].
s s
The two approaches by the truncated Wigner function
A standard CIM based on multiple-pulse DOPO with and the generalised P-representation are equivalent for
all-optical mutual coupling circuits is shown in Fig. 1. highly dissipative systems such as ours. The pump
The system starts with a pulsed master laser at a wave- field is adiabatically eliminated in (3) and (4) by as-
lengthof1.56µm. Asecondharmonicgeneration(SHG) suming that the pump photon decay rate γ is much
p
crystal produces the pulse trains at a wavelength of larger than the signal photon decay rate γ . The term
s
0.78µmwhichinturngeneratemultipleDOPOpulsesat A = (γ γ /2κ2)1/2 is the DOPO field amplitude at a
s s p
a wavelength of 1.56 µm inside a fiber ring resonator. If normalized pump rate p=F /F =2, and κ is the sec-
p th
theroundtriptimeofthefiberringresonatorisproperly ond order nonlinear coefficient associated with the de-
adjusted to N times the pump pulse interval, we can si- generate optical parametric amplification. The variable
multaneously generate N independent DOPO pulses in- t = γ τ/2 is a normalized time, while τ is a real time
s
side the resonator. Each of these pulses is either in 0- in seconds. The term F is the pump field amplitude
p
phase state or π-phase state atwell abovethe oscillation and F = γ γ /4κ is the threshold pump field ampli-
th s p
threshold and represents an Ising spin of up or down. tude. Finally, dW and dW are two independent Gaus-
1 2
InordertoimplementanIsingcouplingJ inEq. (1), sian noise processes that represent the incident vacuum
ij
a part of each DOPO pulse in the fiber ring resonator is fluctuationsfromtheopenportoftheoutputcouplerand
picked-offandfedintoanopticalphasesensitiveamplifier the pump field fluctuation for in-phase and quadrature-
(PSA), followedby opticaldelay lines with intensity and phasecomponents,respectively. Thevacuumfluctuation
phase modulators. Using such N 1 optical delay lines, ofthesignalchannelcontributestothe 1/2termandthe
(arbitrary) i-th pulse can be coup−led to (arbitrary) j-th quantum noise of the pump field contributes to c2 +s2
pulse with a coupling coefficient J . Such an all-optical in the square-rootbracket in (3) and (4).
ij
coupling scheme has been demonstrated for N = 4 and When the i-th signal pulse is incident upon the out-
N =16 CIMs [29, 32]. putcoupler,the output-coupledfieldandremainingfield
3
FIG.1. AcoherentIsingmachinebasedonthetime-divisionmultiplexedDOPOwithmutualcouplingimplementedbyoptical
delaylines. Apartofeachpulseispickedofffromthemaincavitybytheoutputcouplerfollowed byanopticalphasesensitive
amplifier (PSA) which amplifies the in-phase amplitude c˜i of each DOPO pulse. The feedback pulses, which are produced
by combining the outputs from N −1 intensity and phase modulators, are injected back to the main cavity by the injection
coupler.
inside a cavity are written as (8)togetherwiththepumpnoise. Weconductedthe nu-
merical simulation of the coupled CSDE (8) to evaluate
f
c =√Tc √1 T i (5) the performance of the CIM.
i,out i
− − A
s
f
c =√1 Tc +√T i , (6)
i,cavity − i A III. BENCHMARK STUDIES ON MAX-CUT
s
PROBLEMS
where T is the power transmissioncoefficient of the out-
put coupler and f is the incident vacuum fluctuation
i A. MAX-CUT problems on cubic graphs
fromthe openport of the coupler. The out-coupledfield
and the signal amplitude after PSA can be
The MAX-CUT problem on cubic graphs, in which
c 1 T f each vertex has exactly three edges, is called MAX-
i,out i
c˜i ≡ √T =ci−r −T As. (7) CUT-3 problem and also belongs to NP-hard class [35].
The smallest simple MAX-CUT-3 problem is defined on
From these out-coupled pulse stream, the intensity and the complete graph K with four vertices and six edges
4
phase modulators placed in the N 1 delay lines pro- withidenticalweightJ = 1,whereanti-ferromagnetic
− ij −
duce the mutual coupling pulse jξijc˜j, which is actu- couplings have frustration so that the ground states
allyaddedtothei-thsignalpulsePbyaninjectioncoupler. are highly degenerate. The solution to this problem
Here,ξij istheeffectivecouplingcoefficientfromthej-th are the set of two-by-two cuts, which contains six de-
pulse to the i-th pulse, determined by the transmission generate ground states of the Ising Hamiltonian, i.e.,
coefficient√T′oftheinjectioncoupler. Inthehighlydis- , , , , , .
{|↑↑↓↓i |↑↓↑↓i |↑↓↓↑i |↓↑↑↓i |↓↑↓↑i |↓↓↑↑i}
sipativelimitofamutualcouplingcircuit, suchasinour Figure 2 shows the time evolution of c (i = 1,...,4)
i
scheme, we can use the CSDE supplemented with the when p = 1.1 and ξ = 0.1. A correct solution spon-
−
noisy coupling term. Since the transmission coefficient taneously emerges after several tens of round trips. The
√T′oftheinjectioncouplershouldbemuchsmallerthan statistics of obtaining different states against 1000 ses-
one, we do not need to consider any additional noise in sionsofsucha numericalsimulationareshowninFig. 3,
the injected feedback pulse. The CSDE (3) can be now inwhichsixdegenerategroundstatesappearwithalmost
rewritten to include the mutual coupling terms equal probabilities with no errors found.
dc =[( 1+p c2 s2)c + ξ c˜ ]dt
i − − i − i i ij j
Xj B. Many-body interaction problem
1 1
+ c2+s2+ dW . (8)
A r i i 2 i If the interaction is not a standard two-body Ising in-
s
teraction type but rather a four-body interaction such
The summation in Eq. (8) represents the quantum as
measurement-feedback term including the measurement
= J σ σ σ σ , (9)
error given by Eq. (7). The vacuum fluctuation coupled H − 1234 1 2 3 4
to the i-th pulse in the output coupler is already taken where J R, the coupled field into the i-th pulse is
1234
∈
into account in the last term of right-hand side of Eq. no longer given by ξ c˜ but by ξc˜c˜c˜ (j,k,l = i).
j ij j j k l 6
P
4
5, in which eight degenerate ground states are obtained
1.0
ci 1.10 with no errors found.
e
plitud 0.5 11..0005 1.5
m ci
Osignala -00..50 000...899505 ABC mplitude 01..50
P a 0.0
DO -1.0 0.80 33 34 35 36 nal
sig -0.5
O
0 20 40 60 80 100 P -1.0
O
NormalizedtimeHNumberofroundtripsL D -1.5
FIG.2. NormalizedDOPOsignalamplitudesasafunctionof 0 50 100 150 200
normalized time (in unit of cavity round trip numbers) for a Numberofroundtrips
N =4simpleMAX-CUT-3problem. Each colorcorresponds
tothefourdifferentDOPOsindexedwith i=1,...,4. Small FIG.4. Normalized DOPOpulseamplitudesci (i=1,...,4)
window isenlarged toindicatethestatusof signal amplitude under the interaction between four-body Ising coupling ex-
inside a cavity at three components (as in Fig. 1); A: OPA pressed byEq. (9).
gain medium (PPLN waveguide), B: out-coupler, and C: in-
jection coupler of the mutual coupling pulse. The two flat
regions between B and C and between C and A are the pas-
sive propagation in a fiber. 0.14
0.12
¯
0.10 ¯ ¯ ¯
0.15 ¯ y ¯ ¯ ¯
¯ ¯ ¯ bilit 0.08 ¯¯ ¯ ¯¯
a
bility 0.10 ¯¯ ¯ ¯¯ ¯¯ ¯ Prob 00..0046 ¯¯ ¯ ¯¯¯¯ ¯
a
ob ¯ ¯¯ ¯¯ ¯¯¯
Pr 0.05 ¯ ¯ ¯¯ ¯ 0.02 ¯¯¯ ¯¯¯¯
¯¯ ¯¯¯ 0.00
¯¯¯ ¯¯¯¯ SpinconfigurationsΣ
0.00
SpinconfigurationsΣ FIG. 5. Distribution of output spin configurations in 1000
trialsofnumericalsimulationagainstafour-bodyIsingmodel
FIG. 3. Distribution of output spin configurations in 1000 of N = 4. All trials were successful to find one of the eight
trialsofnumericalsimulations against asimpleMAX-CUT-3 degenerate ground states.
problem of graph order N = 4. All trials were successful to
find oneof thesix degenerate ground states.
C. Algorithm description
In this case, the CSDE (8) can be rewritten to include
the four-body coupling term
In this subsection, we will review the four representa-
dc =[( 1+p c2 s2)c +ξc˜ c˜ c˜]dt tive approximation algorithms for MAX-CUT problems.
i − − i − i i j k l
The Goemans-Williamson algorithm (GW) based on
1 1
+ c2+s2+ dW . (10) SDP has a 0.87856-performance guarantee for NP-hard
A r i i 2 i
s MAX-CUT problems [14]. It achieves the optimal ap-
When the four-body coupling coefficient J proximation ratio for MAX-CUT problems under the
1234
is 1 (multi-body anti-ferromagnetic coupling), assumptions of P = NP and the unique games conjec-
there− are eight degenerate ground states, i.e., ture[36]. TheSDP6 relaxationofthe originalMAX-CUT
, , , and their inverse spin problemisavector-valuedoptimizationproblemasmax-
c|↑o↑n↑fi↓giur|a↑t↑i↓o↑nis.|↑F↓i↑g↑uire|↓4↑↑sh↑oiws the time evolution of ci imizing 12 i<jwij(1−~vi·~vj), ~vi ∈Sk−1,whereSk−1 is
(i = 1,...,4) when p = 1.1 and ξ = 0.1. One of the a unit sphePre in Rk and k N (or #V: number of ver-
− ≤
eight degenerate ground states emerges spontaneously tices). Thereexistpolynomialtimealgorithmstofindthe
after several tens of round trips. The statistics of optimal solution of this relaxation problem (with error
observingdifferentstatesin1000independentsessionsof ε>0), and its value is commonly called the SDP upper
the numerical simulation of Eq. (10) are shown in Fig. bound. A final solution to the originalMAX-CUT prob-
5
lem is obtained by projecting the solution vector sets to antee [21], and SG3 is a modified version of it [22]. In
randomly chosen one-dimensionalEuclidean spaces (i.e., this algorithm, nodes V are divided into two disjoint
dividing the sphere by random hyperplanes). subsets S ,S sequentially. For each iterative pro-
1 2
{ }
There are three types of computational complexities cess, the node with the maximum score is selected, and
of the best-known algorithms for solving the SDP relax- it is put into either set S1 or S2 so as to earn larger
ationproblem. IfagraphwithN verticesandmedgesis cuts. Here, the score function of SG3 is defined as
regular,theSDPproblemcanbeapproximatelysolvedin xi = | j∈S1wij − j∈S2wij| (i = 1,...,N). It stops
almostlineartimeasO˜(m)=O(m(logN)2ε−4)usingthe when aPll the edges aPre evaluated to calculate the score
matrix multiplicative weights method [37], where ε rep- function, thus SG3 scales as O(m).
resents the accuracy of the obtained solution. However, The power of breakout local search (BLS) appears in
slower algorithms are required for general graphs. If the the benchmark result for G-set graphs [23]. It updated
edgeweightsofthegraphareallnon-negative,thefastest almosthalfofthebestsolutionsinG-setwiththespecial-
algorithm runs in O˜(Nm) = O(Nm(logN)2ε−3) time ized data structure for sorting and dedicated procedure
based on the Lagrangian relaxation method [38]. For to escape from local minima. The algorithm is combi-
graphswithbothpositiveandnegativeedgeweights,the nation of steepest descent and forced spin flipping: after
SDPproblemiscommonlysolvedusingtheinterior-point being trapped by a local minima as a result of steep-
method, which scales as O˜(N3.5) = O(N3.5log(1/ε)) est descent procedure, three types of forced spin flipping
[39]. Besides, low rank formulation of SDP is effec- (single, pair, and random) are probabilistically executed
tive when the graph is sparse [40–42]. In our compu- according to the vertex influence list (i.e., which vertex
tational experiments, the COPL SDP based on the inte- will increase the number of cut most when it’s flipped)
rior point method was used as a solver for MAX-CUT on each subset of partition.
problems [43]. The SDP upper bound U and the so- These algorithms are coded in C/C++ and run on a
SDP
lution C were obtained using the following parame- single thread of a single core on a Linux machine with
GW
ters: interior point method was used until the relative two 6-core Intel Xeon X5650 (2.67 GHz) processors and
gap r = 1 P /D reached 10−3, where P and 94GBRAM.TheCIMissimulatedbasedonthecoupled
gap obj obj obj
−
D aretheobjectivefunctionsoftheprimalanddualof CSDE(8)onthesamemachine. Notethatthecomputa-
obj
the SDP problem, respectively [44]. Random projection tion time of CIM does not mean the simulation time on
onto one-dimensional space was executed N times. the Linux machine but corresponds to the actual evolu-
tion time of a physical CIM.
For many practical applications, heuristic algorithms
aremoreconvenienttouse,sincetheGWalgorithmgen-
erallyrequireslongcomputationtimeO˜(N3.5). Metropo-
lis et al. introduced a simple algorithmthat can be used D. Computational accuracy on G-set instances
to provideanefficientsimulationofacollectionofatoms
in equilibrium at a given temperature [45]. Kirkpatrick The performance of a CIM with DOPO network
et al. applied the algorithm to optimization problems was tested on the NP-hard MAX-CUT problems on
by replacing the energy of the atomic system to the cost sparse graphs, so-called G-set [47]. These test instances
function of optimization problems and using spin config- wererandomlyconstructedusingamachine-independent
urations σ, which is called the simulated annealing algo- graphgenerator written by G. Rinaldi, with the number
rithm (SA) [15]. Ineachstep ofthis algorithm,a system of vertices ranging from 800 to 20000,edge density from
isgivenwitharandomspinflipandthe resultingchange 0.02%to 6%, and topologyfrom random,almostplanar,
∆E in the energy is computed. If ∆E 0, the spin to toroidal.
≤
flip is always accepted, and the configuration with the The output cut values of running the CIM, SA, GW,
flippedspinisusedasthestartingpointofthenextstep. andthebestknownsolutionssofarwecouldfind[23,48]
If ∆E > 0, the spin is treated probabilistically, i.e., the forsomeofG-setgraphsaresummarizedin TableI.The
probability that the new configuration will be accepted results for CIM are obtained in 50 ms, which correspond
is P(∆E) = exp( ∆E/kBT) with a control parameter to the performanceofanexperimentalsystemafter 5000
−
of system temperature T. This choice of P(∆E) results DOPO cavity roundtrips. The best result andensemble
in the system evolving into an equilibrium Boltzmann averagevalue for 100 trials are shown. Here, the param-
distribution. Repeating this procedure, with the tem- eters are set to be p = 1.6, ξ = 0.06 and the coupling
perature T gradually lowered to zero in sufficiently long constant ξ = ξw / k is no−rmalized by the square
ij ij
time, leads spins σ to convergence to the lowest energy root of the graph averpagheidegree k . The hysteretic op-
state. In practical case, with the finite time, the anneal- timization method, in which thehswiinging and decaying
ing schedule affects the quality of output values. Here in Zeeman term that flips the signal amplitude (spin) back
our numerical simulations, the temperature was lowered and forth, is implemented four times after 10 ms initial
according to the logarithmic function [46]. Note that 1 freeevolution[49]. Eachhystereticoptimizationtakes10
Monte Carlo step corresponds to N trials of spin flip. ms so that the total search takes 50 ms. The result of
SahniandGonzalezconstructedagreedyalgorithmfor SA is also obtained in 50 ms for each graph. For GW,
MAX-CUT problems, which has 1/2-performance guar- thecomputationtimerangedbetween2.3sand1.1 105
×
6
s, depending on N. The best outputs of the CIM were time of CIM also depends on the system configuration
1.62 0.58 % better than GW but 0.38 0.40 % worse and can be made faster when we use the higher clock
± ±
than SA, and CIM found better cut against GW except frequency. In this sense, the ratios between time of CIM
for a toroidal graph (g50) and a disconnected random and that of the other algorithms are arbitrary. Thus we
graph (g70). should study the computation time scaling as a function
As the sizeofoptimizationproblemsincreases,the av- of the problem size.
erage accuracy is important for practical applications. Figure 7 (a) shows the computation time versus prob-
Table I shows that for all G-set graphs, the average ac- lem size (number of vertices). The computation time is
curacy in 100 trials is 0.94148 to the SDP upper bound defined as the CPU time to solve a given MAX-CUT
USDP, i.e., the CIM can find a cut value larger than problem in complete graph for GW; as the CPU time to
0.94148oftheoptimalvaluefortheMAX-CUTproblems reach the same accuracy as GW for SA, SG3, and BLS;
on average, whereas the average accuracy of the GW is andasthetimeestimatedbythe(numberofroundtrips)
0.93025 and that of the SA is 0.94692. Note that USDP (cavity round trip time) to obtain the same accuracy
is always greater than or equal to the optimal value for ×as GW for CIM. The preparation time needed to input
each MAX-CUT problem. J into the computing system, i.e., the graph I/O time,
ij
is not included. For complete graphs of N 20000, the
≤
CIM exhibits a problem-size independent computation
E. Computation time on MAX-CUT problems time of less than 10−3 s if we assume the fixed cavity
circulation frequency of 100 kHz and pulse interval of
10 cm. This means the target accuracy is obtained in
In the previous section, the running time of CIM and
the constant number of round trips. It indicate that the
SA are fixed to estimate the computational accuracy.
computation time of CIM is determined by the turn-on
These two algorithms explored the solutions as good as
delay time of the DOPO network oscillation, which in
possible in 50 ms. Although, if we finish the computa-
turn depends onthe round triptime andthe pump rate.
tion at a certain accuracy,more reasonable computation
time can be defined. Here, the GW solution was used InFigure7(b),thecomputationtimeoftheCIMwith
as the mark of sufficient accuracy because it ensures the different system clock frequency and pulse spacing are
87.856% of the ground states. The CIM and SA then shown (see the Sec. IIA for the definition). Since the
competedthecomputationtimetoreachthesamevalues solutions are obtained in a constant number of cavity
obtained by GW. The time and temperature scheduling round trips, the computation time is pulse spacing in-
parametersoftheSAweresetasfollows: Inversetemper- dependent but linearly depends on the clock frequency,
ature increased with logarithmic function. The number i.e., cavity circulation frequency. The number of pulses
of spin flipping was optimized to be 10l times for some accommodated in the fiber can be changed to vary the
l N, which requires the minimum computation time to pulse spacing under the fixed clock frequency. On the
∈
achieve the same accuracy as with the GW. otherhand, whenthe pulse spacingis fixedandthe fiber
Computational experiments were conducted on fully length is varied, the maximum number of pulses should
connected complete graphs, denoted by K , where the be increased in proportional to the fiber length.
N
number of vertices N ranging from 40 to 20000 and the Going back to the Figure 7 (a), the time complexity
edges are randomly weighted 1. Figure 6 shows the O(N3.5) for the GW is dominated by the interior-point
±
Ising energy in Eq. (1) as a function of running time. method in the Goemans-Williamson algorithm. The SA
Both CIM and SA run stochastically due to quantum seems to scale in O(N2), which indicates that it requires
and thermal noise, respectively, the ensemble average of the number of spin flips to be proportional to N (i.e.,
energies are calculated as follows: For the CIM, the en- constant Monte Carlo steps) to achieve the optimal per-
ergy of all 100 runs was averaged at each round trip. formance. Each spin flip costs a computation time pro-
For the SA, the averaged energy was calculated at each portional to the degree k , where k is equal to N 1
i i
−
point on the time axis with an interpolated value from for all i V in the case of complete graphs. Thus,
real time sampling. The parameters for CIM are chosen the comp∈utation time scales as O(N k ) = O(N2) for
h i
to be p=0.2, ξ = 0.03 (for N =800), and ξ = 0.003 the SA in the complete graphs. Note that CIM and SA
− −
(for N = 4000). In Fig. 6 (a), where N = 800, the GW didn’t always reach the energy obtained by GW for the
achievedanenergyequalto 15624in22.96s, while the graph of N = 40, half of the 100 runs of stochastic al-
CIMandSAreachedthesam−eenergyin6.1 10−4sand gorithms were post-selected to reach that value. SG3
×
in 0.671 s. In Fig. 6 (b), the GW achieved an energy of scalesas O(m)=O(N2) inFig. 7 (a), but the values for
168160 in 6646.25 s, while the CIM reached the same N = 40,160,800 are not shown because it didn’t reach
e−nergy in 5.5 10−4 s and the SA did so in 10.1 s. the accuracy reached by the GW solution. BLS exhibits
×
NotethatthisresultofSAandGWcomesfromaspe- competitiveperformanceagainstSA.Besides,theDOPO
cific computer configuration as mentioned in Sec. IIIC. amplitudesinCIMevolveasinFig. 8(a)whenN =800.
There is room for an improvement in the computation The distribution of computation time for 100 randomly
time in constant factor due to cases like using faster weighted complete graphs of N = 800 is also shown in
CPUs or parallelized codes. Similarly, the computation Fig. 8 (b).
7
TABLEI.PerformanceofthecoherentIsingmachine,simulatedannealingandGoemans-Williamson SDPalgorithminsolving
theMAX-CUTproblemsonsparseG-setgraphs. #V isthenumberN ofverticesinthegraph,#E isthenumbermofedges,
USDP is the optimal solution to the semidefinite relaxation of the MAX-CUT problem when the duality gap reaches to 10−3,
and Cbest is the best known result so far we could find. CGW is the best solution obtained by N projections after SDP. CSA
and hCSAi are the best and average values obtained by SA in 100 trials of 50 ms. CCIM and hCCIMi are the best and average
values in CIM in 100 runs of 50 ms (= 5000 DOPO cavity round trips), respectively. To make comparisons with each other,
every cutvalueC generated from each algorithm is normalized according to (C+Eneg)/(USDP+Eneg),where Eneg ≥0 isthe
numberof negative edges. In thebottom of this table, theaverage and worst values of all 71 G-set graphs are shown.
Graph #V #E USDP Cbest CGW CSA hCSAi CCIM hCCIMi
g1 800 19176 12083 0.9620 0.9457 0.9620 0.9597 0.9614 0.9570
g6 800 19176 2656 0.9607 0.9448 0.9606 0.9592 0.9601 0.9559
g11 800 1600 629 0.9540 0.9327 0.9526 0.9478 0.9455 0.9370
g14 800 4694 3191 0.9602 0.9336 0.9580 0.9544 0.9514 0.9472
g18 800 4694 1166 0.9500 0.9282 0.9492 0.9439 0.9434 0.9372
g22 2000 19990 14136 0.9450 0.9191 0.9445 0.9409 0.9405 0.9361
g27 2000 19990 4141 0.9435 0.9174 0.9422 0.9400 0.9390 0.9356
g32 2000 4000 1567 0.9559 0.9272 0.9508 0.9478 0.9424 0.9384
g35 2000 11778 8014 0.9588 0.9292 0.9551 0.9523 0.9471 0.9438
g39 2000 11778 2877 0.9464 0.9226 0.9431 0.9399 0.9364 0.9318
g43 1000 9990 7032 0.9471 0.9292 0.9471 0.9439 0.9458 0.9396
g48 3000 6000 6000 1.0000 1.0000 1.0000 0.9919 1.0000 0.9747
g51 1000 5909 4006 0.9606 0.9333 0.9583 0.9544 0.9506 0.9468
g55 5000 12498 11039 0.9325 0.9006 0.9264 0.9215 0.9193 0.9160
g57 5000 10000 3885 0.9561 0.9237 0.9496 0.9473 0.9419 0.9384
g59 5000 29570 7312 0.9440 0.9148 0.9376 0.9356 0.9308 0.9288
g60 7000 17148 15222 0.9313 0.8989 0.9231 0.9201 0.9191 0.9152
g64 7000 41459 10466 0.9440 0.9143 0.9347 0.9324 0.9320 0.9299
g67 10000 20000 7744 0.9554 0.9215 0.9480 0.9459 0.9411 0.9388
g70 10000 9999 9863 0.9674 0.9633 0.9523 0.9479 0.9515 0.9482
g81 20000 40000 15656 0.9551 0.9195 0.9187 0.9125 0.9393 0.9376
Average 0.9540 0.9303 0.9502 0.9469 0.9464 0.9415
Worst 0.9313 0.8989 0.9187 0.9125 0.9185 0.9149
HaL HbL
0 0
H -5000 H -50000
y y
g g
ner ner -100000
e -10000 e
g g
n n
Isi CIM SA GW Isi -150000 CIM SA GW
-15000
10-5 10-4 0.001 0.01 0.1 1 10 10-5 0.001 0.1 10 1000
TimeHsL TimeHsL
FIG. 6. Performance comparison of CIM, SA, and GW in solving complete graphs (a) K800 and (b) K4000, where each edge
was randomly weighted ±1. Each time, bundle of curves depicted average energy (solid line) ± standard deviations (dashed
line) in 100 runs. Dotted line is the accuracy obtained by GW algorithm, whose computation time is shown with dot. The
number of spin flip in SA algorithm were 105 for K800 and 106 for K4000, respectively, to optimize the computation time. SA
andGWarerunningon2.67GHzIntelXeonwithoutparallelization andCIMhasthecavitycirculation frequencyof100kHz.
8
HaL HbL
HLutationtimes 1110001046 ìà ìòôà BSGSGLWAS3 ìòôà HLutationtimes 00.00.00.1111 11000kkHHzzcclloocckk 110ccmmmiiinnnttteeerrrvvvaaalll
p ì p m
Com 0.01ìæà ôæà ô CIM æ Com 10-4 1MHzclock 1
10-4ô æ æ 10-5
40 160 800 4000 20000 100 1000 104 105 106 107
ProblemsizeN ProblemsizeN
FIG. 7. (a) Computation time of coherent Ising machine (CIM) empirically scales as O(1), while that of simulated annealing
(SA),andGoemans-Williamson SDP(GW)fittedwelltolinesindicatingO(N2)andO(N3.5),respectively. Computationtime
of SG3 and BLS also fit to O(N2) with a difference of constant factor. The computation time of CIM, SA, SG3, and BLS is
defined as the time to reach the same accuracy achieved by GW (without I/O time). Data points of CIM, SA, SG3, and BLS
werecalculatedbyaveraging100runs. OnlyforN =800case,100randomly±1weightedcompletegraphsaregenerated. The
results for N =800 show the average computation time for 100 graphs with error bars indicating standard deviations (which
aredistributedasinFig. 8(b)). NotethatthecomputationtimeofCIMisevaluatedby(numberofroundtrips)×(10µs)as
in Sec. IIA. (b) The computation time of CIM when the DOPO cavity circulation frequency {10 kHz, 100 kHz, 1 MHz} and
pulse interval {10 cm (solid line), 1 cm (dashed line), 1 mm (dotted line)} are changed. Computation time of CIM scales as
O(N) if we vary thefiberlength proportional to N with fixedpulse interval.
HaL HbL
2
0.5
Žci
e
ud 1 0.4
plit y
am 0 bilit 0.3 CIM
Osignal -1 Proba 0.2 SA
P
O 0.1
D GW
-2
0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.001 0.01 0.1 1 10
TimeHmsL ComputationtimeHsL
FIG.8. DetaileddataforN =800randomlyweightedcompletegraphs(appearedinFig. 6(a)). (a)TimeevolutionofDOPO
amplitudes of the graph of N =800. Only 100 pulses out of N =800 signal pulses are shown. (b)Histogram of computation
times for solving 100 randomly weighted complete graphs of size N =800 bythe CIM, SA,and GW.
Computation time for the random graphs in G-set in- edge weights so that we must use the slowest interior
stancesisalsostudied. Herethesubsetofgraphsinwhich point method.) Then the computation time is almost
MAX-CUT problems can be solved in polynomial time constantforbothSAandCIM.Thecomputationtimeof
(i.e., planar graphs, weakly bipartite graphs, positive SA with constant Monte Carlo step is expected to scale
weighted graphs without a long odd cycle, and graphs O(N k ) O(1) (here for the random graphs in G-set,
withintegeredgeweightboundedbyN andfixedgenus) k h Oi(N∼−1.09)). The computation time of CIM here
h i ∼
areexcluded. TheexecutiontimeofCIMisevaluatedun- is governed by a turn-on delay time of the DOPO net-
derthemachinespecdescribedinSec. IIAwithp=0.2, work to reach a steady state oscillation condition, which
ξ = 0.06, and ξ = ξw / k . Again, the compu- is constant for varying values of N as mentioned above
ij ij
− h i
tation time of SA and CIM ips the actual time (without [25].
graphfileI/O)toobtainthesameaccuracyofsolutionas
GW. Figure 9 shows the computation time as functions
oftheproblemsizeN. Thecomputationalcostofinterior
point method dominates the GW algorithm. (Note that
G-set contains graphs with both positive and negative
9
105 GW age and found better cut compared to the GW for 69
Ls ì ì outof71graphsinG-set. Thecomputationtime forthis
He 1000 ì sparse graph set, including few sessions of hysteretic op-
tim ì à timization, is estimated as 50 ms. The time scaling was
on 10 ì ì also tested on complete graphs of number of vertices up
ati to 2 104 and number of edges up to 108. The results
mput 0.1 à à à à àSA æ iimngpilny×athfiaxtedCIsMystaecmhicelvoecskefmrepqiuriecnaclyly,ic.eo.n,stthaenfitxteimdecasvciatly-
o
C
CIM circulation frequency (fiber length), while SA, SG3, and
0.001 æ æ æ æ æ BLS scale as O(N2) and GW scales as O(N3.5). Those
800 2000 5000 700010000
results suggest that CIM may find applications in high-
ProblemsizeN speed computation for various combinatorial optimiza-
tion problems, in particular for temporal networks.
FIG. 9. Computation time of coherent Ising machine, sim- The present simulation results do not mean that the
ulated annealing algorithm, and Goemans-Williamson SDP CIMcangetareasonablyaccuratesolutionbyaconstant
algorithm on random graphsin G-set instances. Thecompu-
time for arbitrary large problem size. As mentioned al-
tation time for CIM and SA is defined as the time to reach
ready, in the CIM based on a fiber ring resonator, the
thesameaccuracyachievedbyGW(withoutI/Otime). Note
number of DOPO pulses is determined by the the fiber
that the computation time of CIM is evaluated by (number
length and the pulse spacing. In order to implement
of round trips) × (10 µs) as in Sec. IIA.
2 105 and 2 106 DOPO pulses in the 20 km fiber
× ×
ring cavity, we must use a pulse repetition frequency to
2 GHz and 20 GHz, respectively. This is a challenge for
both optical components and electronic components of
IV. SUMMARY AND DISCUSSION
CIM, but certainly within a reach in current technolo-
gies.
The potential for solving NP-hard problems using a ACKNOWLEDGMENTS
CIM was numerically studied by conducting computa-
tional experiments using the MAX-CUT problems on TheauthorswouldliketothankK.Inoue,H.Takesue,
sparseG-setgraphsandfullyconnectedcompletegraphs K. Aihara, A. Marandi, P. McMahon, T. Leleu, S. Ta-
of order up to 2 104. With the normalized pump rate mate, K. Yan, Z. Wang, and K. Takata for their useful
×
and coupling coefficient p=0.2 and ξ = 0.06,the CIM discussions. This project is supported by the ImPACT
−
achieved a good approximation rate of 0.94148 on aver- programof the Japanese Cabinet Office.
[1] R. M. Karp, in Complexity of Computer Computations, Szegedy, J. ACM 45, 501 (1998).
edited by R. E. Millera and J. W. Thatcher (Plenum, [13] J. H˚astad, J. ACM 48, 798 (2001).
NewYork, 1972), pp.85–103. [14] M.X.GoemansandD.P.Williamson, J.ACM42,1115
[2] M.M´ezard, G.Parisi, andM.Virasoro, Spin Glass The- (1995).
ory and Beyond (World Scientific, Singapore, 1987). [15] S. Kirkpatrick, C. D. Gelatt Jr., and M. P. Vecchi, Sci-
[3] D.B.Kitchen,H.Decornez,J.R.Furr,andJ.Bajorath, ence 220, 671 (1983).
Nat.Rev. DrugDiscov. 3, 935 (2004). [16] T. Kadowaki and H. Nishimori, Phys. Rev. E 58, 5355
[4] H.Nishimori,Statistical Physics of Spin Glasses and In- (1998).
formation Processing (Oxford University Press, Oxford, [17] G. E. Santoro, R. Martonˇ´ak, E. Tosatti, and R. Car,
2001). Science 295, 2427 (2002).
[5] F. Barahona, J. Phys. A 15, 3241 (1982). [18] E.Farhi,J. Goldstone,S.Gutmann,J.Lapan,A.Lund-
[6] M. R. Garey and D. S. Johnson, Computers and In- gren, and D. Preda, Science 292, 472 (2001).
tractability: A Guide to the Theory of NP-Completeness [19] W. van Dam, M. Mosca, and U. V. Vazirani, in Pro-
(Freeman, San Francisco, 1979). ceedingsofthe42ndIEEESymposiumonFoundationsof
[7] G. I. Orlova, Y. G. Dorfman, Eng. Cybern. 10, 502 Computer Science (IEEE Computer Society, Los Alami-
(1972). tos, CA, 2001), pp. 279–287.
[8] F. Hadlock, SIAMJ. Comput. 4, 221 (1975). [20] D. Aharonov, W. van Dam, J. Kempe, Z. Landau, S.
[9] M. Gr¨otschel and W. R. Pulleyblank, Oper. Res. Lett., Lloyd, and O. Regev, SIAMJ. Comput. 37, 166 (2007).
1, 23 (1981). [21] S. Sahniand T. Gonzalez, J. ACM 23, 555 (1976).
[10] M.Gr¨otschelandG.L.Nemhauser,Math.Program.,29, [22] S.Kahruman,E.Kolotoglu,S.Butenko,andI.V.Hicks,
28 (1984). Int.J. Comput. Sci. Eng. 3, 211 (2007).
[11] A.Galluccio,M.Loebl,andJ.Vondr´ak,Math.Program., [23] U.BenlicandJ.-K.Hao,Eng.Appl.Artif.Intel.26,1162
90, 273 (2001). (2013).
[12] S. Arora, C. Lund, R. Motwani, M. Sudan, and M. [24] S. Utsunomiya, K. Takata, and Y. Yamamoto, Opt. ex-
10
press 19, 18091 (2011). 2007) pp.227–236.
[25] K. Takata, S. Utsunomiya, and Y. Yamamoto, New J. [38] P.KleinandH.-I.Lu,inProceedings ofthetwenty-eighth
Phys.14, 013052 (2012). annual ACMsymposium onTheoryofcomputing(ACM,
[26] K. Takata and Y. Yamamoto, Phys. Rev. A 89, 032319 1996) pp.338–347.
(2014). [39] F. Alizadeh, SIAMJ. Optim. 5, 13 (1995).
[27] S. Utsunomiya, N. Namekata, K. Takata, D. Akamatsu, [40] M. Yamashita, K. Fujisawa, M. Fukuda, K. Kobayashi,
S. Inoue, and Y. Yamamoto, Opt. express 23, 6029 K.Nakata,andM.Nakata,inHandbook onsemidefinite,
(2015). conicandpolynomialoptimization,editedbyM.F.Anjos
[28] Z. Wang, A. Marandi, K. Wen, R. L. Byer, and Y. Ya- and J. B. Lasserre (Springer, 2012) pp. 687–713.
mamoto, Phys.Rev.A 88, 063853 (2013). [41] L. Grippo, L. Palagi, M. Piacentini, V. Piccialli, and G.
[29] A. Marandi, Z. Wang, K. Takata, R. L. Byer, and Y. Rinaldi, Math. Program. 136, 353 (2012).
Yamamoto, Nat. Photonics 8, 937 (2014). [42] K. Fujisawa, T. Endo, Y. Yasui, and H. Sato, in Pro-
[30] P. D. Drummond and C. W. Gardiner, J. Phys. A 13, ceedings of the 28th IEEE International Parallel & Dis-
2353 (1980). tributed Processing Symposium (IEEE, 2014) pp. 1171–
[31] K.Takata,A.MarandiandY.Yamamoto,Phys.Rev.A 1180.
92, 043821 (2015). [43] Y. Ye, “Computational optimization laboratory”,
[32] K. Takata, Ph.D. dissertation, the University of Tokyo, http://web.stanford.edu/~yyye/Col.html (1999).
Tokyo(2015). [44] S.J. Benson, Y.Ye,andX.Zhang, SIAMJ. Optim.10,
[33] P.Kinslerand P.D.Drummond,Phys.Rev.A 43, 6194 443 (2000).
(1991). [45] N.Metropolis, A.W.Rosenbluth,M. N.Rosenbluth,A.
[34] P. D. Drummond, K. J. McNeil, and D. F. Walls, Opt. H. Teller, and E. Teller, J.Chem. Phys.21, 1087 (1953).
Acta28, 211 (1980). [46] B. Hajek, Math. Oper. Res. 13, 311 (1988).
[35] E. Halperin, D. Livnat, and U. Zwick, in Proceedings of [47] C. Helmberg and F. Rendl, SIAM J. Optim. 10, 673
thethirteenthannualACM-SIAMsymposiumonDiscrete (2000).
algorithms (SIAM,2002) pp. 506–513. [48] T. Ikuta, H. Imai and Y. Yano, “Solving MAX-
[36] S.Khot,G.Kindler,E.Mossel,andR.O’Donnell,SIAM CUTbenchmarkbyoptimization solvers”(inJapanese),
J. Comput. 37, 319 (2007). (RIMS Preprint, Kyoto, Japan, 2015) 1941-09.
[37] S. Arora and S. Kale, in Proceedings of the thirty-ninth [49] G. Zar´and, F. Pazmandi, K. F. Pal, and G. T. Zimanyi,
annualACMsymposium onTheory ofcomputing(ACM, Phys. Rev.Lett. 89, 150201 (2002).