Efficient correlation-based discretization of continuous variables for annealing machines

The Ising model is one of the most fundamental models to describe the properties of magnetic materials. Let $\sigma_{i}\in\{-1,1\}$ denote a variable for the $i$-th spin. Then, the energy function of the Ising model with the state vector $\bm{\sigma}$ is defined as follows:

$\displaystyle E(\bm{\sigma})=-\sum_{i<j}J_{ij}\sigma_{i}\sigma_{j}-\sum_{i}h_{i}\sigma_{i},$

(1)

where $J_{ij}\in\mathbb{R}$ corresponds to the two-body interaction between the spins $\sigma_{i}$ and $\sigma_{j}$, and $h_{i}\in\mathbb{R}$ is the external magnetic field on $\sigma_{i}$.

The QUBO formulation has the following cost function for the state vector $\bm{z}$:

$\displaystyle E(\bm{z})=-\sum_{i,j}Q_{ij}z_{i}z_{j},$

(2)

where $z_{i}\in\{0,1\}$ is $i$-th binary variable in $\bm{z}$, and $Q_{ij}\in\mathbb{R}$ is the strength of the interaction between the binary variables $z_{i}$ and $z_{j}$. This formulation is equivalent to the Ising model via the variable transformation with

$\displaystyle z_{i}=\frac{1+\sigma_{i}}{2}.$

(3)

Conversion between $\{J_{ij}\},\{h_{i}\}$ and $\{Q_{ij}\}$ is also possible.

From the computational viewpoint, annealing machines are domain-specific, and their role is simple; to find the ground state that minimizes Eqs. (1) or (2). Despite the domain-specific characteristic, we can solve various combinatorial optimization problems with annealing machines. The reason is that various optimization problems are converted into the QUBO formulation [4]; we formulate the Ising model or QUBO formulation so that the ground state coincides with the optimal solution of the combinatorial optimization problem. Hence, we can solve combinatorial optimization problems with the annealing machines.

As stated above, the Ising model and the QUBO formulation are equivalent. However, the QUBO formulation is more suitable for computation, and then we employ the QUBO formulation below.

The QUBO formulation is widely applicable to different types of annealing machines.

There are two main types of annealing machines currently under development. One uses simulated annealing (SA) and the other uses quantum annealing (QA). Although the mechanisms of the search process differ, both types accept the QUBO formulation as the input for specific hardware. We briefly review why the QUBO formulation is suitable for both types of machines.

SA is an algorithm based on thermal fluctuations. When the temperature is high, the thermal fluctuations are large. The algorithm searches the large area in the state spaces. A gradual decrease in temperature indicates the settling down to the ground state. Numerically, the probability $P$ of a state change is often defined as follows:

$\displaystyle P=\min\left[1,\exp\left(-\frac{\Delta E}{T}\right)\right],$

(4)

where $\Delta E$ is the energy difference when a certain state is changed, and $T$ is temperature. The ground state is adequately obtained if the temperature decreases slowly enough [17]. Note that $\Delta E$ in Eq. (4) is evaluated using Eq. (2). The algorithm uses the cost function of the QUBO formulation. It indicates that the QUBO formulation is suitable as the input for machines.

QA was proposed by Kadowaki and Nishimori and uses quantum fluctuation [18]. Instead of the classical binary spin variables, we need quantum spins denoted by the Pauli matrix. Introducing the state vector with the $z$ component of the Pauli matrices, $\hat{\bm{\sigma}}^{z}$, and the $x$ component for $i$-th spin, $\hat{\sigma}_{i}^{x}$, the Hamiltonian is defined as follows:

$\displaystyle\widehat{H}=E(\hat{\bm{\sigma}}^{z})-\Gamma\sum_{i}\hat{\sigma}^{x}_{i},$

(5)

where $\Gamma$ is a control parameter. The second term in the right-hand side (r.h.s.) of Eq. (5) corresponds to the quantum fluctuations. The initial state is prepared as the ground state under quantum fluctuations. By gradually weakening the quantum effect, i.e., $\Gamma$, QA seeks the ground state determined by the first term on the r.h.s. of Eq. (5). QA can use more rapid scheduling than SA [19]. Such a computation method is also called adiabatic quadratic computing [20].

The annealing mechanism of QA differs from that of SA. However, the final ground state is defined only by the first term on the r.h.s. of Eq. (5). It is straightforward to obtain the first term from the QUBO formulation in Eq. (2). Hence, practical QA machines, such as the D-Wave Advantage, accept the QUBO formulation as the input.

The main topic of this paper is the discretization of continuous variables. Although the annealing mechanisms vary between SA and QA, we can use the following discussions with the QUBO formulation for both SA and QA machines. Because of the simpler notations, only the classical cases, i.e., SA cases, are considered in the following sections; we do not consider any quantum effects.

Linear regression is one of the fundamental methods in data analysis with which we seek the relationship between a response variable and explanatory variables. Although matrix operations are sufficient to solve simple linear regression problems, further computational efforts are necessary in large-scale cases or when there are regularization terms. Hence, the use of an annealing machine is desirable. With these problems in mind, Date and Potok discussed the implementation and increasing the speed of simple linear regression problems using an annealing machine [15]. Linear regression problems also yield an example to investigate discretization methods of continuous variables. We briefly review Date and Potok’s discussion.

A linear regression model is defined as

$\displaystyle y=w_{1}+w_{2}x_{2}+\cdots+w_{D}x_{D},$

(6)

where $w_{d}\in\mathbb{R}\,\,(d=1,\dots,D)$ is the $d$-th parameter, and $x_{d}\in\mathbb{R}\,\,(d=2,\dots,D)$ is the $d$-th input variable. They are summarized in the following forms:

	$\displaystyle\bm{w}$	$\displaystyle=[w_{1},w_{2},\cdots,w_{D}]^{\mathrm{T}}\in\mathbb{R}^{D},$		(7)
	$\displaystyle\bm{x}$	$\displaystyle=[1,x_{2},\cdots,x_{D}]^{\mathrm{T}}\in\mathbb{R}^{D},$		(8)

where the first component in $\bm{x}$ is a dummy variable; the dummy variable is convenient for the vector expression of Eq. (6) as follows:

$\displaystyle y=\bm{w}^{\mathrm{T}}\bm{x}.$

(9)

In linear regression problems, we seek the model parameters $\bm{w}$ in Eq. (6) that are suitable for our obtained data. Let $\{(\bm{x}_{i},y_{i}):i=1,2,\cdots,N\}$ be a set of $N$ training data, where $\bm{x}_{i}=[1,x_{i2},\cdots,x_{iD}]^{\mathrm{T}}\in\mathbb{R}^{D}$ is a $D$-dimensional input vector, and $y_{i}\in\mathbb{R}$ is an output for $i$-th data. We introduce an output vector $\bm{y}=[y_{1},y_{2},\cdots,y_{N}]^{\mathrm{T}}$ and matrix $X\in\mathbb{R}^{N\times D}$ defined as

$\displaystyle X=\begin{bmatrix}\,1&x_{12}&\cdots&x_{1D}\\ \,1&x_{22}&\cdots&x_{2D}\\ \,\vdots&\vdots&\ddots&\vdots\\ \,1&x_{N2}&\cdots&x_{ND}\\ \end{bmatrix}.$

(10)

Then, a conventional linear regression problem has the following cost function:

$\displaystyle\widetilde{E}(\bm{w})=\|\bm{y}-X\bm{w}\|^{2}=\bm{y}^{\mathrm{T}}\bm{y}-2\bm{w}^{\mathrm{T}}X^{\mathrm{T}}\bm{y}+\bm{w}^{\mathrm{T}}X^{\mathrm{T}}X\bm{w}.$

(11)

In Eq. $\eqref{eq:pre_reg_cost}$, $\bm{y}^{\mathrm{T}}\bm{y}$ is meaningless in optimizing $\bm{w}$. Hence, Eq. $\eqref{eq:pre_reg_cost}$ comes down to the following problem:

$\displaystyle E(\bm{w})=\bm{w}^{\mathrm{T}}X^{\mathrm{T}}X\bm{w}-2\bm{w}^{\mathrm{T}}X^{\mathrm{T}}\bm{y}.$

(12)

Of course, it is possible to find the minimum for Eq. (12) via simple matrix operations. However, the matrix operations include a matrix inversion with a high computational cost. Annealing machines directly solve the optimization problem in the form of Eq. (12).

As denoted in Sec. I, Date and Potok showed that the QUBO formulation and the annealing machines enable faster computation than conventional classical computers even in simple linear regression problems without any regularization term [15]. We now explain the QUBO formulation for simple linear regression problems.

First, we introduce $\bm{b}=[b_{1},b_{2},...,b_{K}]^{\mathrm{T}}\,\,(K\in\mathbb{N})$ as a basis vector. For later use, the components in $\bm{b}$ are in ascending order of the absolute value. An example of this vector is $\bm{b}=[\frac{1}{2},-\frac{1}{2},1,-1,2,-2]^{\mathrm{T}}$.

Second, we introduce a $D\times K$ binary matrix $\widetilde{\bm{z}}\in\{0,1\}^{D\times K}$. Then, the discretization of the linear regression weight $w_{i}$ is

$\displaystyle w_{i}=\sum_{k=1}^{K}b_{k}\widetilde{z}_{ik}\quad(i\in 1,2,\dots,D),$

(13)

which indicates that the binary variable $\widetilde{z}_{ik}$ acts like a flag to determine whether the basis $b_{k}$ is used for the expression. For the above example of $\bm{b}=[\frac{1}{2},-\frac{1}{2},1,-1,2,-2]^{\mathrm{T}}$, we can express the values with $\{-\frac{7}{2},-3,\dots,3,\frac{7}{2}\}$. Note that there may not be a one-to-one correspondence between $w_{i}$ and $\widetilde{\bm{z}}$; we can use redundant expressions for original variables. Although the redundant basis could yield good annealing performance, this point is not our main research topic.

Third, let $\bm{z}=[\widetilde{z}_{11},\cdots,\widetilde{z}_{1K},\widetilde{z}_{21},\cdots,\widetilde{z}_{2K},\cdots,\widetilde{z}_{D1},\cdots,\widetilde{z}_{DK}]^{T}$ denote the flattened one-dimensional vector derived from the binary matrix $\widetilde{\bm{z}}$. Then, we rewrite Eq. (13) in a matrix form as follows:

$\displaystyle\bm{w}=B\bm{z},$

(14)

where $B=I_{D}\otimes\bm{b}$ is the basis matrix derived from the Kronecker product of the $D$-dimensional unit matrix $I_{D}$ and $\bm{b}$. An example of the matrix $B$ is given later.

Finally, we obtain the QUBO formulation of linear regression from Eqs. (12) and (14) as follows:

$\displaystyle E(\bm{z})=\bm{z}^{\mathrm{T}}B^{\mathrm{T}}X^{\mathrm{T}}XB\bm{z}-2\bm{z}^{\mathrm{T}}B^{\mathrm{T}}X^{\mathrm{T}}\bm{y}.$

(15)

As mentioned in Sec. III.2, the size of the binary vector $\bm{z}$ is $D\times K$ because we assign $K$ binary variables to each continuous variable with the naive discretization method. However, for large-scale problems, this method is likely to exceed the size limit of the annealing machine. How can we reduce the number of binary variables within the size limit of the annealing machine?

The simplest method is to reduce the components of the basis assigned in the binary expansion. However, this discretization reduces the precision of the original continuous variables and the ranges that the variables can represent.

One may employ the discretization method for each continuous variable. The flexible discretization may yield good performance. However, we cannot determine the required expression for each original continuous variable before solving the problem.

It is also possible to choose randomly continuous variables with reduced bases; this randomized method might yield a good accuracy on average. However, it could degrade the accuracy of a realization.

In summary, the reduction in the basis degrades the precision of each original continuous variable and yields narrow ranges. Hence, this reduction could lead to poor prediction accuracy for the linear regression.

To avoid the degradation in precision and narrower ranges, we exploit the same binary variables for discretization. For example, $\bm{b}=[\frac{1}{2},-\frac{1}{2},1,-1,2,-2]^{\mathrm{T}}$ expresses $-2.5$ and $-1.5$ as

	$\displaystyle-2.5$	$\displaystyle=\begin{bmatrix}0,1,0,0,0,1\end{bmatrix}\bm{b}$
		$\displaystyle=0\times b_{1}+1\times b_{1}+0\times b_{2}+\cdots+1\times b_{6},$
	$\displaystyle-1.5$	$\displaystyle=\begin{bmatrix}1,0,0,0,0,1\end{bmatrix}\bm{b}$
		$\displaystyle=1\times b_{1}+0\times b_{1}+0\times b_{2}+\cdots+1\times b_{6}.$

When we consider two large negative values, $b_{6}$ will commonly be used. Hence, the coefficient of $b_{6}$ is always $1$. If two variables are strongly correlated, the coefficients of $b_{6}$ for them will be the same. That is, the change in the coefficients of $b_{6}$ largely affects these two variables, which indicates the strong correlation between them. Hence, we can use the same binary variable for $b_{6}$ in this example.

Note that using the same binary variable does not degrade the precision of approximation, and the range is unchanged. These are different from the naive methods in Sec. IV.1. Note that we must find pairs of continuous variables with strong correlations. Hence, we need a preliminary evaluation of correlations in the original continuous variables. Although this evaluation takes some computational effort, only rough estimations for the correlations are sufficient; the computational cost is not high.

The basic idea of the proposed method is summarized as follows: evaluate correlations, find strongly correlated pairs, and reduce redundant binary variables. In the following sections, we first explain how to reduce the number of binary variables and then give the entire procedure to derive the QUBO formulation.

As reviewed in Sec. III.2, continuous parameters $\bm{w}$ are discretized using Eq. (14). For example, we here consider $D=2$, $K=3$, $\bm{b}=[b_{1},b_{2},b_{3}]^{\mathrm{T}}$. Then, Eq. (14) yields the following discretization:

$\displaystyle\begin{bmatrix}w_{1}\\ w_{2}\\ \end{bmatrix}=\begin{bmatrix}b_{1}&b_{2}&b_{3}&0&0&0\\ 0&0&0&b_{1}&b_{2}&b_{3}\\ \end{bmatrix}\begin{bmatrix}\widetilde{z}_{11}\\ \widetilde{z}_{12}\\ \widetilde{z}_{13}\\ \widetilde{z}_{21}\\ \widetilde{z}_{22}\\ \widetilde{z}_{23}\\ \end{bmatrix}.$

(16)

That is, $3$ bits are assigned to each of the continuous variables $w_{1}$ and $w_{2}$, and the total number of binary variables is $6$.

We assume that the continuous variables $w_{1}$ and $w_{2}$ are strongly correlated. We then modify matrix $B$ by considering the strong correlation between the continuous variables. In accordance with the modification, the size of $\bm{z}$ is also reduced. The result of the modification is as follows:

$\displaystyle\begin{bmatrix}w_{1}\\ w_{2}\\ \end{bmatrix}=\begin{bmatrix}b_{1}&b_{2}&0&0&b_{3}\\ 0&0&b_{1}&b_{2}&b_{3}\\ \end{bmatrix}\begin{bmatrix}\widetilde{z}_{11}\\ \widetilde{z}_{12}\\ \widetilde{z}_{21}\\ \widetilde{z}_{22}\\ \widetilde{z}_{23}\\ \end{bmatrix}.$

(17)

Note that the redundant binary variable $\widetilde{z}_{13}$ is deleted by sharing the binary variable $\widetilde{z}_{23}$ with $w_{1}$ and $w_{2}$. If further reduction is necessary, we can use the following expression:

$\displaystyle\begin{bmatrix}w_{1}\\ w_{2}\\ \end{bmatrix}=\begin{bmatrix}b_{1}&0&b_{2}&b_{3}\\ 0&b_{1}&b_{2}&b_{3}\\ \end{bmatrix}\begin{bmatrix}\widetilde{z}_{11}\\ \widetilde{z}_{21}\\ \widetilde{z}_{22}\\ \widetilde{z}_{23}\\ \end{bmatrix}.$

(18)

In summary, we seek a reduced matrix $B^{\prime}$ and the corresponding vector for the binary variables $\bm{z}^{\prime}$ with the following form:

$\displaystyle\bm{w}=B^{\prime}\bm{z}^{\prime}.$

(19)

Note that the binary variable with a large absolute value is shared first. As mentioned in Sec. III.2, $\bm{b}$ is in ascending order. Therefore, the order of sharing is $b_{K},b_{K-1},\dots,b_{1}$. For a linear regression problem, there is a change in the optimization variables, and Eq. (15) is modified as

$\displaystyle E(\bm{z}^{\prime})=(\bm{z}^{\prime})^{\mathrm{T}}(B^{\prime})^{\mathrm{T}}X^{\mathrm{T}}X(B^{\prime})\bm{z}^{\prime}-2(\bm{z}^{\prime})^{\mathrm{T}}(B^{\prime})^{\mathrm{T}}X^{\mathrm{T}}\bm{y}.$

(20)

Note that the cost function is still in the QUBO formulation.

To reduce the redundant binary variables, we should find strongly correlated pairs in advance. Hence, we carry out a short sampling procedure for the original problem with the continuous variables with a conventional computer. We then choose correlated pairs that share the same binary variables. We determine the correlated pairs using the correlation matrix with a certain threshold.

Note that the threshold is determined by how many pairs we want to make; the limitation of the number of spins in the annealing machines affects this.

The procedure is summarized as follows:

1. (i)

   Sample short-time series data of continuous variables using Monte Carlo methods.

2. (ii)

   Calculate the correlations between continuous variables from the time series data.

3. (iii)

   Define a basis matrix $B^{\prime}$ and binary vector $\bm{z}^{\prime}$ using the information on the correlations. That is, reduce matrix $B$ and $\bm{z}$ for the variable pairs whose correlation values exceed a certain threshold.

4. (iv)

   Derive the QUBO formulation by using matrix $B^{\prime}$ and the reduced $\bm{z}^{\prime}$.

We give some comments on this procedure.

First, we need additional samplings in step (i). If the sampling step takes a long time, this procedure could not be practical even if the optimization is fast in annealing machines. However, we do not need highly accurate information on the correlations; only a rough estimation is sufficient for the pairing in step (ii). Hence, the short-time sampling is sufficient, and its computational cost is low.

Second, one may obtain the correlations analytically in some cases, for example, in Gaussian random Markov fields. In such cases, we can skip the Monte Carlo samplings in step (i).

Third, the number of reduced binary variables depends on the threshold with which we define the strongly correlated pairs. As stated above, we simply set the threshold in accordance with how much we want to reduce the number of variables.

We use the following linear function to generate artificial data:

	$\displaystyle f(\bm{x})=$	$\displaystyle 15.5+15.5x_{1}+10.0x_{2}+10.0x_{3}+5.0x_{4}+5.0x_{5}$
		$\displaystyle-0.5x_{6}-0.5x_{7}-15.5x_{8}-15.5x_{9}.$		(21)

Hence, the number of parameters is $D=10$. First, we generate $x_{d}\,(d=1,\dots,9)$ from a uniform distribution $[-1,1]$ and we set it as the $i$-th input data $\bm{x}_{i}$. Observation noise $\eta_{i}\sim\mathcal{N}(0,1)$ is then sampled from the standard normal distribution. By adding $\eta_{i}$ to the true output $f(\bm{x}_{i})$ as

$\displaystyle y_{i}=f(\bm{x}_{i})+\eta_{i},$

(22)

we have the $i$-th data $(\bm{x}_{i},y_{i})$.

We repeat this process $1000$ times and split the generated data into two sets, i.e., training and test. The size of the training data set is $100$, and the remaining $900$ data pairs are used as the test data set. Since the problem is not difficult to solve, we use a small amount of the training data.

Once we fix the training data, the cost function in Eq. (12) is determined. As the preparation step, we then obtain short-time series data from the Monte Carlo sampling for the linear regression cost function in Eq. (12), which has the original continuous variables. Samples are generated from the Gibbs-Boltzmann distribution with the aid of the Metropolis rule. We generate the candidates for the Metropolis rule by adding a random variable generated from $\mathcal{N}(0,0.5)$ to a randomly selected parameter. The sampling interval is $2\times D=20$ steps, and the length of the time-series data is $100$. The sampling temperature is $0.1$, although we confirmed that the other temperatures are also suitable to obtain a rough estimation of the correlations.

Figure 1 shows an example of correlation matrices. There are several correlated pairs. The highest correlated one is $(w_{1},w_{3})$. Although $(w_{1},w_{2})$ and $(w_{2},w_{3})$ are candidates, we do not use them because $w_{1}$ and $w_{3}$ appear in the pair $(w_{1},w_{3})$. We then select the next pairs $(w_{0},w_{5})$ and $(w_{8},w_{9})$. We set the threshold to $0.8$. Hence, we reduced the binary variables from these three pairs, i.e., $(w_{1},w_{2})$, $(w_{0},w_{5})$, and $(w_{8},w_{9})$ for the case in Fig. 1.

Although we do not know detailed information about parameters $\{w_{i}\}$ in advance, it is necessary to define a basis vector in Eq. (13). Since our aim was to show the effectiveness of the proposed method, we choose the basis vector $\bm{b}$ as follows:

$\displaystyle\bm{b}=\begin{bmatrix}0.5&-0.5&1&-1&2&-2&4&-4&8&-8\end{bmatrix}^{\mathrm{T}}.$

(23)

Hence, $K=10$ in Eq. (13). Note that the choice of the basis vector $\bm{b}$ is sufficient to represent the target function $f(\bm{x})$ in Eq. (21). The $\bm{b}$ yields the basis matrix $B$ via $B=I_{D}\otimes\bm{b}$.

Next, we redefine matrix $B^{\prime}$ with the aid of the above-mentioned highly correlated pairs. For the selected pairs, we share certain binary variables, as exemplified in Sec. IV.3; the order is from highest absolute value to lowest. In the experiment, we changed the number of shared binary variables to each pair and investigated the decrease in the prediction accuracy.

Finally, we redefine $\bm{z}^{\prime}$ in accordance with the redefined matrix $B^{\prime}$. We then obtain the final QUBO formulation in Eq. (20) for the annealing machines.

We use a simple SA algorithm to find optimal solutions $\bm{z}^{\prime}$ of the derived QUBO formulation. The SA simulations ran on a MacBookAir with Inter Core i5. In SA, the temperature schedule was as follows:

$\displaystyle T_{t}=T_{0}\times\gamma^{t},$

(24)

where $T_{t}$ is the temperature at the $t$-th iteration, $T_{0}$ is the initial temperature, and $\gamma$ is the decay rate. One iteration is defined as $2\times\textrm{(\# of spins)}$ updates in accordance with the conventional Metropolis rule. After the iterations, we output the state that minimizes the QUBO formulation during the search procedures. Table 1 shows the SA parameters.

The obtained $\bm{z}^{\prime}$ leads to the final estimation for the continuous variables $\{w_{i}\}$ via Eq. (19).