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

Combinatorial optimization is of general interest for both theoretical study and real-world applications. Fast-developing quantum algorithms provide a different perspective on solving combinatorial optimization problems. In this paper, we propose a quantum-inspired tensor-network-based algorithm for general locally constrained combinatorial optimization problems. Our algorithm constructs a Hamiltonian for the problem of interest, effectively mapping it to a quantum problem, then encodes the constraints directly into a tensor network state and solves the optimal solution by evolving the system to the ground state of the Hamiltonian. We demonstrate our algorithm with the open-pit mining problem, which results in a quadratic asymptotic time complexity. Our numerical results show the effectiveness of this construction and potential applications in further studies for general combinatorial optimization problems.


I. INTRODUCTION
Combinatorial optimization is the process of finding an optimal object from a discrete and finite set of objects.
Combinatorial optimization has extensive applications in almost every field of industry, such as supply chain optimization [1], transportation networks and power grids [2], and finance [3].The search space of a combinatorial optimization problem increases rapidly with the problem size.Problems like the Boolean satisfiability problem can have an exponentially large solution space, making an exhaustive search inapplicable for large-scale problems.From the complexity theory perspective, many combinatorial problems fall into the class of NP-hard, which is generally believed to be unsolvable in polynomial time on classical computers.Classical algorithms often employ heuristics and approximations to find nearly optimal solutions [4].Quantum algorithms [5], on the other hand, harness the power of randomness, superposition, entanglement, and interference from quantum mechanics, which might lead to an advantage in exploring the solutions of a combinatorial optimization problem [6][7][8].The implementation of quantum algorithms is currently limited by small-scale, noisy, and error-prone contemporary hardware [9]; nevertheless, they view the problems from a different perspective, motivating many quantum-inspired classical algorithms to appear [10,11].
Tensor networks (TNs) have undergone rapid development in the last two decades, gaining tremendous success in quantum many-body physics, quantum information sciences, statistical physics, and so on.Tensor network algorithms based on matrix product states (MPS) [12,13], projected entangled pair states (PEPS) [14,15], and variational renormalization group methods [16][17][18] are very efficient in simulating a large class of quantum many-body systems.
The tensor network structure can encode the combinatorial optimization problem with local constraints, providing an idea of utilizing the tensor network to solve combinatorial optimization problems.This paper presents a quantum-inspired tensor network algorithm to solve constrained combinatorial optimization problems and demonstrates the algorithm with a particular problem with numerical results.The paper is structured as follows: The general quantum-inspired tensor network algorithm is proposed in Section 2, and the open-pit mining problem is provided as an example in Section 3. Section 4 shows our numerical results with the quantum-inspired tensor network algorithm and concludes with a discussion of open questions and directions for future work.

II. A GENERAL TENSOR NETWORK ALGORITHM FOR COMBINATORIAL OPTIMIZATION
In this section, we divide our algorithm into four components and describe each in detail.Section II A explains how to construct the Hamiltonian for a classical combinatorial problem to transform it to a quantum problem.Section II B serves as inspiration for our core idea, Section II C, where the former studies how to satisfy the constraints without introducing a penalty term, and the latter details how to construct a tensor network state that represents the superposition of all feasible solutions.Finally, in Section II D, we show how to find the optimal solution by evolving the tensor network state to the ground state of the problem Hamiltonian.

A. Hamiltonian construction/Problem mapping
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 The goal is to find the solution x in all feasible solutions to maximize the objective function f , i.e. x = arg max x∈{0,1} ⊗n f (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 where a x are real values representing the elements of the Hamiltonian and 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 Ẑi on the i-th qubit, that is, Thus, a general quantum state can be expressed as Here b x represents the linear coefficient that meets the normalized condition, i.e.
where I is the identity.Such transformations take linear time in the number of variables.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 P |Ψ = |Ψ v , to eliminate states that violate constraints.|Ψ is an arbitrary state in the Hilbert space, and |Ψ v is the projected state that belongs to the subspace denoted by V satisfying the constraints.The projector is then constructed as It satisfies P 2 = P .The ground state |ψ g under the projection is written as follows: We can then extract the optimal solution to the combinatorial optimization problem encoded in |ψ g .It is worth noting that the solver only works for a normalized state, i.e.Ψ |Ψ = 1, while the projected state are not normalized since a general projector P is not unitary.However, with prior knowledge that the ground state is always a product state following the constraints, and thus Ψ 0 | P |Ψ 0 = 1, we could simply apply the solver for the normalized ground state of P Ĥ P .

C. Tensor network construction of the projected state
Directly constructing the projector 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 There are two types of tensors, T [x] and R [c] .Each T [x] has a single physical index labeled α with a bond dimension d = 2 reflecting the binary variable, that is, α = 0 and 1, and multiple virtual indices labeled β 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, 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 ∈ {0, 1} and a single constraint c 1 : β1β2β3 to encode this particular constrained problem, with nonzero elements of the tensors specified as: The four possibilities to satisfy the constraint c 1 have been included in R [c1] ; for example, R [c1] 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.

D. Finding the optimal solution
The optimal solution is one of the basis vectors in Eq. ( 8), denoted by 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 −τ Ĥp to find the ground state of Ĥ.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 where τ is referred to as the imaginary time.If the Hamiltonian only contains the summation of commuting terms, then e −τ Ĥ can be directly rewritten into the product of the local evolution operators, i.e. e −τ Ĥ = Π i e −τ ĥ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 where Ẑ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 where µ 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 where ν is a threshold value related to µ used to detect the slight difference.In practice, µ can be slightly larger than 1 and ν is a small positive number approaching zero, such as µ = 3 and ν = 0.04, both depending on how many degenerate states exist.We can adjust µ until the degenerate states are properly separated, then we can find our preferred one by setting an appropriate threshold ν.
We calculate the expectation values Ψ| Zi |Ψ 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.

III. THE OPEN-PIT MINING PROBLEM A. Problem description
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 • slope constraint: j ∈ {(i x − 1, i y − 1), (i x − 1, i y ), (i x − 1, i y + 1)}, as illustrated in subject to j∈children(i) where 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 max( where λ is a hyperparameter introduced to regularize the penalty for constraint violations.Then the problem Hamiltonian is constructed using the method described in Section II A.
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.

B. The tensor network framework
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 × 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): 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 α = 0 represents an unexcavated block and α = 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 −τ Ĥ in Eq.( 13) into the product of the local evolution operators, i.e. e −τ Ĥ = Π i e −τ ĥi , with ĥi = wi 2 ( Ẑi − Î).Then directly set the imaginary time τ as a relatively large number, such as τ = 10, to achieve the ground state in the complexity of O(n).We implement our algorithm for the mining problem as an open source Python library, available at [30].Specifically, the tensor network construction described in Section III B is a 2D tensor network, which can be viewed as projected entangled pair states (PEPS) [15], with trivial physical indices on the R tensors.We employ Koala [22], a highperformance PEPS simulation library, to perform the imaginary time evolution and calculate the expectation values by contractions.
We perform numerical experiments on open-pit mining problems with different size scales, as defined in Section III A. Problem instances are generated randomly, where the value of each site follows a normal distribution with a mean of 0.1 and a standard deviation of 1.Each problem size is repeated five times with different mine instances to account for the possible fluctuating running condition of the computing device.In an ideal world, the runtime should be the same given a fixed problem size since the algorithm is purely deterministic.The ITE is conducted with τ = 6, which is large enough to achieve the ground state, but not too large to cause numerical instability, such as exceeding the maximum allowed value of the complex type in Python.Note that one can also set τ with a heuristic related to problem size to better achieve the goals.We obtain the solutions by calculating the expectation values of the Pauli Z operators following Eq.( 14).Fig. 5 shows a comparison of our algorithm's runtime with different contraction approaches.If no bond dimension truncation is desired, the "PEPS Exact" method in the figure contracts a PEPS in a generally optimal order [23].
The boundary matrix product state (BMPS) method has poorly scaled performance on its own, but can be executed with bond dimension truncation to significantly lower time and memory complexities [22].Following the references, we can calculate the asymptotic time complexities of the three contraction approaches.While "PEPS Exact" and "BMPS no truncation" both take 2 O( √ n) , BMPS with truncation only takes O(d 6 n) = O(n), where n is the number of sites.Since we need to do n full contractions, the total time to calculate the expectation values with truncation is O(n 2 ).As other parts of the algorithm all run in linear time for the open-pit mining problem, the overall time complexity is O(n 2 ).In comparison, the state-of-the-art classical algorithm Pseudoflow finds the optimal solution in O(|V ||E| log |V |) = O(n 2 log n) [27].Note that although our algorithm is asymptotically faster, its constant factor is larger, and it makes approximations when truncating bond dimensions.In Fig. 5, we see that truncating the bond dimension to d = 2 incurs an additional overhead but has a polynomial runtime in contrast to the exponential runtime of both exact methods, which agrees well with our theoretical analyses.We compare the results produced with the optimal solution given by Pseudoflow.For all problem instances that we run, given an appropriate ITE evolution time, our algorithm can always find the optimal solution with maximized profit and zero constraint violations.Thus, we conclude that truncating the bond dimension to as small as d = 2 does not affect the accuracy of our algorithm.
It is worth noting that the open-pit mining problem has a few desirable properties for our algorithm.First, the system has a large energy gap between the ground state and the first excited state, allowing ITE to efficiently separate the optimal solution.Second, local and structured constraints decrease the complexity of constructing the tensor networks and make them reducible to well-studied tensor network forms.Third, the Hamiltonian contains only local interactions, resulting in well-controllable long-range entanglement, thus allowing aggressive bond-dimension truncations.It would be intriguing to see how the algorithm performs for more complicated and less suitable problems.Furthermore, by increasing the initial bond dimension, our algorithm can be extended beyond binary variables.A study on the performance sacrifice due to the increased initial bond dimension would be of interest.Moreover, as most of the tensor elements are zero in our construction, sparse tensor techniques could be applied, in addition to parallel computing implementations, to further improve performance.
V. ACKNOWLEDGEMENT

FIG. 1 .
FIG. 1. Tensor network construction for a three-variable constrained problem.(a) The schematic for the threevariable constrained problem with x1, x2, x3 ∈ {0, 1} and the constraint c1.(b) The tensor definitions in the general form of the tensor network for the three-variable constrained problem, with α as the physical index and β as the virtual index.

FIG. 2 .
FIG. 2. Two-dimensional open-pit mining problem and the equivalent graph representation.(a) A schematic of two-dimensional open-pit mining problem.The blocks in light brown are excavatable.The blocks in dark brown are unexcavatable due to 45 • slope constraint.ix and iy represent the vertical and horizontal indices, respectively.(b) The schematic for a directed graph representing a problem equivalent to the two-dimensional open-pit mining in panel (a).Each node i = (ix, iy) represents a mining block.Each directed edge (i, j) ∈ V represents the physical constraint: if block i is excavated, block j should be excavated too.

FIG. 3 .FIG. 4 .
FIG. 3. Tensor network construction for open-pit mining problem.(a) Directed graph representation of the openpit mining problem, with the excavatable blocks expressed as nodes in light brown.(b) A corresponding tensor network construction, where the light brown blocks show the physical nodes and the gray blocks show the virtual nodes.(c) The definition of nonequivalent tensors in our tensor network construction.α and β represent the physical and virtual indices, respectively.
FIG. 5. Computational time for two-dimensional open-pit mining problems with mine side length ranging from 3 to 13, using three different PEPS contraction approaches.
The ground state of this Hamiltonian, denoted as |ψ g , is the state that minimizes the energy defined as E(|Ψ ) = Ψ| Ĥ |Ψ / Ψ |Ψ .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.