Quantum-Inspired Algorithm for Classical Spin Hamiltonians Based on Matrix Product Operators

Abstract: We propose a tensor-network (TN) approach for solving classical optimization problems that is inspired by spectral filtering and sampling on quantum states. We first shift and scale an Ising Hamiltonian of the cost function so that all eigenvalues become non-negative and the ground states correspond to the the largest eigenvalues, which are then amplified by power iteration. We represent the transformed Hamiltonian as a matrix product operator (MPO) and form an immense power of this object via truncated MPO-MPO contractions, embedding the resulting operator into a matrix product state for sampling in the computational basis. In contrast to the density-matrix renormalization group, our approach provides a straightforward route to systematic improvement by increasing the bond dimension and is better at avoiding local minima. We also study the performance of this power method in the context of a higher-order Ising Hamiltonian on a heavy-hexagonal lattice, making a comparison with simulated annealing. These results highlight the potential of quantum-inspired algorithms for solving optimization problems and provide a baseline for assessing and developing quantum algorithms.

• The paper introduces a Trotter-free spectral-amplification algorithm that shifts a diagonal cost Hamiltonian, repeatedly powers its MPO, embeds it into an MPS, and samples near-optimal configurations.
• The method improves with iteration count and bond dimension, reaching approximation ratios up to 0.999 on tested Edwards–Anderson instances while avoiding the stagnation observed in uniform-state DMRG.
• The benchmarks show that linear MPO powering can outperform simulated annealing on higher-order Ising problems with up to roughly 580 variables, although bond-dimension growth and spectral-gap dependence limit scalability.

The paper introduces a classical tensor-network algorithm for combinatorial optimization problems (COPs) formulated as diagonal spin Hamiltonians. Rather than variationally minimizing the cost function via the density-matrix renormalization group (DMRG), Watanabe, Tindall, Miyakoshi, and Ueda perform spectral amplification directly on the cost Hamiltonian itself: after a shift and rescale that renders it positive semidefinite, they compute a high power of its matrix product operator (MPO) representation through truncated MPO–MPO contractions, embed the result into a matrix product state (MPS), and sample from it in the computational basis (2602.05224). The method is Trotter-free, avoids gate decompositions required in quantum imaginary time evolution, and is benchmarked against DMRG on Edwards–Anderson spin glasses and against simulated annealing (SA) on higher-order Ising problems defined on IBM's heavy-hexagonal lattice.

Spectral filtering via powered cost operators

The starting point is a cost function $C(\bm z)$ over $d_i$-ary strings, encoded as a Hamiltonian $\hat H$ diagonal in the computational basis and expanded in an orthonormal operator basis built from traceless, Hermitian, diagonal generators of $\mathrm{SU}(d_i)$ (Pauli-$Z$ strings for qubits). The authors transform this into

$\hat G = -\hat H + \Lambda\,\mathbb I,$

with $\Lambda$ exceeding the spectral radius of $\hat H$, so that $\hat G$ is positive semidefinite and its largest eigenvalues correspond to the ground states of $\hat H$. Since $d_i$0 remains diagonal, it can be applied to the uniform product state $d_i$1 to define a pure state whose sampling probabilities are proportional to $d_i$2. In the large-$d_i$3 limit this distribution concentrates entirely on the optimal configurations $d_i$4.

A key theoretical contribution is the entanglement analysis of the resulting state. Because $d_i$5 is positive semidefinite, growth of the bipartite entanglement entropy is suppressed compared to naive power methods; asymptotically $d_i$6 converges to a sum of $d_i$7 rank-one projectors onto ground-state configurations, implying a bond dimension bounded by the degeneracy $d_i$8. For a worst-case two-level cost spectrum separated by gap $d_i$9, the number of iterations needed to sample the optimum with probability $\hat H$0 satisfies

$\hat H$1

so the method scales inversely with the spectral gap and only logarithmically with the configuration space dimension — but this scaling estimate rests on a deliberately pessimistic flat-spectrum assumption, and the practical dependence of the truncation-induced bond dimension on problem structure remains uncharacterized for generic instances.

Algorithmic construction

The MPO bond dimension of the shifted operator is tightly bounded: for arbitrary QUBO instances with coupling-matrix rank $\hat H$2, the maximum bond dimension obeys $\hat H$3, derived from the operator Schmidt rank of the cross-bipartition coupling block $\hat H$4. For geometrically local lattices mapped to a chain, $\hat H$5 grows as $\hat H$6, reflecting the hypercube boundary rather than the volume.

Two contraction schedules are considered. A "linear" schedule applies $\hat H$7 repeatedly ($\hat H$8 multiplications) at cost $\hat H$9 per step, matching DMRG's scaling; a "doubling" schedule reaches $\mathrm{SU}(d_i)$0 in $\mathrm{SU}(d_i)$1 squarings at cost $\mathrm{SU}(d_i)$2. All contractions use the zip-up algorithm with full SVD truncation to a maximum bond dimension $\mathrm{SU}(d_i)$3 and a relative threshold $\mathrm{SU}(d_i)$4, followed by renormalization in the mixed canonical form. The final MPO is embedded into an MPS in $\mathrm{SU}(d_i)$5 and sampled exactly at $\mathrm{SU}(d_i)$6 per draw.

Benchmarks against DMRG

On bimodal Edwards–Anderson spin glasses on hypercubic lattices up to $\mathrm{SU}(d_i)$7 (343 spins), evaluated over five random instances per geometry with exact ground states obtained from CPLEX, the powered-MPO method shows systematic improvement of the approximation ratio (AR) with both iteration count and bond dimension. With a linear schedule at $\mathrm{SU}(d_i)$8, mean peak ARs reach 0.999 on $\mathrm{SU}(d_i)$9, 0.997 on $Z$0, and 0.989 on $Z$1, with exact solutions attained on smaller systems. By contrast, DMRG initialized from the uniform state $Z$2 converges to the same suboptimal solution regardless of whether $Z$3, 64, or 128 — a strong demonstration that variational sweeping stagnates in these rugged landscapes while power iteration does not. Random-initial-state DMRG partially mitigates this, reaching ARs around 0.977–0.997, but still trails the powered-MPO peaks on several geometries, and its improvement with increasing $Z$4 remains marginal.

The two schedules exhibit complementary behavior: doubling converges quickly but saturates or degrades past its peak near $Z$5–11 due to heightened sensitivity of squared-power iteration to truncation error, whereas the linear schedule sustains gradual improvement indefinitely. On the $Z$6 instance, the linear schedule identifies more than 850 unique ground-state bitstrings among 1,000 samples at $Z$7, versus fewer than 200 for doubling — indicating that aggressive squaring collapses sampling diversity prematurely. The observed trade-off motivates hybrid schedules combining rapid initial filtering with stable late-stage refinement.

Higher-order problems on heavy-hexagonal lattices

The method is further tested on HUBO instances with cubic terms on heavy-hexagonal lattices of the type used to compare QAOA against quantum annealing on IBM hardware, with coefficients drawn uniformly from $Z$8. Despite the presence of $Z$9 three-body terms, sparsity keeps the effective MPO bond dimension of $\hat G = -\hat H + \Lambda\,\mathbb I,$0 bounded by roughly $\hat G = -\hat H + \Lambda\,\mathbb I,$1. At $\hat G = -\hat H + \Lambda\,\mathbb I,$2 (about 580 variables) with $\hat G = -\hat H + \Lambda\,\mathbb I,$3, the powered-MPO method already outperforms SA run both natively on the HUBO (via OpenJij) and on an order-reduced QUBO formulation (via dimod with penalty strength 2.0). At $\hat G = -\hat H + \Lambda\,\mathbb I,$4, the linear schedule becomes the only method to surpass the SA solutions, doing so around $\hat G = -\hat H + \Lambda\,\mathbb I,$5. This robustness against the rugged landscapes characteristic of higher-order interactions bears directly on the quantum-computing narrative surrounding these benchmarks: a classical tensor-network baseline achieves solution quality competitive with or superior to SA on exactly the problem class proposed as favorable for QAOA.

The authors also describe a warm-starting heuristic compatible with their framework: applying single-qubit rotations $\hat G = -\hat H + \Lambda\,\mathbb I,$6 to identified configurations generates Hamming-distance-localized superpositions without introducing entanglement, though its quantitative benefit is not evaluated in this work.

Limitations and open questions

The paper concedes several limitations explicitly. Reaching exact solutions requires overcoming substantial entanglement barriers, and the scaling of the required bond dimension across problem classes is left as an open question requiring systematic clarification. Although the linear schedule matches DMRG's per-step cost, it requires considerably more steps to converge, so DMRG remains preferable when walltime or resources are tightly constrained. Extension to two-dimensional tensor networks is explored only preliminarily: projected entangled-pair operator (PEPO) representations were found to have constant bond dimension on square and heavy-hexagonal lattices, but direct PEPO–PEPO contraction proved fragile under belief-propagation-based truncation, leaving more accurate truncation schemes as a necessary development before PEPS-based versions of the method become viable. Finally, comparisons against SA rely on default hyperparameters; the authors acknowledge that tuned or specialized SA variants could narrow the reported gap.

Conclusion

This work reframes combinatorial optimization on classical hardware as spectral amplification of a positive semidefinite MPO, with sampling performed on the induced MPS. Its principal empirical findings are that power iteration systematically improves with bond dimension where DMRG stagnates, and that it outperforms simulated annealing on higher-order Ising instances designed to stress quantum approaches. The theoretical bounds on MPO bond dimensions and iteration counts provide a concrete complexity picture, while the entanglement-growth barrier under truncation defines the central unresolved issue for scaling the method to larger, more connected instances. Beyond its value as a classical solver, the approach establishes a demanding baseline against which claims of quantum advantage on COPs must be measured.