A quantum-inspired tensor network method for constrained combinatorial optimization problems

We use an objective function $f$ and a discrete set of feasible solutions x to specify the combinatorial optimization problem. In many such problems, each solution involves a binary selection of the individual components under certain constraints, denoted as $\textbf{x}=(x_{1},\cdots,x_{n})$ with binary variables $x_{n}\in\{0,1\}$, satisfying constraints $\textbf{c}=(c_{1},\cdots,c_{m})$. The goal is to find the solution $\textbf{x}^{\star}$ in all feasible solutions to maximize the objective function $f$, i.e. $\textbf{x}^{\star}=\mathop{\mathrm{arg\,max}}_{\textbf{x}\in\{0,1\}^{\otimes n}}f(\textbf{x})$. This goal can be achieved by mapping the problem to a Hamiltonian and looking for the ground state in the Hilbert space. The Hamiltonian is written as

$\displaystyle\hat{H}=\sum_{\mathbf{x}\in\{0,1\}^{\otimes n}}a_{\mathbf{x}}\ket{\mathbf{x}}\bra{\mathbf{x}},$

(1)

where $a_{\textbf{x}}$ are real values representing the elements of the Hamiltonian and $|\textbf{x}\rangle=|x_{1}\rangle\otimes|x_{2}\rangle\cdots\otimes|x_{n}\rangle$. Note that in this paper, we are interested in problems that require $n$ qubits, where $n$ is the number of the binary variables. Each variable $x_{i}$ of the set of solutions x is assigned on the basis of the Pauli operator $\hat{Z}_{i}$ on the $i$-th qubit, that is, $\hat{Z}_{i}=\ket{0}\bra{0}-\ket{1}\bra{1}$ with $\ket{x_{i}}\in\{\ket{0},\ket{1}\}$. Thus, a general quantum state can be expressed as

$\displaystyle\ket{\Psi}=\sum_{\mathbf{x}\in\{0,1\}^{\otimes n}}b_{\mathbf{x}}\ket{\mathbf{x}}.$

(2)

Here $b_{\mathbf{x}}$ represents the linear coefficient that meets the normalized condition, i.e. $\sum_{i=1}^{n}|b_{x_{i}}|^{2}=1$. Consider a general unconstrained binary linear optimization problem:

$\displaystyle\min\{\mathbf{w}^{T}\mathbf{x}:\mathbf{x}\in\{0,1\}^{\otimes n}\}$

(3)

with weights $\mathbf{w}=(w_{1},\cdots,w_{n})$. Then the Hamiltonian is transformed as

$\displaystyle\hat{H}$

$\displaystyle=\sum_{i}\frac{w_{i}}{2}(\hat{Z}_{i}-\hat{I}),$

(4)

where $I$ is the identity. Such transformations take linear time in the number of variables.
The ground state of this Hamiltonian, denoted as $\ket{\psi_{g}}$, is the state that minimizes the energy defined as $E(\ket{\Psi})=\bra{\Psi}\hat{H}\ket{\Psi}/\langle\Psi\ket{\Psi}$. With an appropriate mapping, the original optimization problem is equivalent to solving the ground state of the Hamiltonian. Then the ground-state solving techniques in quantum many-body physics can be utilized to achieve an optimized solution of the combinational optimization problem.

For constrained problems, additional terms are usually needed in the Hamiltonian to penalize constraint violations. These penalty terms should be conditioned by multiplying with appropriate penalty factors to ensure that the overall Hamiltonian behaves as intended. Optimization of such a Hamiltonian containing penalty terms inevitably involves tuning of additional hyperparameters, which greatly increases the difficulty of the overall optimization process.

In our study, we present a hyperparameter-free algorithm that directly projects a randomized initial quantum state into a subspace that satisfies the constraints. We use a projector, defined as $\hat{P}\ket{\Psi}=|\Psi_{v}\rangle$, to eliminate states that violate constraints. $\ket{\Psi}$ is an arbitrary state in the Hilbert space, and $|\Psi_{v}\rangle$ is the projected state that belongs to the subspace denoted by $V$ satisfying the constraints. The projector is then constructed as

$\displaystyle\hat{P}=\sum_{\mathbf{x}\in V}|\mathbf{x}\rangle\langle\mathbf{x}|.$

(5)

It satisfies $\hat{P}^{2}=\hat{P}$. The ground state $\ket{\psi_{g}}$ under the projection is written as follows:

	$\displaystyle\ket{\psi_{g}}$	$\displaystyle=\min_{\ket{\Psi}}\ E(\ket{\Psi}),$		(6)
	$\displaystyle E(\ket{\Psi})$	$\displaystyle=\frac{\bra{\Psi_{v}}\hat{H}\ket{\Psi_{v}}}{\braket{\Psi_{v}}{\Psi_{v}}}=\frac{\bra{\Psi}\hat{P}\hat{H}\hat{P}\ket{\Psi}}{\braket{\Psi}{\hat{P}}{\Psi}}.$		(7)

We can then extract the optimal solution to the combinatorial optimization problem encoded in $\ket{\psi_{g}}$. It is worth noting that the solver only works for a normalized state, i.e. $\langle\Psi\ket{\Psi}=1$, while the projected state are not normalized since a general projector $\hat{P}$ is not unitary. However, with prior knowledge that the ground state is always a product state following the constraints, and thus $\braket{\Psi_{0}}{\hat{P}}{\Psi_{0}}=1$, we could simply apply the solver for the normalized ground state of $\hat{P}\hat{H}\hat{P}$.

Directly constructing the projector $\hat{P}$ incurs an exponential cost in time and memory. However, by encoding the constraints in a tensor network state, we can significantly reduce the cost to $O(mn)$, or $O(m+n)$ if the constraints are “local”, i.e. the number of variables involved in each constraint is limited by a constant and vice versa. The general form of the tensor network state can be written as

$\displaystyle\ket{\Psi}=\sum_{\{\alpha,\beta\cdots\}}\prod_{i=1}^{n}T^{[x_{i}]}_{\alpha_{i}\beta_{1}\cdots\beta_{x}}R^{[c_{i}]}_{\beta_{1}\cdots\beta_{c}}|\alpha_{1}\cdots\alpha_{i}\cdots\alpha_{n}\rangle.$

(8)

There are two types of tensors, $T^{[x]}$ and $R^{[c]}$. Each $T^{[x]}$ has a single physical index labeled $\alpha$ with a bond dimension $d=2$ reflecting the binary variable, that is, $\alpha=0$ and $1$, and multiple virtual indices labeled $\beta$ also with the same bond dimension. Each auxiliary tensor $R^{[c]}$ encodes the constraint $c$ and connects to $T^{[x]}$ tensors whose variable $x$ is involved in the constraint $c$. To map the selection of the binary variable from the physical index to other indices, $T^{[x]}_{0\cdots 0}$ and $T^{[x]}_{1\cdots 1}$ are set to $1$. The $R^{[c]}$ tensors are constructed by traversing the allowed assignments of the constrained variables and setting the elements located in the corresponding indices to $1$. We initialize the $T^{[x]}$ and $R^{[c]}$ tensors as zeros. For example, if we have variables $x_{1},x_{2},x_{3}\in\{0,1\}$ and a single constraint $c_{1}:x_{1}+x_{2}+x_{3}\mod 2=1$, then we can use four sparse tensors, i.e. $T^{[x_{1}]}_{\alpha\beta_{1}},T^{[x_{2}]}_{\alpha\beta_{2}},T^{[x_{3}]}_{\alpha\beta_{3}}$ and $R^{[c_{1}]}_{\beta_{1}\beta_{2}\beta_{3}}$ to encode this particular constrained problem, with nonzero elements of the tensors specified as:

	$\displaystyle T^{[x_{1}]}_{00}=T^{[x_{1}]}_{11}=1$			(9)
	$\displaystyle T^{[x_{2}]}_{00}=T^{[x_{2}]}_{11}=1$			(10)
	$\displaystyle T^{[x_{3}]}_{00}=T^{[x_{3}]}_{11}=1$			(11)
	$\displaystyle R^{[c_{1}]}_{001}=R^{[c_{1}]}_{010}=R^{[c_{1}]}_{100}=R^{[c_{1}]}_{111}=1$	.		(12)

The four possibilities to satisfy the constraint $c_{1}$ have been included in $R^{[c_{1}]}$; for example, $R^{[c_{1}]}_{001}=1$ with indices $001$ represents $x_{1}=x_{2}=0$ and $x_{3}=1$. The construction of these tensors is visualized in Fig.1.

By this construction, we encode all feasible variable assignments in the initialized tensor network. Since traversing the allowed assignments of each constraint requires exponential time in the number of variables involved, the overhead of our method is much lower for problems with local constraints, compared to directly building the projector. Additionally, existing tensor algebra methods and tensor network algorithms are applicable to structured problems, helping us solve problems efficiently.

The optimal solution is one of the basis vectors in Eq. (8), denoted by $|\alpha_{1}\cdots\alpha_{i}\cdots\alpha_{n}\rangle$, with $\alpha_{i}=0,1$ on the $i$-th binary variable $x_{i}$. Having the superposition of all allowed variable assignments, we can screen out the optimal solution by imaginary time evolution (ITE), a technique to project the initial state to the ground state of an objective Hamiltonian. ITE effectively performs a power iteration by repetitively applying the ITE operator $e^{-\tau\hat{H}_{p}}$ to find the ground state of $\hat{H}$. The ground state minimizes the energy, which is equivalent to minimizing or maximizing the objective function. The simulation of the ITE takes the form of

$\displaystyle\ket{\Psi}=\lim_{\tau\rightarrow\infty}\frac{e^{-\tau\hat{H}}\ket{\Psi_{0}}}{||e^{-\tau\hat{H}}\ket{\Psi_{0}}||},$

(13)

where $\tau$ is referred to as the imaginary time. If the Hamiltonian only contains the summation of commuting terms, then $e^{-\tau\hat{H}}$ can be directly rewritten into the product of the local evolution operators, i.e. $e^{-\tau\hat{H}}=\Pi_{i}e^{-\tau\hat{h}_{i}}$. Otherwise, the Trotter-Suzuki decomposition[19, 20], an approximation for doing the ITE if containing non-commuting terms in the Hamiltonian, should be involved.

To extract the optimal solution from the tensor networks, we take inspiration from measuring the spin magnetic moment. The assignment of the variable $x_{i}$ is calculated as

$\displaystyle x^{\star}_{i}=\mathbbm{1}(\frac{\braket{\Psi}{\hat{Z}_{i}}{\Psi}}{\braket{\Psi}{\Psi}}<0),$

(14)

where $\hat{Z}_{i}$ is the Pauli matrix. That is, if the expectation value is negative, $x_{i}$ is $1$; otherwise, it is $0$. Since we need to do this calculation for each variable, there is an $O(n)$ factor multiplied with the complexity of computing a single expectation value, which we will discuss shortly.

Note that there could be certain variables whose either assignment gives the same final objective value and obeys the constraints, known as the degeneracy. In this case, if we prefer $x_{i}$ to be assigned a specific value, for example, zero, then the measurement operator can be adjusted as

$\displaystyle\hat{\tilde{Z}}_{i}=\begin{bmatrix}1&0\\ 0&-\mu\end{bmatrix},$

(15)

where $\mu$ is a value larger than $1$ so as to create a slight difference between the expectation values of the $0$ and $1$ variable assignments. Then we should measure

$\displaystyle\tilde{x}^{\star}_{i}=\mathbbm{1}(\frac{\braket{\Psi}{\tilde{\hat{Z}}_{i}}{\Psi}}{\braket{\Psi}{\Psi}}<\nu),$

(16)

where $\nu$ is a threshold value related to $\mu$ used to detect the slight difference. In practice, $\mu$ can be slightly larger than $1$ and $\nu$ is a small positive number approaching zero, such as $\mu=3$ and $\nu=0.04$, both depending on how many degenerate states exist. We can adjust $\mu$ until the degenerate states are properly separated, then we can find our preferred one by setting an appropriate threshold $\nu$.

We calculate the expectation values $\bra{\Psi}\tilde{\hat{Z}}_{i}\ket{\Psi}$ by performing tensor network contractions [21, 22, 23], which sometimes take exponential time and space to obtain an exact result. Tensor decomposition methods, such as singular value decomposition (SVD), are regularly used in tensor network contractions to limit the rapidly increasing bond dimension. Specifically, we can preserve a predetermined number of columns and rows in the decomposed matrices with the largest singular values. An alternative common practice is to preserve columns and rows whose singular values are above a threshold. Using these techniques cleverly can reduce the complexity of the whole contraction to be polynomial in $m$ and $n$, while still resulting in an accurate enough value. Moreover, for problems that have a structured tensor network construction, we can utilize existing tensor network algorithms to further optimize performance.

We use a real-world combinatorial optimization problem, the open-pit mining problem, as an example to demonstrate our algorithm. The goal of designing optimal open-pit mines is to maximize ore production while avoiding unnecessary mining of rocks. The planning is also subject to a variety of constraints on the size, shape, and form of the mine, making it a computationally taxing process. Theoretical studies of the problem often apply simplifying assumptions, converting the problem into a simpler and more well-understood form [24].

In particular, the 2D open-pit mining problem can be formally stated as a combinatorial optimization problem. Consider a 2D square lattice of mining blocks, where each block $i$ has an associated profit $w_{i}$. The coordinate of block $i$ is denoted as $(i_{x},i_{y})$. A feasible solution should follow physical constraints, i.e. if the block $i$ is excavated, then all its child blocks $j$ should be excavated as well. In this work, we consider the $45^{\circ}$ slope constraint: $j\in\{(i_{x}-1,i_{y}-1),(i_{x}-1,i_{y}),(i_{x}-1,i_{y}+1)\}$, as illustrated in Fig. (2) (a). Equivalently, we can consider an $n$-node directed graph $G=(V,E)$ with node weight $w_{i}$, as shown in Fig. (2) (b). $G$ is structured into levels, where each level contains two more nodes than the previous level, starting from one node in the first level. Each node before the last level has exactly three child nodes at the next level, and each node after the first level has one to three parent nodes.

By assigning $x_{i}=0$ (unexcavated) or $x_{i}=1$ (excavated) to each node $i$, we can write the open-pit mining problem as

$\displaystyle\max\sum_{i=1}^{n}w_{i}x_{i}$

(17)

subject to

	$\displaystyle x_{i}\in\{0,1\}$	$\displaystyle\text{\ \ \ for }i=1,2,\cdots,n$		(18)
	$\displaystyle\sum_{j\in\mathit{children}(i)}x_{i}(1-x_{j})=0$	$\displaystyle\text{\ \ \ for }i=1,2,\cdots,n,$		(19)

where $\mathit{children}(i)$ denotes the set of child nodes of node $i$.

Traditionally, the open-pit mining problem is solved by reducing to the maximum closure problem or the maximum flow problem and utilizing efficient graph algorithms. In particular, the Lerchs-Gorssman (LG) algorithm [25, 26] was the most widely used algorithm in the mining industry, giving a provably optimal solution in polynomial time. In recent years, it has been surpassed by the more efficient Pseudoflow algorithm [27, 28], a $O(|V||E|\log|V|)$ algorithm for the maximum flow problem. Recently, a quantum computing approach was proposed as the first attempt to solve this problem with quantum computers [29]. It modifies the objective function to

$\displaystyle\max(\sum_{i}w_{i}x_{i}-\lambda\sum_{i,j\in p(i)}x_{i}(1-x_{j})),$

(20)

where $\lambda$ is a hyperparameter introduced to regularize the penalty for constraint violations. Then the problem Hamiltonian is constructed using the method described in Section II.1.

Our algorithm takes inspiration from the quantum computing approach but is intrinsically different. Using tensor networks to represent quantum states, our algorithm can employ powerful non-unitary operations, which are unavailable to quantum algorithms. This fact allows us to directly construct the superposition state of all feasible solutions and avoids having to optimize the regularization of the penalty term. Our algorithm is also completely different from the graph algorithms. It provides a new perspective on such problems, utilizing the rapidly developing tensor network methods. Moreover, the core ideas can be applied to general combinatorial problems, not limited to just this one.

We construct the configurations of all allowed states obeying the smoothness constraints as a tensor network state. Fig.3 visualizes the construction process using the $5\times 3$ mine as an example. Referring to the smoothness constraints that the walls of the pit should not exceed a maximum steepness, the sites marked dark brown in Fig.2(a) would, if excavated, violate the smoothness constraints, which therefore should be excluded from the pit profile as well as the initialized tensor network construction. The values of the tensors as defined in Fig.3(c) are as follows (here we group the tensors corresponding to different variables or constraints but have the same form together):

	$\displaystyle T^{[1]}_{11}=T^{[1]}_{00}=1$		(21)
	$\displaystyle T^{[2]}_{1111}=T^{[2]}_{0000}=1$		(22)
	$\displaystyle T^{[3]}_{111}=T^{[3]}_{000}=1$		(23)
	$\displaystyle T^{[4]}_{111}=T^{[4]}_{000}=1$		(24)
	$\displaystyle T^{[5]}_{11111}=T^{[5]}_{00000}=1$		(25)
	$\displaystyle R^{[1]}_{1111}=R^{[1]}_{1110}=R^{[1]}_{1101}=R^{[1]}_{1100}=R^{[1]}_{1000}=R^{[1]}_{0100}=R^{[1]}_{0000}=1$		(26)
	$\displaystyle R^{[2]}_{111}=R^{[2]}_{110}=R^{[2]}_{100}=R^{[2]}_{010}=R^{[2]}_{000}=1.$		(27)

All remaining elements are zero, which means that the corresponding projections are forbidden. The first index of the $T$ tensors is the physical index, where $\alpha=0$ represents an unexcavated block and $\alpha=1$ represents an excavated block. The virtual indices other than the physical one of $T$ are used to transfer the status of neighboring tensors. As in the illustration of the tensors shown in Fig.3(c), the arrow on each index is used to distinguish the directions of the transferring status of the geometrical bonds, i.e. the indices with incoming arrows are carrying the status originated from their parent blocks, and the indices with outgoing arrows are carrying the status, which should be the same as their physical indices and will head to their child blocks. The constraints of the excavated ore have been reflected in our tensor network definition with very limited nonzero elements of the tensor. Under the constraint that if a block $(i_{x},i_{y})$ is excavated, so must its parent blocks $(i_{x}-1,i_{y}-1)$, $(i_{x}-1,i_{y})$ and $(i_{x}-1,i_{y}+1)$, but an excavated block itself does not have to have child blocks.

Since there is only a very limited number of nonzero elements in the initialized tensor network, there is a lot of room for optimization in terms of computational memory and time cost. In addition, the three-dimensional open-pit mine can also be defined in a similar way, just by adjusting the local tensors $T$ and the auxiliary tensors $R$ to higher-order tensors and designing the nonequivalent tensor for boundary conditions in the three-dimensional mine.

The schematic of ITE is shown in Fig.4a. As the Hamiltonian for the open-pit mining problem, the same as defined in Eq.(4), only contains single site commuting terms, we can rewrite the evolution operator for the whole system, i.e. $e^{-\tau\hat{H}}$ in Eq.(13) into the product of the local evolution operators, i.e. $e^{-\tau\hat{H}}=\Pi_{i}e^{-\tau\hat{h}_{i}}$, with $\hat{h}_{i}=\frac{w_{i}}{2}(\hat{Z}_{i}-\hat{I})$. Then directly set the imaginary time $\tau$ as a relatively large number, such as $\tau=10$, to achieve the ground state in the complexity of $O(n)$.