J. Phys. Soc. Jpn. 88, 061010 (2019) [10 Pages]
SPECIAL TOPICS: Quantum Annealing: Recent Development and Future Perspectives

Application of Ising Machines and a Software Development for Ising Machines

+ Affiliations
1Recruit Communications Co., Ltd., Chuo, Tokyo 104-0054, Japan2Green Computing Systems Research Organization, Waseda University, Shinjuku, Tokyo 162-0042, Japan3JST PRESTO, Kawaguchi, Saitama 332-0012, Japan

An online advertisement optimization, which can be represented by a combinatorial optimization problem is performed using D-Wave 2000Q, a quantum annealing machine. To optimize the online advertisement allocation optimization, we introduce a generalized version of the Markowitz mean-variance model which is a basic model of portfolio optimization. The obtained optimization performance using D-Wave 2000Q is higher than that using the greedy method which is a conventional method. Additionally, to conveniently use Ising machines including a quantum annealing machine, new software called PyQUBO is developed. The first half of the paper gives a review of several combinatorial optimization problems and how to represent them using the Ising model or the quadratic unconstrained binary optimization (QUBO) form. We show the results of the online advertisement allocation optimization and the explanation of PyQUBO in the last half of the paper.

©2019 The Author(s)
This article is published by the Physical Society of Japan under the terms of the Creative Commons Attribution 4.0 License. Any further distribution of this work must maintain attribution to the author(s) and the title of the article, journal citation, and DOI.
1. Introduction

Ising machines have actively developed since the birth of commercial quantum annealing machine, D-Wave.16) Recently developed Ising machines specialize in searching better solutions to combinatorial optimization problems represented by an Ising model or a quadratic unconstrained binary optimization (QUBO) form.7,8) Some earlier studies reported the advantages of using Ising machines for specific combinatorial optimization problems.5,912) However, it is still unclear whether Ising machines are more efficient than conventional computing technology. Thus, it is necessary to do physics-based theoretical studies and develop hardware and software for various applications.

Combinatorial optimization problems exist in our daily life, e.g., optimization of shift-planning, traffic flow optimization, and delivery optimization. As the amount of data traffic increases due to the development of more sensing devices, better high-performance computing technology for combinatorial optimization processing is required. Additionally, new combinatorial optimization problems are expected to emerge in the internet of things (IoT) society because of improved global network connectivity. Thus, Ising machines have attracted much attention as a new computing technology for high-performance combinatorial optimization processing.

Kadowaki and Nishimori in 1998 proposed quantum annealing as an alternative to simulated annealing.13) They investigated the dynamic behavior of the Ising model with the time-dependent transverse field using the Schrödinger equation. In the quantum annealing, we prepare the ground state of the Ising model with a large transverse field as the initial state. After that, we gradually decrease the transverse field to zero, and we observe the final state. When the transverse field is decreased with a sufficiently slow schedule so that adiabatic condition is satisfied, we obtain the ground state. Reference 13 investigates the difference between the dynamic behavior of simulated annealing and that of quantum annealing. The authors concluded that the quantum annealing could obtain the ground state with higher probability than the simulated annealing when we use the same scheduling function which is a faster schedule that satisfies the Geman–Geman theorem for the simulated annealing.14)

In the original version of quantum annealing proposed by Kadowaki and Nishimori,13) the Ising model is used as a classical Hamiltonian, while the time-dependent transverse field is used as a quantum driver Hamiltonian. Quantum Monte Carlo simulation can be used to simulate the Ising model with transverse field by applying the Suzuki–Trotter decomposition,15,16) where the Ising model with transverse field (quantum Ising model) in d-dimension is transformed to the Ising model without transverse field (classical Ising model) in \(d+1\)-dimension. Then, the comparison of the performances of quantum annealing and simulated annealing can be made using Monte Carlo methods or other similar methods such as mean-field approximation based methods.1720) Moreover, the concept of quantum annealing has been extended for a system which cannot be easily represented using the Ising model. For example, Kurihara et al. proposed a method for quantum annealing that can be applied to clustering by introducing a new type of quantum driver Hamiltonian.21) Since then, similar extensions of quantum annealing have been proposed.2224)

Before the birth of commercial quantum annealing machine, theoretical and numerical studies on the performance of quantum annealing were done. After quantum annealing machines are available, the comparisons of experimental quantum annealing performance and that of other methods using a digital computer have been made. Recently, Maezawa et al. proposed and developed a new type of quantum annealing machine.25,26) Additionally, several Ising machines with operating principles not based on quantum annealing have been developed.26)

Toward the practical usage of Ising machines, it is necessary to develop software for Ising machines and discover various applications. A high-level programming language for the D-Wave machine based on Julia,27) Python-based software tools for annealing machines,28) and virtual hardware embedding suite for adiabatic quantum computing29) are examples of software developed recently.

In this paper, we review approaches to use Ising machines. Furthermore, we show our result of an application of Ising machines and software for Ising machines developed by one of the authors (K.T.). The organization of the paper is as follows. In Sect. 2, we briefly review combinatorial optimization problems. In Sect. 3, we explain the relation of the combinatorial optimization problem and the Ising model. In Sect. 4, we show the experimental result of online advertisement allocation optimization using D-Wave 2000Q. In Sect. 5, we introduce software called PyQUBO. In Sect. 6, we provide a summary of this paper.

2. Combinatorial Optimization Problem

Optimization problems pertain to finding conditions that minimize a given function, called cost function, under some given constraints. A mathematical expression of these problems are given by \begin{equation} \mathbf{x}^{*} = \mathop{\mathrm{arg\,min}}\nolimits_{\mathbf{x}} f(\mathbf{x}),\quad \mathbf{x}\in \mathcal{S}, \end{equation} (1) where \(\mathbf{x}\) is the vector representing the decision variables, \(f(\mathbf{x})\) is the cost function that is real-valued, and \(\mathcal{S}\) is the set of decision variables satisfying the given constraints, in other words, the set of feasible solutions. When \(\mathbf{x}\) is a discretized vector, the above optimization problem is called a combinatorial optimization problem.

We assume that equalities and inequalities represent the given constraints. Suppose that \(\mathbf{x}\) is a real-valued vector, the above optimization problems are rewritten as \begin{align} &\mathbf{x}^{*} = \mathop{\mathrm{arg\,min}}\nolimits_{\mathbf{x}} f(\mathbf{x}),\quad \mathbf{x}\in \mathcal{R}^{n}, \end{align} (2) \begin{align} &\begin{cases} g_{\ell}(\mathbf{x}) = 0 &\text{($\ell = 1, \ldots, L$)},\\ h_{m}(\mathbf{x}) \leq 0 &\text{($m = 1, \ldots, M$)}, \end{cases} \end{align} (3) where \(g_{\ell}(\mathbf{x})\) and \(h_{m}(\mathbf{x})\) are called the equality constraints and the inequality constraints, respectively. Here we show a method, called the penalty function method, to perform the optimization only when an equality constraint exists, \(g(\mathbf{x})=0\). In the penalty function method, we consider the following optimization problem without any constraints \begin{equation} \mathbf{x}^{*} = \mathop{\mathrm{arg\,min}}\nolimits_{\mathbf{x}} \{f(\mathbf{x}) + \lambda [g(\mathbf{x})]^{2}\},\quad \mathbf{x}\in \mathcal{R}^{n}. \end{equation} (4) Practically, we start with a small positive value of λ and gradually increase it to obtain solutions satisfying the given constraint. Note that feasible solutions are obtained with high probabilities for a large enough value of λ. Also, the above scheme can be adopted even when two or more equality constraints exist.

Next, for the case when only an inequality constraint exists, \(h(\mathbf{x})\leq 0\). In the penalty function method, we consider the following optimization problem without any constraints \begin{equation} \mathbf{x}^{*} = \mathop{\mathrm{arg\,min}}\nolimits_{\mathbf{x}} \{f(\mathbf{x}) + \lambda \max[h(\mathbf{x}),0]\},\quad \mathbf{x}\in \mathcal{R}^{n}. \end{equation} (5) Similarly to the case when an equality constraint exists, the above scheme can be used even when two or more inequality constraints exist.

Herein, we concentrate on the case where the decision variables are integers, i.e., \(\mathbf{x}\in \mathcal{Z}^{n}\). Hence, we refer to the optimization problem as the combinatorial optimization problem as stated above. Additionally, we focus on the case where the cost function and the given constraints are polynomial functions. The above assumptions are valid for many combinatorial optimization problems, e.g., traveling salesman problem and knapsack problem. The cost functions of the combinatorial optimization problems under the above assumption are represented by the Hamiltonian of the Ising model as will be described in the next section.

3. Ising Model or QUBO Corresponding to Combinatorial Optimization Problems

To perform combinatorial optimization using Ising machines, we should prepare the Hamiltonian of the Ising model or the QUBO which corresponds to the cost function of the combinatorial optimization problem. In this section, we explain methods to construct the Hamiltonian of the Ising model or the QUBO from the given combinatorial optimization problems. Before showing methods to represent the cost function of the combinatorial optimization problem using the Hamiltonian of the Ising model or the QUBO, we explain the relations between the Ising model and the QUBO. The Ising model and the QUBO are defined on an undirected graph. Consider an undirected graph \(G=(V,E)\), where V and E are the sets of vertices and edges on G, respectively. The Hamiltonian of the Ising model on \(G=(V,E)\) is given by \begin{equation} \mathcal{H}(\mathbf{s}) = \sum_{i \in V} h_{i} s_{i} + \sum_{(ij) \in E} J_{ij} s_{i} s_{j},\quad s_{i} = \pm 1, \end{equation} (6) where \(s_{i}\) is the spin at site \(i\in V\), \(h_{i}\) is the magnetic field at site \(i\in V\), and \(J_{ij}\) is the interaction on edge \((ij)\in E\); both of values \(h_{i}\) and \(J_{ij}\) are real-valued. The argument of LHS of Eq. (6), \(\mathbf{s}\), corresponds to the array of spins, i.e., \(\mathbf{s}=(s_{1},s_{2},\ldots,s_{N})\), where N is the number of spins which is equal to the number of vertices \(|V|\).

Let \(x_{i}\) be a binary variable at site i such that \(x_{i}=0,1\). The cost function of 0–1 integer programming problem with only linear and quadratic terms is given by \begin{align} f(\mathbf{x}) &= \sum_{i \in V} a_{i} x_{i} + \sum_{(ij) \in E} b_{ij} x_{i} x_{j} \end{align} (7) \begin{align} &= \sum_{i \in V} \sum_{j \in V} Q_{ij} x_{i} x_{j}. \end{align} (8) The above formulation is called the QUBO and the matrix Q is called the QUBO matrix. The relations of the coefficients \(a_{i}\) and \(b_{ij}\) and the QUBO matrix are given by \begin{align} Q_{ii} &= a_{i}\quad [\forall i \in V], \end{align} (9) \begin{align} Q_{ij} &= \begin{cases} b_{ij} &\text{[$\forall(ij) \in E$, $i \neq j$]}\\ 0 &\text{[$\forall(ij) \notin E$, $i \neq j$]}\end{cases}. \end{align} (10) Using the relation \(x_{i} = (s_{i} + 1)/2\), the transformation from the QUBO to the Ising model can be derived as \begin{align} &h_{i} = \frac{a_{i}}{2} + \sum_{j \in \partial_{i}} \frac{b_{ij}}{2},\quad \forall i \in V, \end{align} (11) \begin{align} &J_{ij} = \frac{b_{ij}}{4},\quad \forall (ij) \in E, \end{align} (12) where \(\partial_{i}\) is the set of nearest neighbor vertices of the vertex i on \(G=(E,V)\).

From the above, when the QUBO represents the cost function, we can reformulate the combinatorial optimization problem using the Ising model. In the penalty function method (see Sect. 2), the equality/inequality constraints which are imposed in the original combinatorial optimization problem are included in another combinatorial optimization problem [see Eq. (4)]. Thus, if we represent the equality/inequality constraints in the original combinatorial optimization problem using only linear and quadratic terms of 0–1 binary variables, we can reformulate the combinatorial optimization problem using the Ising model or the QUBO. In principle, even if k-body (\(k\geq 3\)) interactions among 0–1 binary variables or spins exist in the cost function or the term corresponding to the constraints in the original combinatorial optimization problem, we can reconstruct the Ising model or the QUBO. For example, if we have a Hamiltonian \(\mathcal{H} = xyz\), where x, y, and z are 0–1 binary variables, we can introduce a new auxiliary 0–1 binary variable a to express the product of x and y. To ensure that the relations \(a=xy\) is maintained, we also introduce an additional penalty term \(\text{AND}(a, x, y)\equiv xy - 2a(x + y) + 3a\) to the Hamiltonian. Thus, the final Hamiltonian is \(\mathcal{H} = az +\alpha \text{AND}(a, x, y)\), where α is a parameter of the penalty strength. The software tool named PyQUBO which will be introduced in Sect. 5, can automatically reduce the degrees of freedom of the Hamiltonian by introducing auxiliary variables in this way.

Typical expressions of combinatorial optimization problems

Here, we explain how to express some types of combinatorial optimization problem using the Ising model or QUBO. According to the existence of equality/inequality constraints, we categorize combinatorial optimization problems into three types. References 7 and 8 provide more detailed information.

Case I: No constraints

Here, we consider combinatorial optimization problems without constraints and represent them using the Ising model or QUBO. A typical combinatorial optimization problem without constraints is the number partitioning problem described as follows. Given a set of N positive integers, \(\mathcal{S}=\{n_{1}, n_{2},\ldots, n_{N}\}\), the problem aims to find two disjoint subsets \(\mathcal{S}_{1}\) and \(\mathcal{S}_{2}=\mathcal{S} -\mathcal{S}_{1}\) so as to minimize the difference between the sum of the elements in \(\mathcal{S}_{1}\) and the sum of the elements in \(\mathcal{S}_{2}\). Here, the difference between the sums of both sets is regarded as the cost function. Thus the cost function is given by \(|\sum_{i\in \mathcal{S}_{1}} n_{i} -\sum_{i\in \mathcal{S}_{2}} n_{i}|\). To express the cost function using the Ising model, let \(s_{i}\) be the spin variable such that \(s_{i}=\pm 1\). The spin variable \(s_{i}\) indicates the subset where the number \(n_{i}\) is in, i.e., if \(n_{i}\in \mathcal{S}_{1}\), \(s_{i}=1\), whereas if \(n_{i}\in \mathcal{S}_{2}\), \(s_{i}=-1\). Hence, we can construct the Ising model as follows: \begin{equation} \mathcal{H} = \left(\sum_{i=1}^{N} n_{i} s_{i} \right)^{2}. \end{equation} (13) The Hamiltonian is minimized when the cost function of the original problem is minimized.

Case II: Equality constraints

Here, we consider combinatorial optimization problems with equality constraints and represent them using the Ising model or QUBO. In this section, we show two examples, i.e., the graph partitioning problem and traveling salesman problem. The graph partitioning problem has an equality constraint in the original form of the combinatorial optimization problem. Conversely, in the traveling salesman problem, equality constraints do not appear in the original form of the combinatorial optimization problem but are introduced when we use the QUBO expression.

The graph partitioning problem is described as follows. Consider an undirected graph \(G=(V,E)\), where V and E are the sets of vertices and edges on G, respectively; here, we assume that the number of vertices \(|V|\) is even. The problem aims to find a partition of G into two disjoint subgraphs \(G_{1}=(V_{1},E_{1})\) and \(G_{2}=(V_{2},E_{2})\) so as to minimize the number of edges connecting the two subgraphs under the constraint, \(|V_{1}|=|V_{2}|\). The cost function of this problem is the number of edges connecting two subgraphs \(G_{1}\) and \(G_{2}\), while the equality constraint is \(|V_{1}|=|V_{2}|\).

To represent the cost function using the Ising model, let \(s_{i}\) be the spin variable, such that \(s_{i}=\pm 1\), similarly to the number partitioning problem described in the previous subsection. The spin variable \(s_{i}\) indicates the subgraph where the vertex i is in, i.e., if \(i\in G_{1}\), \(s_{i}=1\) whereas if \(i\in G_{2}\), \(s_{i}=-1\). Hence, we can represent the cost function as \begin{equation} \mathcal{H}_{\text{cost}} = \sum_{(ij) \in E} \frac{1-s_{i}s_{j}}{2}. \end{equation} (14) The argument of the RHS in Eq. (14) represents the number of edges connecting the two subgraphs, since \((1-s_{i}s_{j})/2\) is equal to unity when the vertices i and j are in different subgraphs, otherwise, it is equal to zero. Thus, Eq. (14) is the cost function of the problem. Next, we represent the equality constraint, \(|V_{1}|=|V_{2}|\), using the Ising model. Using the spin variables \(s_{i}\), we can express the constraint as \begin{equation} \mathcal{H}_{\text{constraint}} = \left(\sum_{i \in V} s_{i} \right)^{2}. \end{equation} (15) When \(|V_{1}|=|V_{2}|\), the above Hamiltonian becomes zero, which is the minimum value.

From the above discussion, we prepared the Hamiltonians representing the cost function and the equality constraint. As in the penalty function method [Eq. (4)], we introduce a new Hamiltonian as \begin{equation} \mathcal{H} = \mathcal{H}_{\text{cost}} + \lambda \mathcal{H}_{\text{constraint}}. \end{equation} (16) Generally, the probability of obtaining a feasible solution increases but the cost function becomes increases as the value of λ increases. Thus, λ should be tuned to an appropriate value by repeating the calculation. As will be explained in Sect. 5, by using a Placeholder class of the software called PyQUBO,30) we can save the preparing time of QUBO.

Next, we consider the traveling salesman problem. The traveling salesman problem for a graph \(G=(V,E)\) is described as follows. The weight \(\ell_{ij}\) [\((i,j)\in E\)] is defined by the distance from the vertex i to the vertex j. When the given graph is an undirected graph, \(\ell_{ij} =\ell_{ji}\), otherwise \(\ell_{ij}\neq \ell_{ji}\). Note, however, that the representation which will be shown is robust for either undirected or directed graphs. The problem aims to find the minimum of the sum of the weights with constraints that all vertices are only visited once and that the traveler needs to return to the starting point at the final step. Herein we assume that G is a complete graph.

We can express the cost function as \begin{equation} L = \sum_{t=1}^{N} \ell_{v(t),v(t+1)},\quad v(N+1) = v(1), \end{equation} (17) where N is the number of vertices, i.e., \(N=|V|\), and \(v(t)\) is the vertex where the traveler visits at time step t. Next, we represent the cost function by using the QUBO form. Let \(x_{i,t}\) (\(1\leq i\leq N\) and \(1\leq t\leq N\)) be a 0–1 binary variable. When the traveler visits the vertex i at time step t, \(x_{i,t} = 1\), otherwise \(x_{i,t}=0\). Using the 0–1 binary variables defined above, we can redefine the cost function in the QUBO form as \begin{equation} \mathcal{H}_{\text{cost}} = \sum_{t=1}^{N} \sum_{j=1}^{N} \sum_{i=1}^{N} \ell_{i,j} x_{i,t} x_{j,t+1},\quad x_{i,N+1} = x_{i,1}\ (\forall i). \end{equation} (18) However, the above equation is not enough to express the traveling salesman problem. The following conditions for the 0–1 binary variables need to be satisfied: \begin{align} \sum_{i=1}^{N} x_{i,t} &= 1,\quad \forall t, \end{align} (19) \begin{align} \sum_{t=1}^{N} x_{i,t} &= 1,\quad \forall i. \end{align} (20) The former condition means that the traveler should only be at a single place at each time step t. The latter condition means that the traveler should visit vertex i only once. To add the above conditions to the QUBO, we consider the following Hamiltonian: \begin{align} &\mathcal{H} = \mathcal{H}_{\text{cost}} + \lambda_{1} \mathcal{H}_{\text{constraint 1}} + \lambda_{2} \mathcal{H}_{\text{constraint 2}}, \end{align} (21) \begin{align} &\mathcal{H}_{\text{cost}} = \sum_{t=1}^{N} \sum_{j=1}^{N} \sum_{i=1}^{N} \ell_{i,j} x_{i,t} x_{j,t+1}, \end{align} (22) \begin{align} &\mathcal{H}_{\text{constraint 1}} = \sum_{t=1}^{N} \left(\sum_{i=1}^{N} x_{i,t} - 1 \right)^{2}, \end{align} (23) \begin{align} &\mathcal{H}_{\text{constraint 2}} = \sum_{i=1}^{N} \left(\sum_{t=1}^{N} x_{i,t} - 1 \right)^{2}. \end{align} (24) When the constraints given in Eqs. (19) and (20) are satisfied, Eqs. (23) and (24) become zero, which is the minimum value.

Case III: Inequality constraints

Here, we consider combinatorial optimization problems with inequality constraints and represent them using the Ising model or the QUBO. Here, we explain how to represent the knapsack problem using the QUBO, which has an inequality constraint in the original form of the combinatorial optimization problem.

The knapsack problem is described as follows. There are N items and the weight and price of the i-th (\(1\leq i\leq N\)) item are \(w_{i}\) and \(p_{i}\), respectively. The problem is to find the item set that maximizes the total price P with a constraint that the total weight W is smaller than or equal to the given weight \(\mathcal{W}\). Let \(x_{i}\) be a 0–1 binary variable. When the i-th item is selected, \(x_{i}=1\), otherwise \(x_{i}=0\). Using the 0–1 binary variables \(\{x_{i}\}_{i=1}^{N}\), the total price and weight are respectively given by \begin{align} P &= \sum_{i = 1}^{N} p_{i} x_{i}, \end{align} (25) \begin{align} W &= \sum_{i = 1}^{N} w_{i} x_{i}. \end{align} (26) The knapsack problem aims to find the combination of 0–1 binary variables \(\{x_{i}\}_{i=1}^{N}\) such that P is maximized under the condition \(W\leq \mathcal{W}\). The cost function of the knapsack problem is given by \begin{equation} \mathcal{H}_{\text{cost}} = - \sum_{i=1}^{N} p_{i} x_{i}. \end{equation} (27) The minus sign in the RHS of the above equation is needed to change the maximization problem to a minimization problem.

Next, we consider a method to represent the inequality constraint, \(W\leq \mathcal{W}\). Let \(y_{j}\) (\(1\leq j\leq \mathcal{W}\)) be 0–1 another binary variable. Using the 0–1 binary variables \(\{y_{j}\}_{j=1}^{\mathcal{W}}\), the QUBO form of the inequality constraint is given by \begin{equation} \mathcal{H}_{\text{constraint}} = \left(1 - \sum_{j=1}^{\mathcal{W}} y_{j} \right)^{2}{}+\left(\sum_{j=1}^{\mathcal{W}} j y_{j} - \sum_{i=1}^{N} w_{i} x_{i} \right)^{2}. \end{equation} (28) The above Hamiltonian becomes zero if and only if \(y_{W}=1\) and \(y_{i}=0\) for arbitrary i except for \(i=W\). Thus, the total Hamiltonian corresponding to the knapsack problem is given by \begin{equation} \mathcal{H} = \mathcal{H}_{\text{cost}} + \lambda \mathcal{H}_{\text{constraint}}. \end{equation} (29)

Embedding

In the previous subsection, we explained standard methods to represent typical combinatorial optimization problems by the Ising model or the QUBO depending on whether equality/inequality constraints exist. At the stage of thinking about the representation of the combinatorial optimization problem using the Ising model or the QUBO, we do not take care of the graph structure of Ising machines. Let \(G_{\text{L}}=(V_{\text{L}},E_{\text{L}})\) be the undirected graph on which the Ising model or the QUBO is defined, which is referred to as the logical graph. Here \(V_{\text{L}}\) and \(E_{\text{L}}\) are the sets of vertices and edges on the logical graph \(G_{\text{L}}\), respectively. Each Ising machine has a specific graph structure, for example, the chimera structure for the current version of D-Wave and the king's graph for CMOS annealing machine. Then, let \(G_{\text{P}}=(V_{\text{P}}, E_{\text{P}})\) be the undirected graph of which the graph structure of Ising machine we want to use, which is referred to as the physical graph. Here \(V_{\text{P}}\) and \(E_{\text{P}}\) are the set of vertices and edges on the physical graph \(G_{\text{P}}\), respectively. To utilize Ising machines, we have to efficiently perform the transformation from \(G_{\text{L}}\) to \(G_{\text{P}}\).29,3133)

4. Online Advertisement Allocation Optimization from the Time-dependent Click-through Rate Data

In this section, we show the result of our study on the optimization of online advertisement allocation using a quantum annealing machine. In general, advertisements need to be served to appropriate users both from the viewpoints of the clients and users. In the case of online advertising, an indicator of the performance of serving online advertisement is the click-through rate (CTR) which is the rate at which displayed advertisements in websites are clicked. Here, the CTR fluctuates with time. It is necessary to increase the average value of CTR and to reduce the variance of CTR. The situation is similar to the portfolio optimization in economics. Thus, we construct a model based on the Markowitz mean-variance portfolio theory which is a prototypical problem in economics.

Before showing our study on online advertisement allocation optimization, we explain the Markowitz mean-variance portfolio model (Markowitz model).34) We assume that there are n financial assets in the financial market. Let \(r_{i}\) (\(1\leq i\leq n\)) be the return rate of the i-th financial asset, where \(r_{i}\) is real-valued. The Markowitz model focuses on only the average of and the variance of the expected rate of return. Let \(a_{i}\) (\(1\leq i\leq n\)) be the ratio of investment of the i-th asset. Here \(a_{i}\) is non-negative real-valued and the sum of \(a_{i}\) should be unity. The expected rate of return is defined by \begin{equation} E(\mathbf{a}) = \sum_{i=1}^{n} a_{i} r_{i}, \end{equation} (30) where \(\mathbf{a}\) is an n-dimensional vector, \(\mathbf{a}= (a_{1}, a_{2},\ldots, a_{n})\). Let \(c_{ij}\) be the covariance of assets i and j. For \(i=j\), \(v_{i}:=c_{ii}\) is the variance of return rate of the i-th asset. By using \(c_{ij}\), the variance of return rate is given by \begin{equation} V(\mathbf{a}) = \sum_{i=1}^{n} \sum_{j=1}^{n} a_{i} a_{j} c_{ij}. \end{equation} (31) The Markowitz model can consider the case that \(E(\mathbf{a})\) is maximized and \(V(\mathbf{a})\) is minimized, i.e., \begin{equation} \mathbf{a}^{*} = \mathop{\mathrm{arg\,max}}\nolimits_{\mathbf{a}} \left[E(\mathbf{a}) - \frac{\gamma}{2} V(\mathbf{a})\right],\quad \sum_{i=1}^{n} a_{i} = 1, \end{equation} (32) where γ is a hyperparameter that specifies how much risk an investor can tolerate.

The above equation is the most basic form of the Markowitz model. Reflecting that the minimum unit of investment exists, however, a discretized version of the Markowitz model, where the non-negative real variable \(a_{i}\) (\(1\leq i\leq n\)) is changed to a non-negative integer variable \(k_{i}\) (\(1\leq i\leq n\)), is often used. Let K be the total asset the investor has. The discretized version of the Markowitz model is given by \begin{align} &\mathbf{k}^{*} = \mathop{\mathrm{arg\,max}}\nolimits_{\mathbf{k}} \left[E(\mathbf{k}) - \frac{\gamma}{2} V(\mathbf{k})\right], \end{align} (33) \begin{align} &\sum_{i=1}^{n} k_{i} = K, \end{align} (34) where \(\mathbf{k}\) is an n-dimensional vector, \(\mathbf{k}= (k_{1}, k_{2},\ldots, k_{n})\).

To use Ising machines, we have to represent the non-negative integer variable \(k_{i}\) using binary variables, that is, \begin{equation} k_{i} = \sum_{d=1}^{D} f(d) y_{d}^{(i)}, \end{equation} (35) where \(y_{d}^{(i)}\) (\(1\leq d\leq D\)) are binary variables and D is the number of binary variables. Here D depends on the type of choice of \(f(d)\). For example, \(D=\lceil\log_{2} (k_{i})\rceil\) for \(f(d)=2^{d-1}\) called binary method and \(D=k_{i}\) for \(f(d)=1\) called unary method explained in Sect. 3.1.3. Additionally, the equality constraint given by Eq. (34) is represented as \begin{equation} P(\mathbf{k}) = \left(\sum_{i=1}^{n} k_{i} - K \right)^{2}, \end{equation} (36) where we use the technique explained in Sect. 3.1.2. Thus, the QUBO formulation of the Markowitz model is given by \begin{equation} \mathbf{y}^{*} = \mathop{\mathrm{arg\,max}}\nolimits_{\mathbf{y}} \left[E(\mathbf{y}) - \frac{\gamma}{2} V(\mathbf{y}) - \lambda P(\mathbf{k})\right],\quad \lambda \geq 0, \end{equation} (37) where λ is the penalty coefficient. Rosenberg et al. studied the model based on the above equation including the trade cost by using D-Wave Two.35)

Based on the previous study,35) we constructed a model to optimize online advertisement allocation. The optimization of online advertisement allocation aims to find the matching of advertisements displayed to users so as to maximize the CTR and minimize the variance of the CTR. A smaller variance of the CTR is preferred because the future profit can be more predictable when the fluctuation of the profit keeps low. The optimization we consider is regarded as the matching problem on a bipartite graph (see Fig. 1).


Figure 1. Schematic picture of our model of online advertisement allocation optimization given by Eq. (38).

Here the number of advertisements is m. Let \(r_{ij}\) be the expected CTR when the advertisement j (\(1\leq j\leq m\)) is displayed to the user i (\(1\leq i\leq n\)). We introduce the covariance of the CTR of the advertisement j to the user i and that of the advertisement \(j'\) to the user \(i'\), \(c_{ij,i'j'}\). We use the binary variable \(x_{ij}\) to represent whether the advertisement j is displayed to the user i, i.e., when advertisement j is displayed to user i, \(x_{ij}=1\), otherwise \(x_{ij}=0\). In addition, we impose two types of constraints. First, the number of advertisement displayed to the i-th user, \(k_{i}\), is specified. Second, the numbers of all advertisement displays are the same. Thus, the QUBO formulation of online advertisement allocation optimization is given by \begin{equation} \begin{split} &\mathbf{x}^{*} = \mathop{\mathrm{arg\,max}}\nolimits_{\mathbf{x}}F(\mathbf{x})\\ &F(\mathbf{x}) = E(\mathbf{x}) - \alpha V(\mathbf{x}) - \lambda_{1} \sum_{i=1}^{n} \left(\sum_{j=1}^{m} x_{ij} - k_{i} \right)^{2}\\ &\quad - \lambda_{2}\sum_{j=1}^{m} \left(\sum_{i=1}^{n} x_{ij} - \frac{\sum_{\ell=1}^{n} k_{\ell}}{m} \right),\\ &E(\mathbf{x}) = \sum_{i=1}^{n} \sum_{j=1}^{m} x_{ij} \bar{r}_{ij},\\ &V(\mathbf{x}) = \sum_{i=1}^{n} \sum_{j=1}^{m} \sum_{i'=1}^{n} \sum_{j'=1}^{m} x_{ij}x_{i'j'} c_{ij,i'j'}, \end{split} \end{equation} (38) where \(\bar{r}_{ij}\) is the normalized expected CTR given by \(r_{ij} -\frac{1}{mn}\sum_{k=1}^{n}\sum_{l=1}^{m}r_{kl}\), α is a parameter describing the coefficient of the term of variance of CTR, and \(\lambda_{1}\) and \(\lambda_{2}\) are penalty coefficients for two types of constraints. Note that we minimize \(-F(\mathbf{x})\) when we use Ising machines for the optimization.

The first term in Eq. (38) is the sum of the normalized expected CTR. The formulation with the normalized expected CTR \(\bar{r}_{ij}\) is equivalent to the one with \(r_{ij}\) when the total number of advertisement displays is fixed. We use the normalized expected CTR \(\bar{r}_{ij}\) instead of \(r_{ij}\), since the absolute value of \(\bar{r}_{ij}\) is smaller than that of \(r_{ij}\). A smaller absolute value of each term is preferred when we use Ising machines because it does not require higher resolution of the elements of QUBO matrix implemented on Ising machines. The second term in Eq. (38) is the expected variance of the profit rate. The third term in Eq. (38) is the penalty function for the constraint that the number of displays to user i is \(k_{i}\). The fourth term in Eq. (38) is also the penalty function for the constraint that the number of displays of each advertisement is \(\sum_{\ell=1}^{n} k_{\ell}/m\), i.e., each advertisement is displayed the same number of times. The first and second terms are the same as those of the discretized version of Markowitz model. Here, we have two additional constraint terms for the matching of advertisements. The meaning of the third term is similar to the third term of Eq. (37). Note that the “user” in this formulation can be replaced with a cluster of users having a similar liking.

We investigated our model with artificial data reflecting practical situation. We denote the CTR which is the advertisement j to the user i at time step t as \(\mathrm{CTR}(i,j,t)\), where t is discretized. We calculate the expected CTR in Eq. (38): \(r_{ij}\) and the covariance matrix in Eq. (38): \(c_{ij,i'j'}\) using the following equations. \begin{align} r_{ij}(t_{\text{i}},t_{\text{s}}) &= \frac{1}{t_{\text{s}}}\sum_{t=t_{\text{i}}}^{t_{\text{i}}+t_{\text{s}}-1}\mathrm{CTR}(i,j,t) \end{align} (39) \begin{align} c_{ij,i'j'}(t_{\text{i}},t_{\text{s}}) &= \frac{1}{t_{\text{s}}}\sum_{t=t_{\text{i}}}^{t_{\text{i}}+t_{\text{s}}-1} (\mathrm{CTR}(i,j,t)-r_{ij}) (\mathrm{CTR}(i',j',t)-r_{i'j'}), \end{align} (40) where \(t_{\text{i}}\) is the initial time step and \(t_{\text{s}}\) is the time steps for the optimization. We determine the matching of the advertisements at time step \(t=t_{\text{i}}+t_{\text{s}}\) by optimizing Eq. (38). The solution of this optimization is denoted as \(\mathbf{x}^{*}(t_{\text{i}},t_{\text{s}})\). Average CTR at time step \(t=t_{\text{i}}+t_{\text{s}}\) obtained by the solution \(\mathbf{x}^{*}(t_{\text{i}},t_{\text{s}})\) is called “Test CTR”. Test CTR at time step \(t=t_{\text{i}}+t_{\text{s}}\) is defined as \begin{align} &\text{Test CTR}(t_{\text{i}}+t_{\text{s}}) \notag\\ &\quad\equiv \frac{1}{\sum_{i=1}^{n}k_{i}}\sum_{i=1}^{n}\sum_{j=1}^{m}x_{ij}^{*}(t_{\text{i}},t_{\text{s}})\mathrm{CTR}(i,j,t_{\text{i}}+t_{\text{s}}). \end{align} (41) From the above equation, we calculate Test CTR values at time steps from \(t=t_{\text{s}}+1\) to \(t=T\). By using them, average and standard deviation of Test CTR for all time steps are defined by \begin{align} \text{Avg. Test CTR}&\equiv \frac{1}{T-t_{\text{s}}}\sum_{t=t_{\text{s}}+1}^{T}\text{Test CTR($t$)}, \end{align} (42) \begin{align} \text{STDEV. Test CTR}&\equiv \Biggl[\frac{1}{T-t_{\text{s}}}\sum_{t=t_{\text{s}}+1}^{T}(\text{Test CTR($t$)}\notag\\ &\quad - \text{Avg. Test CTR})^{2} \Biggr]^{\frac{1}{2}}. \end{align} (43) Avg. Test CTR represents the profit obtained by the matching optimization, and STDEV. Test CTR represents the fluctuation of the profit.

Since the definition of CTR is a ratio of the number of clicks to the number of advertisement displays, CTR can be modeled using the binomial distribution \(B(n, p)\), where n is the number of advertisement displays, and p is the CTR. When n is large enough, the binomial distribution is approximated by the normal distribution. The normal distribution can approximately model CTR since the typical number of daily advertisement displays of an advertisement is more than a thousand. In this experiment, we assume that a CTR obeys the normal distribution, i.e., the CTR values \(\mathrm{CTR}(i,j,t)\) is sampled from the multivariate normal distribution \(N(\mu,\Sigma)\), where μ and Σ are the \(nm\)-dimensional mean-vector and the \(nm\times nm\)-dimensional covariance matrix. Since μ is a real-valued vector whose each element is in \([0,1]\), we generate μ from the Beta distribution \(\text{Beta}(a,b)\), where a and b are positive real-valued parameters. To fit the actual situation where the typical value of CTR is 0.5%, we use the parameters \(a=100\) and \(b=20000\) so that the ratio is \(a/b=5\times 10^{-3}\). The standard derivation of μ is about 10% of the average value of μ, which is valid to consider the CTR of the actual advertisement. In addition, we generate the covariance matrix Σ from the inverse Wishart distribution \(IW(\Psi, s)\) which is used for the prior distribution of covariance of multivariate normal distribution, where s is the degree of freedom parameter and Ψ is the scale matrix which is an \(nm\times nm\) matrix. Here we use the parameters, \(s=nm+1\) and \(\Psi=1\times 10^{-6}\mathbb{E}\), where \(\mathbb{E}\) is the \(nm\times nm\) identity matrix. The square roots of diagonal elements of the covariance matrix Σ represent the fluctuation of CTR per unit time width. The square root of diagonal elements is about 30% of the CTR value, which is valid compared to the actual CTR fluctuation.

We compare the performance of D-Wave 2000Q and that of the greedy method with the artificially generated data as mentioned above. In our experiment, the following parameters are used. The number of advertisements is \(n=8\), the number of users is \(m=8\), the number of advertisement displays is \(k_{i}=1\) for all users i, the total time steps is set to \(T=40\), and the number of the time steps for optimization is set to \(t_{\text{s}} = 30\). The value of the penalty parameters \(\lambda_{1},\lambda_{2}\) are set to a larger value than the value of the normalized profit rates and the covariance matrix to obtain feasible solutions with a high probability. Specifically, we set \(\lambda_{1}=\lambda_{2}=0.6A\), where \(A=\max_{k,l}Q_{k,l}^{\text{CTR}}\) and \(\mathbf{x}^{\top}Q^{\text{CTR}}\mathbf{x}=E(\mathbf{x})-\alpha V(\mathbf{x})\). The value of A is the largest absolute value of the QUBO matrix of the profit rate and the covariance matrix. In the greedy method, we initialize the state so that \(x_{i,j}=0\) for all i and j. Next, we select the edge \(x_{i,j}\) which maximizes \(\mathbf{x}^{\top}Q^{\text{CTR}}\mathbf{x}\), and set \(x_{i,j}=1\). We repeat this process until there is no edges to select. To obtain feasible solutions, we need to select the edge such that the number of connected edges to the user i is less than \(k_{i}\) and the number of connected edges to the advertisement j is less than \(k'_{j}\), where \(k'_{j}\equiv (\sum_{i=1}^{n}k_{i})/m\).

Now, let us look at the optimization with a quantum annealing machine, D-Wave 2000Q. The quantum annealing machine used in our experiment does not have a complete chimera graph because of the qubit yield. Thus, we employed the virtual hardware graph solver named VFYC (Virtual Full Yield Chimera)36) which enables us to solve the problem on the abstract complete chimera graph by applying post-processing to the solution obtained by the quantum annealing machine with missing qubits. We solve the QUBO of our problem with 64 fully connected logical spins realized by an embedding method for complete graphs.37) The coupling strength in the embedding is set to \(2A\). In this experiment, we fix the annealing time at 20 µs, the number of annealing runs at 1000 and the number of spin reversal transformations38) at 4. We also used the optimization postprocessing provided by the D-Wave 2000Q.39)

When we obtain the solution with the optimization mentioned above method, we calculate the Test CTR at each time step by Eq. (41). Furthermore, the average and the standard deviation of the Test CTR for all time steps are calculated using Eqs. (42) and (43). We also sampled the parameters μ and Σ of the multivariate normal distribution from the Beta and inverse Wishart distribution respectively as explained before. The expected value of Avg. Test CTR, STDEV. Test CTR and the objective function of Eq. (38) were calculated using 5 sampled parameters.

We compared the dependencies of Avg. Test CTR and STDEV. Test CTR on the coefficient of the variance α with different optimization methods (Fig. 2). In Fig. 2(a), we observed that the Avg. Test CTR obtained by D-Wave 2000Q is larger than that obtained by the greedy method at \(\alpha\geq 50\) while no significant difference is observed at \(\alpha=0\). In Fig. 2(b), the STDEV. Test CTR decreases as α decreases in both methods, which is consistent with the fact that α is the coefficient which suppresses the fluctuation of CTR. Moreover, the standard deviation obtained by the D-Wave 2000Q is smaller than that obtained by the greedy method at \(\alpha\geq 50\).


Figure 2. Dependencies of (a) the average, (b) the standard deviation and (c) the objective function of Eq. (38) on the variance coefficient α, using the greedy method and the quantum annealing performed using D-Wave 2000Q.

Next, let us look at the objective function \(F(\mathbf{x})\) at \(\alpha\geq 50\) where the Avg. Test CTR obtained by D-Wave 2000Q is larger than that obtained by the greedy method as shown in Fig. 2(c). Note that the optimization problem is to find the condition so as to maximize the objective function \(F(\mathbf{x})\). The objective function of D-Wave 2000Q is larger than that of the greedy method as shown in Fig. 2(c), which means that the solution precision obtained by D-Wave 2000Q is higher than that obtained by the greedy method. The result suggests that the higher Avg. Test CTR and the smaller STDEV. Test CTR obtained by D-Wave 2000Q is owing to the better optimization.

From the above results, we conclude that the profit of the online advertisement allocation and the standard deviation of the CTR is improved by the proposed method with a quantum annealing machine.

5. Software

One of the authors (K.T.) developed software called PyQUBO that helps to abstract codes.30) The advantages of using PyQUBO are readability of the code which leads to the prevention of errors in the code, reduction of compilation time, automatic validation whether the given constraints are satisfied. In this section, we provide a short explanation of PyQUBO.

To utilize Ising machines, we need to represent combinatorial optimization problems using the Ising model or the QUBO as described in Sect. 3. Here, we consider the case where the cost function and the equality/inequality constraints of the combinatorial optimization problem are formulated in terms of polynomial functions with binary variables. If products of three or more binary variables exist in the polynomial function, we need to reduce the higher-order interactions by introducing auxiliary variables, as explained in Sect. 3. After the reduction of higher-order interactions, we obtain the QUBO matrix corresponding to the combinatorial optimization problem we want to solve.

Practically, we need to write a program to prepare the QUBO matrix corresponding to the combinatorial optimization problem. To obtain the QUBO matrix, we rewrite the Hamiltonian such that the linear and quadratic terms of the spin variables or 0–1 binary variables are separate. The code of a program written to achieve such program created by the above concept tends to be complex and unreadable. Let us see an example of this process for the number partitioning problem. As described in Sect. 3.1.1, the number partitioning problem aims to partition a set of given numbers into two groups such that the sums of the numbers in each group are equal. The Hamiltonian of this problem can be written by a format similar to Eq. (13).

To write a program for preparing the QUBO, we need to manually expand the Hamiltonian as \begin{equation} \mathcal{H} = 2\sum_{i<j}n_{i}n_{j}s_{i}s_{j} + \sum_{i} n_{i}^{2}. \end{equation} (44) Using the above equation and dimod developed by D-Wave Systems,40) we can write a code to produce the QUBO in Python (Table I).

Data table
Table I. Python code for building the QUBO for the number partitioning problem.

The code is unreadable since the original form of the Hamiltonian, i.e., before the expansion of the Hamiltonian, cannot be directly seen in the code. No clear correspondence between the Hamiltonian and its expression in the code exists. Thus, it may be difficult to fix the code if one wants to modify the Hamiltonian partially. Additionally, this program is susceptible to software bugs from miscalculations because the expansion of the Hamiltonian needs to be manually calculated when there are many equality/inequality constraints and the form of the cost function is complicated. The serious situation often occurs when we consider a complicated combinatorial optimization problem in the real world. To consider the complicated combinatorial optimization problem, we start with a simple Hamiltonian to obtain the first stage of a proof-of-principle result. Next, we modify the Hamiltonian to be acceptable for real-world applications based on the proof-of-principle result. Since the above strategy of program development has been adopted in many industrial products, we should write readable and extendable code. With PyQUBO, we can define a Hamiltonian with spin or binary objects in a more straightforward way. We demonstrate construction of the QUBO matrix for the number partitioning problem with PyQUBO (Table II).

Data table
Table II. PyQUBO code to produce a QUBO of number partitioning problem.

Firstly, we define the array of spin variables (line 6) and the Hamiltonian of the problem (line 9). Then, we compile the Hamiltonian to get the model (line 12). Finally, the QUBO is generated by calling the to_qubo() method of the model (line 15). The code becomes much straightforward and readable than the conventional one in Table I, which leads to the extendable code.

The Hamiltonian in PyQUBO is internally represented as a tree. For example, the Hamiltonian for the number partitioning problem is represented as a tree shown in the left panel of Fig. 3. The compilation process of PyQUBO is shown in the right panel of Fig. 3. Firstly, the Hamiltonian is converted to a higher-order polynomial by folding the tree. Then, a higher-order polynomial is reduced to a second-order polynomial by introducing auxiliary variables. Finally, the QUBO matrix is created from the coefficients of the second-order polynomial. In the next section, some features of PyQUBO are introduced, including placeholder and automatic validation of constraints.


Figure 3. Tree representation of the Hamiltonian for the number partitioning problem in PyQUBO (left). Compilation process of PyQUBO (right).

Placeholder

To obtain better feasible solutions with high probability, we should update penalty coefficient parameters in the Hamiltonian. If we compile the QUBO every time we change the penalty coefficient, the large number of compilations is needed, and as a result, we spend much compilation time. However, using a placeholder can reduce the number of compilations. The actual value of the parameters defined by the placeholder can be specified after we compile the QUBO, which means that we can obtain QUBOs with various parameter values with only one compilation. Let us see an example of using the placeholder for the graph partitioning problem (Table III). For this case, we need to tune the coefficient λ in the Hamiltonian given by Eq. (14). We define the penalty coefficient λ by the placeholder (line 14) so that we can change the value even after compilation. The actual value of the penalty strength is specified when the to_qubo() method is called (line 27).

Data table
Table III. Example of using Placeholder and Constraint in PyQUBO for the graph partitioning problem.
Automatic validation of constraints

When we study combinatorial optimization problems with equality/inequality constraints using Ising machines, we should prepare additional functions in the code to judge whether the obtained results satisfy the given constraints. However, we do not need to write additional functions in PyQUBO. We can tell the compiler which terms in the Hamiltonian represent equality/inequality constraints using the Constraint class in PyQUBO. By just using Constraint, we can extract the information about unsatisfied constraints in the obtained results. In the example for the graph partitioning problem (Table III), the Hamiltonian for the problem is contained in a Constraint object (line 17). The information about the unsatisfied constraint is obtained when you decode the solution (line 29). If the broken object is empty, it means that all constraints are satisfied.

6. Conclusion

In the first half of this paper, we reviewed combinatorial optimization problems and how to express the combinatorial optimization problem using the Ising model or the QUBO formulation. In Sect. 4, we showed our experimental result for online advertisement allocation optimization using D-Wave 2000Q. We calculated the CTR that measures the benefit of online advertisement. We observed that the average value of CTR obtained using the greedy method is larger than that obtained using D-Wave 2000Q, and the variance of CTR obtained using D-Wave 2000Q is less than that obtained using the greedy method. Thus, we confirmed that quantum annealing outperformed the greedy method for our problem. In Sect. 5, a software developed, called PyQUBO, which conveniently helps to utilize Ising machines was explained.

Acknowledgment

One of the authors (S.T.) was supported by JST, PRESTO Grant Number JPMJPR1665, Japan and JSPS KAKENHI Grant Numbers 15K17720, 15H03699.


References

  • 1 M. W. Johnson, M. H. S. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, E. M. Chapple, C. Enderud, J. P. Hilton, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva, C. J. S. Truncik, S. Uchaikin, J. Wang, B. Wilson, and G. Rose, Nature 473, 194 (2011). 10.1038/nature10012 CrossrefGoogle Scholar
  • 2 M. Yamaoka, C. Yoshimura, M. Hayashi, T. Okuyama, H. Aoki, and H. Mizuno, IEEE J. Solid-State Circuits 51, 303 (2016). 10.1109/JSSC.2015.2498601 CrossrefGoogle Scholar
  • 3 T. Inagaki, Y. Haribara, K. Igarashi, T. Sonobe, S. Tamate, T. Honjo, A. Marandi, P. L. McMahon, T. Umeki, K. Enbutsu, O. Tadanaga, H. Takenouchi, K. Aihara, K. Kawarabayashi, K. Inoue, S. Utsunomiya, and H. Takesue, Science 354, 603 (2016). 10.1126/science.aah4243 CrossrefGoogle Scholar
  • 4 P. L. McMahon, A. Marandi, Y. Haribara, R. Hamerly, C. Langrock, S. Tamate, T. Inagaki, H. Takesue, S. Utsunomiya, K. Aihara, R. L. Byer, M. M. Fejer, H. Mabuchi, and Y. Yamamoto, Science 354, 614 (2016). 10.1126/science.aah5178 CrossrefGoogle Scholar
  • 5 S. Matsubara, H. Tamura, M. Takasu, D. Yoo, B. Vatankhahghadim, H. Yamasaki, T. Miyazawa, S. Tsukamoto, Y. Watanabe, K. Takemoto, and A. Sheikholeslami, Conference on Complex, Intelligent, and Software Intensive Systems, 2017, 432. 10.1007/978-3-319-61566-0_39 CrossrefGoogle Scholar
  • 6 Web [https://annealing-cloud.com]. Google Scholar
  • 7 A. Lucas, Front. Phys. 2, 5 (2014). 10.3389/fphy.2014.00005 CrossrefGoogle Scholar
  • 8 S. Tanaka, R. Tamura, and B. K. Chakrabarti, Quantum Spin Glasses, Annealing and Computation (Cambridge University Press, Cambridge, U.K., 2017). Google Scholar
  • 9 V. S. Denchev, S. Boixo, S. V. Isakov, N. Ding, R. Babbush, V. Smelyanskiy, J. Martinis, and H. Neven, Phys. Rev. X 6, 031015 (2016). 10.1103/PhysRevX.6.031015 CrossrefGoogle Scholar
  • 10 F. Neukart, G. Compostella, C. Seidel, D. von Dollen, S. Yarkoni, and B. Parney, Frontiers in ICT 4, 29 (2017). 10.3389/fict.2017.00029 CrossrefGoogle Scholar
  • 11 K. Terada, D. Oku, S. Kanamaru, S. Tanaka, M. Hayashi, M. Yamaoka, M. Yanagisawa, and N. Togawa, Proc. Int. Symp. VLSI Design, Automation and Test (VLSI-DAT), 2018, 2018. 10.1109/VLSI-DAT.2018.8373233 CrossrefGoogle Scholar
  • 12 K. Kitai, J. Guo, S. Ju, S. Tanaka, K. Tsuda, J. Shiomi, and R. Tamura, arXiv:1902.06573. Google Scholar
  • 13 T. Kadowaki and H. Nishimori, Phys. Rev. E 58, 5355 (1998). 10.1103/PhysRevE.58.5355 CrossrefGoogle Scholar
  • 14 S. Geman and D. Geman, IEEE Trans. Pattern Anal. Mach. Intell. PAMI-6, 721 (1984). 10.1109/TPAMI.1984.4767596 CrossrefGoogle Scholar
  • 15 M. Suzuki, Prog. Theor. Phys. 56, 1454 (1976). 10.1143/PTP.56.1454 CrossrefGoogle Scholar
  • 16 H. F. Trotter, Proc. Am. Math. Soc. 10, 545 (1959). 10.1090/S0002-9939-1959-0108732-6 CrossrefGoogle Scholar
  • 17 T. Kadowaki, Ph. D. thesis, Tokyo Institute of Technology (1999). Google Scholar
  • 18 K. Tanaka and T. Horiguchi, Electron. Commun. Jpn. 83 [3], 84 (2000). 10.1002/(SICI)1520-6440(200003)83:3<84::AID-ECJC9>3.0.CO%3B2-N CrossrefGoogle Scholar
  • 19 R. Martoňák, G. E. Santoro, and E. Tosatti, Phys. Rev. E 70, 057701 (2004). 10.1103/PhysRevE.70.057701 CrossrefGoogle Scholar
  • 20 D. A. Battaglia, G. E. Santoro, and E. Tosatti, Phys. Rev. E 71, 066707 (2005). 10.1103/PhysRevE.71.066707 CrossrefGoogle Scholar
  • 21 K. Kurihara, S. Tanaka, and S. Miyashita, Proc. 25th Conf. Uncertainty in Artificial Intelligence (UAI2009), 2009. Google Scholar
  • 22 I. Sato, K. Kurihara, S. Tanaka, H. Nakagawa, and S. Miyashita, Proc. 25th Conf. Uncertainty in Artificial Intelligence (UAI2009), 2009. Google Scholar
  • 23 I. Sato, S. Tanaka, K. Kurihara, S. Miyashita, and H. Nakagawa, Neurocomputing 121, 523 (2013). 10.1016/j.neucom.2013.05.019 CrossrefGoogle Scholar
  • 24 M. Ohzeki, S. Okada, M. Terabe, and S. Taguchi, Sci. Rep. 8, 9950 (2018). 10.1038/s41598-018-28212-4 CrossrefGoogle Scholar
  • 25 M. Maezawa, K. Imafuku, M. Hidaka, H. Koike, and S. Kawabata, IEEE Conf. Proc. 16th Int. Superconductive Electronics Conf. (ISEC2017), 2017. Google Scholar
  • 26 M. Maezawa, G. Fujii, M. Hidaka, K. Imafuku, K. Kikuchi, H. Koike, K. Makise, S. Nagasawa, H. Nakagawa, M. Ukibe, and S. Kawabata, J. Phys. Soc. Jpn. 88, 061012 (2019). 10.7566/JPSJ.88.061012 LinkGoogle Scholar
  • 27 D. O’Malley and V. V. Vesselinov, IEEE High Performance Extreme Computing Conference (HPEC), 2016, 2016. 10.1109/HPEC.2016.7761616 CrossrefGoogle Scholar
  • 28 Web [https://docs.ocean.dwavesys.com/]. Google Scholar
  • 29 T. D. Goodrich, B. D. Sullivan, and T. S. Humble, Quantum Inf. Process. 17, 118 (2018). 10.1007/s11128-018-1863-4 CrossrefGoogle Scholar
  • 30 Web [https://github.com/recruit-communications/pyqubo]. Google Scholar
  • 31 J. Cai, B. Macready, and A. Roy, arXiv:1406.2741. Google Scholar
  • 32 D. Oku, K. Terada, M. Hayashi, M. Yamaoka, S. Tanaka, and N. Togawa, submitted. Google Scholar
  • 33 N. Dattani and N. Chancellor, arXiv:1901.07676. Google Scholar
  • 34 H. Markowitz, J. Finance 7, 77 (1952). 10.1111/j.1540-6261.1952.tb01525.x CrossrefGoogle Scholar
  • 35 G. Rosenberg, P. Haghnegahdar, P. Goddard, P. Carr, K. Wu, and M. L. de Prado, IEEE J. Sel. Top. Signal Process. 10, 1053 (2016). 10.1109/JSTSP.2016.2574703 CrossrefGoogle Scholar
  • 36 Web [https://docs.dwavesys.com/docs/latest/c_post-processing_5.html]. Google Scholar
  • 37 V. Choi, Quantum Inf. Process. 10, 343 (2011). 10.1007/s11128-010-0200-3 CrossrefGoogle Scholar
  • 38 S. Boixo, T. F. Rønnow, S. V. Isakov, Z. Wang, D. Wecker, D. A. Lidar, J. M. Martinis, and M. Troyer, Nat. Phys. 10, 218 (2014). 10.1038/nphys2900 CrossrefGoogle Scholar
  • 39 Web [https://docs.dwavesys.com/docs/latest/c_post-processing_1.html]. Google Scholar
  • 40 Web [https://github.com/dwavesystems/dimod]. Google Scholar

Author Biographies


Kotaro Tanahashi received the M.Eng. from Kyoto University in 2015. He currently works for Recruit Communications Co., Ltd. as a machine learning engineer. He is also a project manager of MITOU Target Program at Information-technology Promotion Agency (IPA).

Shinichi Takayanagi received the M.Sc. from Hokkaido University in 2006. He is presently an data scientist at LINE Corporation. His research interests are quantum annealing, statistical mechanics, and machine learning. He is a member of the Physical Society of Japan and Information Processing Society of Japan.

Tomomitsu Motohashi received the M.Eng. from Sophia University in 2009. He is presently a CTO at digital medical start-up Susmed Corporation. His research interests are quantum annealing, blockchain, operations research, and machine learning. He is a member of Information Processing Society of Japan. He was awarded the encouraging award in 2009 from The Institute of Systems, Control and Information Engineers of Japan and 2nd prise of KDD Cup 2015.

Shu Tanaka received the B.Sc. degree from Tokyo Institute of Technology in 2003 and the M.Sc. and Dr.Sc. degrees from The University of Tokyo in 2005 and 2008, respectively. He is presently an associate professor at Green Computing Systems Research Organization, Waseda University, a PRESTO researcher at Japan Science and Technology Agency (JST), a project manager of MITOU Target Program at Information-technology Promotion Agency (IPA). His research interests are quantum annealing, Ising machine, statistical mechanics, and materials science. He was awarded the 9th Award for the Encouragement of Young Physicists by the Physical Society of Japan in 2015. He is a member of Physical Society of Japan.

Cited by

View all 32 citing articles

no accessMany-Qudit Representation for the Travelling Salesman Problem Optimisation

114002, 10.7566/JPSJ.90.114002

no accessPulsed Quantum Annealing

094003, 10.7566/JPSJ.89.094003

no accessDerivation of QUBO Formulations for Sparse Estimation

034801, 10.7566/JPSJ.89.034801

free accessDeeper Understanding of Constrained Quantum Annealing from the Perspective of the Localization Phenomena

10, 10.7566/JPSJNC.17.10