QUBO Formulations for a System of Linear Equations

With the advent of quantum computers, many quantum computing algorithms are being developed. Solving linear systems is one of the most fundamental problems in almost all science and engineering. The Harrow-Hassidim-Lloyd algorithm, a monumental quantum algorithm for solving linear systems on gate model quantum computers, was invented and several advanced variations have been developed. For a given square matrix 𝐴 ∈ ℝ 𝑛×𝑛 and a vector 𝑏⃗ ∈ ℝ 𝑛 , we will find unconstrained binary optimization (QUBO) models for a vector 𝑥 ∈ ℝ 𝑛 that satisfies 𝐴𝑥 = 𝑏⃗ . To formulate QUBO models for a linear system solving problem, we make use of a linear least-square problem with binary representation of the solution. We validate those QUBO models on the D-Wave system and discuss the results. For a simple system, we provide a Python code to calculate the matrix characterizing the relationship between the variables and to print the test code that can be used directly in the D-Wave system.


Introduction
Quantum computing has opened up a new paradigm for approaching computing problems where classical (i.e., traditional) computing can provide unparalleled speeds for specific problems. A specific subset of quantum computing is quantum annealing aimed at optimization problems. Quantum annealing processors naturally return low-energy solutions; some applications require the real minimum energy and others require good low-energy samples 1 . One popular model of optimization problems on which quantum annealers are based is the Ising model 2 .
A monumental quantum algorithm for solving linear systems on gate model quantum computers, was invented in 2008, and a variety of advanced variations were developed 3 . At that time, their algorithm and its variants had many limitations, such as the quantum state of the solution, the need for quantum RAM, and the limitations of Matrix A 4 . In particular, the quantum state of the solution limits the use of the algorithm to submodules of the entire system, or reduces the quantum speed improvement in real world implementations, as it requires iterative computations to read the quantum state of each base. In 2019, a method of using linear least squares was proposed to solve a linear problem [5][6][7] . They suggested a QUBO model of linear equations for by linear systems. In addition, when solving a linear problem using the QUBO model, they compared several classic methods such as QR factorization and SVD 8 . There has been remarkable research [9][10][11][12][13][14] since the development of D-Wave quantum annealers.
In this paper, we propose two representative QUBO models for by linear systems. Additionally, we show the overall CPU calculations when performing the new QUBO modeling. The amount of calculation used here consists of simple addition and multiplication. These results will be of great advantage when solving linear equations for large . We also provide Python code to create QUBO models that can be used in D-Wave quantum annealers. To date, the number of qubits that can be used in quantum annealers is approximately 200, so the size of the linear problem that can be solved using our QUBO model is not large. However, with the development of quantum annealers, the method we propose will have advantages. It is also expected that other issues could be approached in a similar way to ours.

Method Background
The Ising model is a mathematical model for ferromagnetism in statistical mechanics. The energy Hamiltonian (the cost function) is formulated as follows:
Quadratic unconstrained binary optimization (QUBO) is a combinatorial optimization problem in computer science. In this problem, a cost function is defined on an -dimensional binary vector space onto ℝ.
where is an upper diagonal matrix, = ( 1 , ⋯ , ) , and ∈ {0,1}. The problem is finding * , which minimizes the cost function . Since we have 2 = , the cost function is reformulated as follows:

Equation 3
In , the diagonal terms , and the off-diagonal terms , represent the linear terms and the quadratic terms, respectively. The unknowns of the Ising Model and the unknowns of the QUBO Model have the linear relation Therefore, we can choose appropriate format depending on the binary unknown.
Given a matrix ∈ ℝ × and a column vector of variables ∈ ℝ and a column vector ⃗ ∈ ℝ . The linear least squares problem is to find the that minimizes ∥ − ⃗ ∥. In other words, it can be described as follows: To solve Eq.

Equation 6
Taking the 2-norm square of the resultant vector of Eq. 6, we get the following:

Equation 10
where positive integer is the number of digits of and negative is the number of digits of fractional terms. Now, we can represent real value is below

Equation 11
To represent both positive and negative numbers, , + and , − are involved. Note that this representation can have the same value with different binary combinations.
To drive QUBO model, we insert Eq. 11 into Eq. 9. We obtain the two summation terms of the first term in Eq. 9 as below: Since is a positive, zero, or negative number, , 1 + × , 2 − for the same can be zero. We can get Eq. 13 from Eq.
12. Additionally, since 2 = , Eq. 15 can be finally obtained. Here, the front summation represents a linear term, and the second summation represents a quadratic term.
The QUBO form of the second term in Eq. 9 can be obtained as below:

Equation 17
When we calculate the QUBO term, we use an upper triangular matrix. Therefore , + , , − , , + , and , − have separate indices, and we can derive Eq. 17 from Eq. 16. This equation is part of the quadratic terms in our QUBO model.
The QUBO form of the third term in Eq. 9 can be obtained as below: The above equation is part of the linear terms in our QUBO model.
The last term in Eq. 9 represents the minimum value when the vector satisfies the linear equation in QUBO form. Therefore, our first QUBO model for linear systems is the summation of Eq. 15, Eq. 17, and Eq. 18.

QUBO modeling 2
To reduce the number of qubits to be used in Eq. 11, the new approximation of was introduced as follows:

Equation 19
Since the coefficient of − represents the amount of translation of the range of the summation part in Eq. 19, any coefficient of − that can represent the range of we want can be used. We select the coefficient of − , which represents the largest range that can have.
To drive QUBO model, we insert Eq. 19 into Eq. 9. We get another QUBO form of the first term in Eq. 9 as below:

Equation 23
Since 2 = , Eq. 23 can be finally obtained. The front summation represents linear terms, and the second summation represents quadratic terms.
The QUBO form of the second term in Eq. 9 can be obtained as below:

Equation 27
The last term in Eq. 9 represents the minimum value when the vector satisfies the linear equation in QUBO form.
Therefore, our second QUBO model for linear systems is the summation of Eq. 23, Eq. 26, and Eq. 27.

Implementation
Let's consider the following example 16 : where is binary and the range of integer is −3 ≤ ≤ 3.
They used the related matrix below to solve Eq. 28.

Equation 29
The diagonal entries and the upper diagonal entries represent coefficients of linear terms and quadratic terms, respectively, in Formula 2 .
Our matrix from QUBO modeling 1 is below:

Equation 30
Two different matrices have the same solution for Eq. 28.

Empirical Test
We tested the QUBO model of Eq. 30 on the D-Wave 2000Q system with 10,000 anneals, and compared the result from Eq. 29 16      For 2 with zero terms in QUBO, Table 2 shows the number of occurrences of the lowest energy −26.
The total occurrence of all combinations of qubits having the lowest energy is 8,400 from 10,000 anneals. Thus, 8,400 anneals achieved the lowest energy state, with an energy equal to -26.   Table 3 shows the number of occurrences of the lowest energy −26. The total occurrence of all combinations of qubits having the lowest energy is 8,725 from 10,000 anneals.
Thus, 8,725 anneals achieved the lowest energy state, with an energy equal to -26.

Discussion
The most computationally intensive part of our first QUBO form is Eq. 17, one of the quadratic terms.
We want to discuss Eq. 17 below. By computing Eq. 17, we calculate variables , , and . Variable can vary from 1 to , and and represent the ( , ) components of the upper triangular matrices without diagonal components, so the total amount of calculation can be expressed as follows: Major computing cost for ( , ) = ( − )

Equation 33
To make the QUBO model solve the linear equation, the total amount of calculation seems to be large as shown in

Equation 34
In Eq. 34, the coefficient of each quadratic term is divisible with respect to . This means that it can be efficiently partitioned when distributing n parallel processes. At this time, each process requires computing cost as in Eq. 32.
To calculate the remaining term the fastest, it can be calculated through ( − 1)/4 parallel systems for each . In this case, the efficiency of the total computational cost is reduced, but the computational cost of the QUBO model is approximately 2 2 . When performing our QUBO modeling, as the total number of processes increases, the computing cost decreases. However, there is a large drawback to calculating linear equations in a quantum annealer.
Our QUBO formula requires that one qubit be connected to all other qubits. However, each quantum annealer has its own topology, which limits the connectivity between qubits. To develop a QUBO model that can be used in a quantum annealer where the connection between qubits is limited, it is necessary to study the method of diagonalizing 2 based on a quantum annealer's topology.
QUBO modeling 1 has a major disadvantage in that it uses approximately twice as many qubits as QUBO modeling 2. However, we mainly discussed the results of QUBO modeling 1 to explain the important results for the constraints of the QUBO models. To solve the system of linear equations, each element of can be divided by positive, negative, or zero. When Eq. 11 is used for the element , and the result can be expressed as combinations of elements in + and − . However, a solution can be made with either ∑2 , + or −∑2 , − . Because of this property, we can assume this problem to be a constrained problem for certain quadratic terms. To make certain quadratic terms 0, it is common to add the coefficients of the terms to a large value. To make certain quadratic terms equal to 0, it is common to add the coefficients of the terms to a large value. However, in this paper, we propose a new method of adding different values to each term so that it is set to 0 when the coefficient of each term is negative.
Therefore, we set each coefficient of , 1 + , 2 − to zero for each . Equation 30 was made by adding the same positive value to the coefficient of each negative constrained term in Eq. 29. In this way, we experimented by adding values of 100, 500, and 1000 to the constraint terms of Q matrices of various sizes with a size of 32 by 32 or less. However, we did not obtain better results than the two methods of setting each constrained term to 0 introduced in this paper.
We obtained these results in D-wave 2000Q. The reason why we obtained the best results by removing the constrained terms from QUBO is probably because we obtained a physical advantage by simplifying the superposition and entanglement of qubits by reducing the terms.