A gentle introduction to optimization

Device: Dirac-3

Entropy quantum computing (EQC) is an analog quantum computing paradigm for optimisation. Optimisation problems arise in many contexts, such as economics, drug discovery, finance, energy, supply chain, artificial intelligence (AI), transportation, and many more. While many users may be familiar with optimisation, some may not be. The aim of this tutorial is to provide a guide for those who are not familiar with the topic, for a more in-depth discussion of combinatorial optimisation, a particularly important subfield of optimisation, please see this lesson.

options

Our 3rd generation of EQC, Dirac-3, makes use qudits variables. Unlike qubits, qudits can take more than two possible values. This means that Dirac-3 can solve problems beyond binary (0,1), including integers and continuous numbers (all positive real rational numbers). For further information on qudits, please read through the Qudit Primer to better understand the benefits of high-dimensional programming.

This tutorial walks through how to program an EQC device, specifically Dirac-3, explains what an optimization problem is, and demonstrates with examples. For a more focused tutorial on how to start using the device, without introductory concepts, see our Dirac-3 quick start.

What Is Optimization

Optimization starts with identifying goals: searching for the best solution from a range of options, such as maximizing profits, minimizing costs, finding the quickest route, or using resources efficiently. These goals are often referred to as the objective function, and less commonly in more physics oriented literature, energy.

options

Constraints are the rules or limitations, such as budget restrictions, time limits, or fuel constraints. Sometimes constraints are even just statements of physical realities, for example that a travelling salesperson can only be in a single city at a time. They play a crucial role in optimization because they determine the maximum or minimum possibilities.

In a nutshell, Dirac-3 is like having a genius mathematician working for you, ensuring that you make good choices, saving you time, energy, and resources.

A rectangular garden is to be constructed using a rock wall as one side of the garden and wire fencing for the other three sides (Figure 1). Given 100ft of wire fencing, determine the dimensions that would create a garden of maximum area. What is the maximum area?

image.png

First, let 𝑥 denote the length of the side of the garden perpendicular to the rock wall and 𝑦 denote the length of the side parallel to the rock wall. Then the area of the garden is

A=xy.\begin{equation*} A = xy. \end{equation*}
We want to find the maximum possible area subject to the constraint that the total fencing is 100ft. From Figure 1, the total amount of fencing used will be 2𝑥+𝑦. Therefore, the constraint equation is
2x+y=100.\begin{equation*} 2x + y = 100. \end{equation*}
Solving this equation for 𝑦, we have 𝑦 = 100 − 2𝑥. Thus, we can write the area as
A(x)=x(1002x)=100x2x2\begin{equation*} A(x) = x(100 - 2x) = 100x - 2x^2 \end{equation*}

Let's go ahead and solve this optimization problem with Dirac-3.

In [6]:

  • # Import the necessary packages
  • import numpy as np
  • import networkx as nx
  • import matplotlib.pyplot as plt
  • import os
  • from qci_client import QciClient

Connect to Dirac-3 through the cloud using your unique token and access the QCi Client API. (if you don't have a token, go sign up for our Free Trial Cloud Access).

In [ ]:

  • # Input your access token
  • token = "your_token"
  • # Connect to the API URL
  • api_url = "https://api.qci-prod.com"
  • qclient = QciClient(api_token=token, url=api_url)

The objective of EQC devices is to find the global minimum, representing the lowest values. However, in the problem we're solving, we aim to maximize the area with a limited amount of fence. Therefore, we need to negate the equation:

minimize100x2x2\mathrm{minimize}\,\, 100x-2x^{2}

becomes

maximize2x2100x.\mathrm{maximize}\,\, 2x^{2}-100x.

Formulating our equation into Dirac-3

Dirac-3 takes an Ising Hamiltonian (up to fifth-order) as input and determines the lowest energy state. An Ising Hamiltonian represents a mathematical expression used to calculate the optimal energy or cost associated with a particular configuration of variables (e.g., configuring magnets or covering a specific area with limited fencing), aiming to minimize this energy or cost.

Hobjective=i=1NCiVi+i,j=1N,NJijViVj+i,j,k=1N,N,NTijkViVjVk+i,j,k,l=1N,N,N,NQijklViVjVkVl+i,j,k,l,m=1N,N,N,N,NPijklmViVjVkVlVmH_\mathrm{objective} = \sum_{i=1}^{N} C_{i}V_{i} + \sum_{i,j=1}^{N, N} J_{ij} V_{i}V_{j} + \sum_{i,j,k=1}^{N, N, N} T_{ijk} V_{i}V_{j}V_{k} + \sum_{i,j,k,l=1}^{N, N, N, N} Q_{ijkl} V_{i}V_{j}V_{k}V_{l} + \sum_{i,j,k,l,m=1}^{N, N, N, N, N} P_{ijklm} V_{i}V_{j}V_{k}V_{l}V_{m}
Rconstraints=i=1ViR_\mathrm{constraints} = \sum_{i=1} V_{i}

Since our problem only involves quadratic terms, we do not need higher order couplings, the expanded out objective function for our two variable problem (the reason for introducing a second variable will become clear later) becomes

Hobjective=C1V1+C2V2+J11V1V1+J12V1V2+J21V2V1+J22V2V2H_{objective} = C_{1}V_{1} + C_{2}V_{2} + J_{11} V_{1}V_{1} + J_{12} V_{1}V_{2} + J_{21} V_{2}V_{1} + J_{22} V_{2}V_{2}
Rconstraints=100R_\mathrm{constraints} = 100

including terms which evaluate to zero, our objective is

Hobjective=2x12+0x1x2+0x2x1+0x22100x1+0x2H_\mathrm{objective}=2\,x^{2}_{1}+0\,x_{1}x_{2}+0\,x_{2}x_{1}+0x^{2}_{2}-100\,x_{1}+0\,x_{2}

Convert polynomial into matrix form:

[x1x12x1x2x1xnx2x2x1x22x2xnxmxmx1xmx2xmxn]C=[x1x2xm]J=[x12x1x2x1xnx2x1x22x2xnxmx1xmx2xmxn]C=[1000],J=[2000] \begin{bmatrix} x_{1} & x_{1}^{2} & x_{1}x_{2} & \cdots & x_{1}x_{n} \\ x_{2} & x_{2}x_{1} & x^{2}_2 & \cdots & x_{2}x_{n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ x_{m} & x_{m}x_{1} & x_{m}x_{2} & \cdots & x_{m}x_{n} \\ \end{bmatrix} \quad \rightarrow \quad C = \begin{bmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{m} \\ \end{bmatrix} \quad J = \begin{bmatrix} x_{1}^{2} & x_{1}x_{2} & \cdots & x_{1}x_{n} \\ x_{2}x_{1} & x^{2}_2 & \cdots & x_{2}x_{n} \\ \vdots & \vdots & \ddots & \vdots \\ x_{m}x_{1} & x_{m}x_{2} & \cdots & x_{m}x_{n} \\ \end{bmatrix} \quad \rightarrow \quad C = \begin{bmatrix} -100 \\ 0 \\ \end{bmatrix} \quad , \quad J = \begin{bmatrix} 2 &0 \\ 0 &0 \\ \end{bmatrix}

An astute reader might be curious why we included x2x_2 at all, since all the terms involving it our zero, the reason relates to an addtional constraint on the problem information from the way the hardware designed known as a "sum constraint". We won't discuss the details here, but they can be found in our Dirac-3 quick start.

In [190]:

  • C = np.array([[-100],[0]])
  • J = np.array([[2, 0], [0, 0]])
  • C, J
  • # Constraint (set to fence limit here, but generally could be set to any larger value as well due to presence of second variable)
  • sum_constraint = 100

Copy the ising_solver class to incorporate the necessary functions.

In [191]:

  • def ising_solver():
  • def sample_hamiltonian(C : np.ndarray, J : np.ndarray, sum_constraint : float, schedule : int, solution_precision : float, client : QciClient):
  • if int(solution_precision) == solution_precision:
  • soltype = "integer"
  • else:
  • soltype = "continuous"
  • n = C.shape[0]
  • H = np.hstack([C.reshape([n, 1]), J])
  • ham_file = {"file_name": "qudit-tutorial-hame", "file_config": {"hamiltonian": {"data": H}}}
  • file_id = client.upload_file(ham_file)["file_id"]
  • job_tags = ["qudit-tutorial"]
  • job_body = client.build_job_body(job_type="sample-hamiltonian", hamiltonian_file_id=file_id,
  • job_params={"sampler_type": "dirac-3", "nsamples": 1, "solution_type": soltype,
  • "sum_constraint": sum_constraint, "relaxation_schedule": schedule}, job_tags=job_tags)
  • response = client.process_job(job_type="sample-hamiltonian", job_body=job_body, wait=True)
  • return response
  • def get_results(response):
  • results_file_config = response["results"]["file_config"]
  • assert len(results_file_config) == 1, "Unknown results format"
  • return list(results_file_config.values())[0]
  • def show_results():
  • results = get_results(response)
  • print("Status:", response["job_info"]["details"]["status"])
  • results = get_results(response)
  • print("Solution")
  • print(results["solutions"][0])
  • solution = np.array(results["solutions"][0])
  • print("Energy:", results["energies"][0])
  • print("Solution Value (should match energy):", C.T@solution + solution.T@J@solution)

Now, you need to choose how precise you want your solutions to be with the solution_precision. This allows a distillation process which can give approximate solution to integer problems, an alternative is to treat the problem as continuous. In this case the optimal solution to the continious problem has fence lengths which take integer values, so either can be used.

For schedule, Schedules 1, 2, 3, and 4 correspond to different time settings. Higher schedule numbers indicate longer runtime and, consequently, a higher probability of obtaining favorable results.

In [192]:

  • schedule = 1
  • solution_precision = 0.1
  • response = sample_hamiltonian(C, J, sum_constraint, schedule, solution_precision, qclient)
  • get_results(response)
  • show_results()

Out [ ]:

Dirac allocation balance = 0
Job submitted job_id='66354d3fd448b017e54f9436'-: 2024/05/03 16:46:55
running: 2024/05/03 16:46:56
completed: 2024/05/03 16:47:30
Dirac allocation balance = 0
Status: completed
Solution
[25, 75]
Energy: -1250
Solution Value (should match energy): [-1250]

Breakdown

Our solution

x1=25,x2=75x_{1}=25, x_{2}=75

evaluates to

2x12+0x1x2+0x2x1+0x22100x1+0x2=2x^{2}_{1}+0x_{1}x_{2}+0x_{2}x_{1}+0x^{2}_{2}-100x_{1}+0x_{2}=
2x12100x1=2x^{2}_{1}-100x_{1}=
2(25)2100(25)=02*(25)^{2}-100*(25) = 0

25 is the optimal solution for 𝑥, the result of the evaluation is -1,250 corresponding to a maximum area of 1,250. A more step-by-step method is to determine 𝑦 to maximize the area.

2x+y=1002(25)+y=100y=502x+y=100 \rightarrow 2*(25)+y=100 \rightarrow y=50
Amax=2550A_{max} = 25*50
Amax=1250.A_{max} = 1250.

This solution for maximum area aligns with the Energy and Solution Value, which are negative because we negated our objective function. The terminology of calling the objective function value "energy" orginates in physics, where optimisation processes can be thought of as analogies to cooling processes which minimise system energy.

Conclusion

This tutorial has been a very gentle introduction into optimization specifically integer/continuous optimization with Dirac-3. If you want to start solving your own problems with this device then please see Dirac-3 quick start. If you want more details about the device itself, please see the Dirac-3 product page. For a much more in-depth understanding, please read our recent preprint about the device. For examples of different kinds of optimisation please try our tutorials on quadratic unconstrianed binary optmization and quadratic linearly constrained binary optimization.