The Quantum Mechanics of Discrete Optimization I

In a previous post, we discussed how to turn a classical optimization problem, such as a MILP, into an Ising model or QUBO. I told you it was because these problems are amenable to quantum computers, but how exactly to go from these quadratic models to the actual quantum systems remained a black box. To unveil this mystery, I plan to write a post series going into the nitty-gritty details of converting a classical discrete optimization problem into a quantum system. As a bonus, it will also act as a “theoretical manual” to the inner workings of TenSolver.jl. No physics background is required, although some intuition or standard terminology may be cast here or there. On the other hand, we’ll be using a hell of a lot of linear algebra, so I am assuming you, dear reader, know your tensors well.

This is the first post, aimed at teaching the necessary concepts and representations. Actual algorithms and performance improvements will come later in the series. Since the theory allows it, instead representing just the Ising model, we do the more general Polynomial Unconstrained Integer Optimization \begin{equation} \tag{PUIO} \begin{array}{rl} \min_{x} & \sum_{\alpha} c_\alpha x_1^{\alpha_1}\cdots x_N^{\alpha_N} \\ \text{s.t.} & x_i \in U_i \end{array} \end{equation} where each U_i is a finite set of real numbers, called the domain of variable x_i.¹

We recover the Ising model from this specification by forcing p to be quadratic, i.e., p(x) = x^\top J x + h x and setting all domains to U_i = \{-1, +1\}. Similarly, this model is a QUBO whenever p(x) = x^\top Q x and all domains are U_i = \{0, 1\}. Although we are doing the general case, these binary models will often pop up for intuition throughout the post.

The whole “Quantum for Optimization” approach consists of converting the discrete problem into a Hermitian operator H whose least eigenvalue equals the original optimal, i.e., \begin{array}{rl} \min_{x} & p(x) \\ \text{s.t.} & x_i \in U_i \end{array} = \begin{array}{rl} \min_{\lambda, \psi} & \lambda \\ \text{s.t.} & H\ket\psi = \lambda \ket \psi \end{array} with the corresponding eigenvectors \ket\psi acting as probability distributions over the optimal solutions. These operators are exactly what a quantum computer expects to get, but they also have additional interest because you can use classical numerical linear algebra to grind through them on GPUs.

Quantum System Reformulation

So, how do you turn a finite nonlinear model into a continuous linear one? The magic answer is with probabilities! The central idea is that instead of working with “values”, we should look at probability distributions over them.²

Variables, Domain, Probabilities and States

Let’s say you throw a coin that can land either up or down (heads or tails). This determines a discrete system with possible outcomes \{\downarrow, \uparrow\}. Still, the probability distributions over all such outcomes are infinite and connected. There are the “certain” distributions \ket\downarrow and \ket\uparrow meaning an event absolutely occurs with no chance left for anything else, but there are also the other (more interesting) distributions having the form³ \ket p = (1 - p) \ket{\downarrow} + p \ket{\uparrow}, \quad p \in [0, 1]. Geometrically, this extends the system from just the points \{\ket \downarrow, \ket \uparrow\} to the whole orange line in the figure below.

For the general case, one could go from a finite set U to the simplex of probabilities \sum_{i \in U} p_i \ket i, where p_i \ge 0 and \sum p_i = 1… And this is almost this post’s strategy! Unfortunately, the simplex can be a strange beast to handle algorithm-wise — too many sharp corners and too strange symmetries. So the plan is to follow the physicist’s lead and do the exact same construction as above but using the sphere — the complex sphere — instead. You can’t go smoother than that!

Consider a finite set U \subseteq \mathbb{R} containing d possible outcomes. A state or wavefunction over it is a complex unit vector \ket\psi = \psi_1 \ket 1 + \ldots + \psi_{d} \ket {d}, where each \psi_i \in \mathbb{C} and \sum \left|\psi\right|^2 = 1. Again we write \ket i for the “certain” distribution for the i-th outcome — also called the computational basis. States are sometimes called 2-norm distributions because, one can recover probabilities from a state by noticing that the squared components \left|\psi_i\right|^2 are real, nonnegative and sum to one. This way, each state still provides a clear probabilistic interpretation over U.

To illustrate the geometric difference, we showcase (with real coefficients) a state representing a coin flip with probability p of landing upwards. Notice that classical distributions and quantum states only intersect at the computational basis \ket \downarrow and \ket \uparrow.

Alright, back to optimization. In the original problem (PUIO), the decision runs through a finite domain U, while, in the quantum setting, this becomes the unit sphere⁴ in the finite-dimensional Hilbert space \mathcal{H}\simeq \mathbb{C}^{U}. With this, we start to build a dictionary between discrete models and quantum systems (or linear algebra). To aid in comprehension, let’s add three columns to it. Original stands for the concept in the optimization problem, Classical for the lifting of these concepts to this new setting (but without any new information), and Quantum for the whole new idea including those without a discrete analogous.

Chaining Multiple Variables

Up to now, we turned a single variable into a state represented by a complex vector. But in a real-world optimization problem there are many variables. How does a state model this?

First, let’s look at our coin flipping example. If you flip three coins, what are the possible outcomes? Well, you can enumerate all 8 = 2^3 possible outcomes as \mathcal{C}^3 = \begin{Bmatrix} (\downarrow, \downarrow, \downarrow) & (\downarrow, \downarrow, \uparrow) & (\downarrow, \uparrow, \downarrow) & (\downarrow, \uparrow, \uparrow) \\ (\uparrow, \downarrow, \downarrow) & (\uparrow, \downarrow, \uparrow) & (\uparrow, \uparrow, \downarrow) & (\uparrow, \uparrow, \uparrow), \end{Bmatrix} Or you can follow the smarter route and notice that since each coin can be either \downarrow or \uparrow, their joint outcome is the Cartesian product (conjunction) of the individual outcomes: \mathcal{C}^3 = \{ \downarrow, \uparrow\} \times \{ \downarrow, \uparrow\} \times \{ \downarrow, \uparrow\}.

This chaining structure lifts up to probabilities — be them real or complex. The only difference is that conjunction changes from Cartesian to tensor products,⁵ thanks to the identity \mathbb{C}^{A \times B} \simeq \mathbb{C}^A \otimes \mathbb{C}^B.

Thus, since for each (PUIO) decision variable x_i \in U_i we associate a state space \mathcal{H}_i \simeq \mathbb{C}^{U_i}, the states over the full possible decisions x live in the equivalent tensor product, \begin{gathered} x \in U = U_1 \times \ldots \times U_N \\[0.7em] {\Large\downarrow} \\[0.7em] \ket\psi \in \mathcal{H}= \mathcal{H}_1 \otimes \ldots \otimes \mathcal{H}_N \end{gathered} Notice that if each U_i has d_i elements, the dimension of \mathcal{H} is \prod_{i = 1}^N d_i. This further increases the dictionary.

The quantum world also introduces notions without an equivalent in the discrete model. Variable independence, for example. A state \ket \psi \in \mathbb{C}^A \otimes \mathbb{C}^B defines independent outcomes over A and B when it is a tensor product, \ket\psi = \ket{\psi_A} \otimes \ket{\psi_B}. By checking the dimensions, independent outcomes have only \mathrm{card}(A) + \mathrm{card}(B) degrees of freedom, while the full tensor space has \mathrm{card}(A) \times \mathrm{card}(B) to account for correlations. The existence of those states is the famous phenomenon of entanglement, which despite being irrelevant for this post, will have important comebacks throughout the series.

Minimization via Eigenvalues

Now that you know how to quantize domains, it is time to handle objective functions. Mimicking the previous section, we begin with an “opaque” domain U before decomposing it using tensor products.

Recall that any state \ket\psi defines a classical distribution where outcome x has probability \left|\psi_x\right|^2. By interpreting a finite function f \colon U \to \mathbb{R} as a random variable, a natural thing to do is to use states to calculate expected values, \mathbb{E}^\psi\left[{f}\right] = \sum_{x \in U} \left|\psi_x\right|^2 f(x). In particular, we recover the evaluations f(x) themselves by taking the average with a basis state \ket x. Now let’s proceed to turn the expression above more “linear algebraic”.

A finite function f \colon U \to \mathbb{R} is essentially a vector carrying \mathrm{card}(U) real numbers. And that is indeed how you’d view it if we were doing real-valued probability. However we’re living in a quantum world here, so, as you soon shall see, it makes more sense to reshape it into a diagonal matrix F \colon \mathcal{H}\to \mathcal{H}.

In this setting, one recovers the average by taking inner products: \begin{align*} \braket{\psi | F | \psi} &= \sum_{x \in U} \sum_{y \in U} \bar{\psi_x} \psi_y \braket{x | F | y} \\ &= \sum_{x \in U} \left|\psi_x\right|^2 \braket{x | F | x} \\ &= \sum_{x \in U} \left|\psi_x\right|^2 f(x) \\ &= \mathbb{E}^\psi\left[{f}\right]. \end{align*}

Again, it is useful to visualize it using a simple example on the real unit circle. Suppose again that we have two values \{\downarrow,\uparrow\}. For those, the function f basically amounts to two values f(\downarrow) and f(\uparrow). The quantum version \braket{ \psi | F | \psi }, on the other hand, “fills in” all values between those two, producing an actually continuous function. The illustration below showcases this property for a simple example. You can click on the “animate” button to see how the value varies with different states.

The min-max theorem states that the least eigenvalue of an operator such as F equals the minimum average evaluation on the unit sphere. Furthermore, by construction we know that the eigenvectors of F are exactly the basis states, i.e. F\ket{x} = f(x) \ket{x}. This way, a discrete optimization problem becomes the well-known problem of finding the least eigenvalue of an operator — also called its ground state.

\begin{aligned} \lambda_1(F) &= \inf_{\left\lVert\psi\right\rVert = 1} \braket{\psi | F | \psi} \\ &= \min_{x \in U } \braket{x | F | x} \\ &= \min_{x \in U} f(x). \end{aligned}

In theory, this fulfills this post’s promise: we’ve turned a discrete minimization problem into an eigenvalue problem. Nevertheless, the construction is too opaque, since it makes no use of the tensor product or polynomial structure. Also, to be fair, for a large enough domain U, explicitly enumerating all f(x) is computationally intractable. Even worse: if we’re able to do that, we can just keep a counter and find the optimum instead of employing all that quantum machinery. We still have to dig deeper.

Also keep in mind that although all operators coming from finite functions are diagonal, the theory also works for Hermitian operators. In quantum mechanics it is common to encounter those operators representing more complicated interactions and we may take advantage of that in further posts about actually searching for eigenvalues. For today, though, all operators are diagonal.

The dictionary thus increases.

An Interlude on the Ising Model

It is appropriate to start with an example. Recall the coin flipping example with outcome space \{\downarrow, \uparrow\} \times \{\downarrow, \uparrow\} \times \{\downarrow, \uparrow\}. Let’s say you go to a casino and learn the following betting game: the dealer flips 3 coins in sequence, every time 2 sequential coins align (land with the same face), the dealer pays you an amount J, but when they differ you pay the dealer an amount J. This is modeled by an extremely simple Ising Model, where we do an isomorphism \{ \downarrow, \uparrow\} \to \{-1, +1 \} and write the reward as h(s) = J s_1 s_2 + J s_2 s_3. You can see an illustration of such a system with J = 1 in the interactive diagram below. Click on a “coin” to flip it. Hovering a node shows its current state, while hovering an edge shows its contribution. The total calculated outcome is in the lower caption.

In a casino, you will want to use a probabilistic model instead of a discrete one. To keep up with this post’s topic, we again model the coin flips as quantum states c_{\downarrow} \ket\downarrow+ c_{\uparrow} \ket\uparrow. Doing as in the previous section, we write the reward operator by enumerating all solutions, i.e., H = \begin{bmatrix} 2J & & & & & & & \\ & 0 & & & & & & \\ & & -2J & & & & & \\ & & & 0 & & & & \\ & & & & 0 & & & \\ & & & & & -2J & & \\ & & & & & & 0 & \\ & & & & & & & 2J \\ \end{bmatrix} but this method does not take into account any structure and escalates terribly for more coins. A more interesting approach is to break it into small parts and then reassemble as a polynomial over local operators.

The Pauli matrix Z = \begin{bmatrix*}[r] +1 & 0 \\ 0 & -1 \end{bmatrix*} turns a single spin (or coin) into a number as \begin{array}{rcr} Z \ket\uparrow&=& \ket\uparrow, \\ Z \ket\downarrow&=& -\ket\downarrow. \\ \end{array} We can think of Z as a function realizing a “spin” domain for a state. It is customary to apply it as a local operator acting on a single variable, e.g., \begin{aligned} Z_1 &= Z \otimes \mathrm{I}\otimes \mathrm{I}, \\ Z_2 &= \mathrm{I}\otimes Z \otimes \mathrm{I}, \\ Z_3 &= \mathrm{I}\otimes \mathrm{I}\otimes Z. \end{aligned} These matrices allow us to rewrite the Hamiltonian H in a much more compact form as the polynomial “h(Z)”. Algebraically, it is almost the same as the discrete version! H = J Z_1 Z_2 + J Z_2 Z_3.

General Polynomial Reformulation

The previous example already showed all steps we are now taking. The ensuing section is just a more abstract and general form of the same idea.

Before proceeding, let’s simplify some notation by writing multivariate monomials using multi-indices. That is, for \alpha \in \mathbb{N}^N, define x^\alpha \coloneqq x_1^{\alpha_1}\cdots x_N^{\alpha_N}. This way, multivariate polynomials admit a very compact form p(x) = \sum_{\alpha} c_\alpha x^\alpha, reminiscent of the single variable case. You’ll see how this pays off in legibility.

The first important idea is that a finite domain A = \{a_1,\ldots, a_d\} \subset \mathbb{R} can itself be ordered into a function a \colon \{1,\ldots, d\} \to \mathbb{R}, which we then reshape into a diagonal operator D \ket{k} = a_k \ket{k}. These are local operators recovering a possible evaluation from a state.

For the N decision variables in a PUIO, we actually have N domains U_i which give rise to N operators D_i \colon \mathcal{H}_i \to \mathcal{H}_i representing the variables’ outcomes in the objective function. We want those to work in the full tensor product \mathcal{H} though. Thankfully, the multi-index notation extends to these matrices in the straightforward way, allowing us to write the action of multiple domains in the compact form \begin{aligned} D^\alpha &\colon \mathcal{H}\to \mathcal{H}\\ D^\alpha &\coloneqq D_1^{\alpha_1} \otimes \ldots \otimes D_N^{\alpha_N}. \end{aligned}

With that, we have all ingredients to quantize the objective function itself! We essentially just turn a polynomial p(x) = \sum_{\alpha} c_\alpha x^\alpha into the operator for “p(D)": \begin{align} H &\colon \mathcal{H}\to \mathcal{H}\nonumber \\ H &= \sum_{\alpha} c_\alpha D^\alpha \end{align}

This is exactly the same huge diagonal operator we built previously, but without the need of enumerating all solutions, as you can see by sandwiching between two basis states \begin{align*} \braket{y | H | x} &= \sum_{\alpha} c_\alpha \braket{y | D^\alpha | x} \\ &= \sum_{\alpha} c_\alpha x^\alpha \braket{y | x} \\ &= \begin{cases} \sum_{\alpha} c_\alpha x^\alpha, & \text{if}\;\; y=x \\ 0,& \text{if}\;\; y \ne x. \end{cases} \end{align*}

Some visualisation might be helpful here. You can think of D_i as a “machine” transforming states from its corresponding space \mathcal{H}_i into other values in \mathcal{H}_i.

If we denote the tensor product by stacking the wires, all that D^\alpha does is to apply the corresponding D_i to each wire without interactions between different components. It is a pure tensor operator and, thus, preserves independence. For example, consider D^\alpha = D_1 \otimes D_2^2 \otimes \mathrm{I}\otimes D_4^3,

The Hamiltonian H is formed by a superposition of such operators and, even though it is diagonal, it’s no longer a pure tensor. You can view H as a machine that, while processing a state, actually intertwines the information from different channels.

Finally, it may seem that we didn’t get anything new with this construction… and that is exactly the idea! The operator H, called our system’s Hamiltonian, carries the exact same information as the original problem, but most important of all, it is describable using the exact same amount of information as the original problem. This is important for both manipulating the system classically and for the communication between classical and quantum computers.

We end this section with a new dictionary where instead of states, only operators take the role of representing the original problem.

Farewell for Now

This post focused on the translation process from the language of optimization problems to that of quantum systems. It might be non-obvious sometimes but I hope it now makes sense to you.

If you have something like a D-Wave computer at hand, and want to use it to solve a QUBO or Ising model, this is all you need. As you feed your model to the machine you can proudly say that you know what transformations happen to your data in the process.⁶

On the other hand, for simulating these systems on actual classical computers — such as the laptop I’m writing from right now — there is still a long way to go. After all, we managed to turn an N-variable problem into a 2^N \times 2^N matrix — not very reassuring, unless you can access one of those magical devices. There is no reason to fear, nonetheless! In the next post we shall explore how to write a Hamiltonian using tensor networks. These are classical CPU and GPU-friendly representations capable of losslessly reducing the necessary coefficients to a polynomial amount. Moreover, future posts will also break the “unconstrained” restriction and explain what constraints become in the quantum setting. Stay tuned!