# Quick Dirac-3 Programming Tutorial

Entropy quantum computing (EQC) is a unique tool for optimization leveraging exotic effects in open quantum systems, including quantum Zeno effect. Most quantum and quantum-inspired programming models, including quantum annealing models are restricted to qubits. Dirac-3 hybrid quantum optimization machine, 3rd generation of EQC, allows the use of qudits, which are units of quantum information similar to qubits, but each taking more than two possible values. Be sure to read the Qudit Basics to further understand the benefit of high dimensional encoding using qudits.

This tutorial serves as a quick starter guide, offering practical examples to help users learn the solver characteristics, submit problem and interpret the results from Dirac-3 effectively.

## Device - Dirac-3

Dirac-3 solves quadratic Hamiltonians of up to 949 variables. It also solves highly complex problems with three and four-body interaction. For more details, see the user guide.

## Tutorial Structure

The goal of this tutorial is to explain how to work with Dirac-3.

The first step is to turn an objective function into a Hamiltonian. It is important that the Hamiltonian is important in Hermitian form. That means the interaction matrices or tensors must be symmetric upon conjugation or tranposition operations.

For example, if the objective function is

`$E=x_1-x_2+x_1^2+bx_1x_2$`

Then the Hamiltonian in the matrix form shall be (assuming $b$ is real)

```
$H=\begin{bmatrix} 1 & -1 \end{bmatrix}\begin{bmatrix}x_1\\ x_2\end{bmatrix}
+ \begin{bmatrix} x_1 & x_2 \end{bmatrix} \begin{bmatrix} 1 & b / 2 \\ b / 2 & 0 \end{bmatrix}
\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}$
```

for state vector $\begin{bmatrix}x_1 \\ x_2 \end{bmatrix}$. Note that cross-term in the second part of the Hamiltonian is $b/2$. To code the Hamiltonian into Dirac-3, the linear and quadratic parts are defined separately as

`$C=\begin{bmatrix} 1 \\ -1\end{bmatrix}$`

and

`$J=\begin{bmatrix} 1 & b/2 \\ b/2 & 0 \end{bmatrix}.$`

These two matrices are input into the `sample_hamiltonian`

method in the examples below.

In the following , we will use three different examples to illustrate how to prepare and submit problems.

- The first is a simple example to help with understanding the use of the qudits.
- The second example is a polynomial which will be used for demonstrating how different values of the summation constraint can be used to change the solution of the problem.
- The third example is a simple problem which, like the first, is minimized by setting one value to the whole quantity of $R$. The difference is that there are three solutions at optimality. Repeated runs will reveal different solutions.

First, we'll get our Python environment set up. We need `numpy`

, two methods to sample and evaluate quadratic Hamiltonians, respectively, and the connection to a device.

In [1]:

`import numpy as npimport matplotlib.pyplot as pltfrom qci_client import QciClientimport osos.environ["QCI_API_URL"] = api_url = "https://api.qci-prod.com"os.environ["QCI_TOKEN"] = token = "" # REPLACE WITH YOUR TOKENqclient = QciClient(api_token=token, url=api_url)`

The `sample_hamiltonian`

method takes six required arguments. `C`

and `J`

are the linear and quadratic terms of the polynomial. These matrices are dense and the method converts into a sparse format for the `process_job`

call. The `sum_constraint`

parameter is a value which indicates the total returned given solution vector must add up to. The `schedule`

parameter takes four different values, 1, 2, 3 or 4.

- This is the quickest execution. It has the lowest probability of obtaining high quality solutions.
- This execution takes more time, likely a few seconds, but has a higher probability of finding good solutions.
- This option may take tens of seconds. The quality of these solutions is expected to be higher.
- This is the longest running relaxation option. It takes up to multiple minutes to run and has the highest probability of returning quality solutions.

The `solution_precision`

parameter is a value which indicates how granular a solution should be. Use a value of 1 for integer solutions. Use decimals such as 0.1 or 0.05 for higher precision with the highest precision accepted being 0.00001. Specify a connected `EqcClient`

in the `client`

parameter. Suppress print statements in the call using a value of `True`

for `suppress_output`

.

In [2]:

`def sample_hamiltonian(C : np.ndarray, J : np.ndarray, sum_constraint : float, schedule : int, solution_precision : float, client : QciClient):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_precision": solution_precision,"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)# total_device_time += resp["runtime"]return responsedef get_results(response):if "results" in response and response["results"] is not None:results_file_config = response["results"]["file_config"]# results file config has one key, named by the job type detailassert len(results_file_config) == 1, "Unknown results format"results = list(results_file_config.values())[0]else:if "job_info" in response and "job_result" in response["job_info"]:details = response["job_info"]["job_result"]else:details = Noneraise RuntimeError(f"Execution failed. See details: {details}")return results`

## Qudit Domain

Dirac-3 uses qudits to solve discrete optimization problems.

The first example is a three-variable problem with uniform two-body interaction, whose solution corresponds to one variable taking on the total value of the summation constraint and the rest taking the value 0.

In the coding please note the `dtype`

of `float32`

. This is the digital precision for floating point numbers used in the gRPC protocol. A warning will be raised if higher precision decimal values are used.

In [3]:

`h = np.array([[-1],[0],[0]], dtype=np.float32)J= np.array([[0, 1, 1],[1, 0, 1],[1, 1, 0]], dtype=np.float32)h, J`

Out [3]:

(array([[-1.], [ 0.], [ 0.]], dtype=float32), array([[0., 1., 1.], [1., 0., 1.], [1., 1., 0.]], dtype=float32))

In [4]:

`response = sample_hamiltonian(h, J, 400, 1, 0.1, qclient)if "job_info" in response:print("Status:", response["job_info"]["details"]["status"])results = get_results(response)print("Energy:", results["energies"][0])print("Solution")print(results["solutions"][0])solution = np.array(results["solutions"][0])print("Solution Value (should match energy)", h.T@solution + solution.T@J@solution)x = np.array([400, 0, 0])print("Known ground state", h.T@x + x.T@J@x)print("With solution", x)else:print(response)`

Out [ ]:

Dirac allocation balance = 0 Job submitted job_id='65dff2f21e7c62f822901851'-: 2024/02/28 19:58:58 running: 2024/02/28 19:58:59 completed: 2024/02/28 19:59:46 Dirac allocation balance = 0 Status: completed Energy: -400 Solution [400, 0, 0] Solution Value (should match energy) [-400.] Known ground state [-400.] With solution [400 0 0]

## Dynamic Range

Solution precision is dictated by the machine sensitivity and precision, including those of the photon-counting modules and optoelectronic devices, and shot-noise limited quantum measurements. Dirac-3 has a dynamic range of at least 23 dB. This means that Hamiltonian coefficients having an absolute value smaller than ~1/200 of the peak coefficient value may not be recognized and will be effectively treated as 0 when encoding the Hamiltonian. In some cases, the dynamic range can reach as high as 40 dB. However, it is not guaranteed.

During computing cycles, each qudit takes discrete numbers of levels, which is at least 200, but can reach 10,000, as limited by the 23 dB or 40 dB dynamic range. At the output, however, the optimal solutions may consist of variables that have fractional components, as a result of normalization to the summation constraint. Thus when the number of levels is 200, each qudit value will be a multiple of R/200.

Device sensitivity is not a controllable parameter and so the number of levels is not a part of problem formulation. More details on this are available in the Dirac-3 User Guide.

## Summation Constraint

The summation constraint is required to be some value between 1 and 10,000. This restriction is a property of Dirac-3.1 and not EQC in general. This constraint is helpful in some models and a formulation hurdle in others. In all formulations of TSP, for instance, the sum of a final solution is known beforehand. A Hamiltonian cycle (another one of those confusing terminologies that shows up when mathematicians and physicists work together) has exactly $n$ edges, where $n$ is the number of nodes in a graph. This means that solving a TSP with Dirac-3 requires setting the sum constraint equal to the number of nodes in the graph.

Other cases where the sparsity or sum of the solution is not known before solving can be solved with a couple of different approaches.

### Machine Slack Qudits

When formulating a model, it can be advantageous to add slack variables. To use this method when setting a summation constraint, formulate the Hamiltonian like normal, then determine a lower bound and an upper bound on the sum of the qudit values. Introduce an additional qudit into the formulated model which can sum to the upper bound minus the lower bound. This additional variable in the model will be ignored in the solution. See some of the more indepth examples, like Max Cut, to understand slack qudits in action.

### Minimizing a simple polynomial

Let $H(x)=-5x_1x_3-4x_2^2-10x_1-10x_3$. This polynomial itself has no lower bound, but as we put a summation constraint, it does. If $\sum_i x_i = R$, then $H(s)$ has a minimum value which changes negatively with $R$ and changes the cardinality of $x$. Below are the results for different $R$ values.

In [5]:

`h = np.array([[-10.],[ 0.],[-10.]], dtype=np.float32)J = np.array([[0., 0., -2.5],[0., -4., 0. ],[-2.5, 0., 0. ]], dtype=np.float32)for S in [2, 3, 4, 5, 6, 7]:print("**********************************")print(f"S={S}")response = sample_hamiltonian(h, J, S, 3, 0.1, qclient)results = get_results(response)print("Status:", response["job_info"]["details"]["status"])results = get_results(response)print("Energy:", results["energies"][0])print("Solution")print(results["solutions"][0])print()`

Out [ ]:

********************************** S=2 Dirac allocation balance = 0 Job submitted job_id='65dff3241e7c62f822901852'-: 2024/02/28 19:59:48 running: 2024/02/28 19:59:50 completed: 2024/02/28 19:59:56 Dirac allocation balance = 0 Status: completed Energy: -24.99830436706543 Solution [0.9815806150436401, 1.0000000116860974e-07, 1.0184195041656494] ********************************** S=3 Dirac allocation balance = 0 Job submitted job_id='65dff32f1e7c62f822901853'-: 2024/02/28 19:59:59 running: 2024/02/28 20:00:00 completed: 2024/02/28 20:00:06 Dirac allocation balance = 0 Status: completed Energy: -41.234031677246094 Solution [1.5565037727355957, 2.0000000233721948e-07, 1.4434959888458252] ********************************** S=4 Dirac allocation balance = 0 Job submitted job_id='65dff3391e7c62f822901854'-: 2024/02/28 20:00:09 running: 2024/02/28 20:00:10 completed: 2024/02/28 20:00:16 Dirac allocation balance = 0 Status: completed Energy: -59.988162994384766 Solution [2.0486626625061035, 1.0000000116860974e-07, 1.9513375759124756] ********************************** S=5 Dirac allocation balance = 0 Job submitted job_id='65dff3431e7c62f822901855'-: 2024/02/28 20:00:19 running: 2024/02/28 20:00:20 completed: 2024/02/28 20:00:26 Dirac allocation balance = 0 Status: completed Energy: -100 Solution [0, 5, 0] ********************************** S=6 Dirac allocation balance = 0 Job submitted job_id='65dff34d1e7c62f822901856'-: 2024/02/28 20:00:29 running: 2024/02/28 20:00:30 completed: 2024/02/28 20:00:36 Dirac allocation balance = 0 Status: completed Energy: -104.99845123291016 Solution [3.017604351043701, 1.0000000116860974e-07, 2.982395648956299] ********************************** S=7 Dirac allocation balance = 0 Job submitted job_id='65dff3571e7c62f822901857'-: 2024/02/28 20:00:39 running: 2024/02/28 20:00:40 completed: 2024/02/28 20:00:46 Dirac allocation balance = 0 Status: completed Energy: -196 Solution [1.0000000116860974e-07, 7, 0]

For reference, here is a job response that is returned in JSON from the REST API

In [6]:

response

Out [6]:

{'job_info': {'job_id': '65dff3571e7c62f822901857', 'job_submission': {'problem_config': {'normalized_qudit_hamiltonian_optimization_continuous': {'hamiltonian_file_id': '65dff356fd7ff43e52669ab3'}}, 'device_config': {'dirac-3': {'num_samples': 1, 'sum_constraint': 7, 'relaxation_schedule': 3}}}, 'job_status': {'submitted_at_rfc3339nano': '2024-02-29T03:00:39.435Z', 'queued_at_rfc3339nano': '2024-02-29T03:00:39.436Z', 'running_at_rfc3339nano': '2024-02-29T03:00:40.249Z', 'completed_at_rfc3339nano': '2024-02-29T03:00:46.68Z'}, 'job_result': {'file_id': '65dff35efd7ff43e52669ab5', 'device_usage_s': 1}, 'details': {'status': 'completed'}}, 'results': {'file_id': '65dff35efd7ff43e52669ab5', 'num_parts': 1, 'num_bytes': 256, 'file_config': {'normalized_qudit_hamiltonian_optimization_continuous_results': {'counts': [1], 'energies': [-196], 'solutions': [[1.0000000116860974e-07, 7, 0]]}}}}

## Degeneracy Demonstration

Dirac-3 will return various solutions to a problem with multiple solutions near the optimal value.

Let's see what happens with a clearly degenerate problem. Take the first example modified to give the same weight to all the linear terms. Each solution: `[S, 0, 0]`

, `[0, S, 0]`

and `[0, 0, S]`

give the same value for the Hamiltonian.

In [7]:

`h = np.array([[-1],[-1],[-1]], dtype=np.float32)J = np.array([[0, 1, 1],[1, 0, 1],[1, 1, 0]], dtype=np.float32)`

In [8]:

`super_test = {}for i in range(10):response = sample_hamiltonian(h, J, 400, 1, 0.1, qclient)results = get_results(response)solution = [round(v, 1) for v in results["solutions"][0]]solution = tuple(solution)if solution in super_test:super_test[solution] += 1else:super_test[solution] = 1print("Solution | Frequency")for solution in super_test:print(solution, " ", "*"*super_test[solution])labels = list(super_test.keys())sizes = [super_test[key] for key in labels]fig, ax = plt.subplots()patches, labels = ax.pie(sizes, labels=labels)`

Out [ ]:

Dirac allocation balance = 0 Job submitted job_id='65dff3611e7c62f822901858'-: 2024/02/28 20:00:49 completed: 2024/02/28 20:00:50 Dirac allocation balance = 0 Dirac allocation balance = 0 Job submitted job_id='65dff3651e7c62f822901859'-: 2024/02/28 20:00:53 completed: 2024/02/28 20:00:54 Dirac allocation balance = 0 Dirac allocation balance = 0 Job submitted job_id='65dff3691e7c62f82290185a'-: 2024/02/28 20:00:57 completed: 2024/02/28 20:00:59 Dirac allocation balance = 0 Dirac allocation balance = 0 Job submitted job_id='65dff36d1e7c62f82290185b'-: 2024/02/28 20:01:01 completed: 2024/02/28 20:01:03 Dirac allocation balance = 0 Dirac allocation balance = 0 Job submitted job_id='65dff3711e7c62f82290185c'-: 2024/02/28 20:01:05 completed: 2024/02/28 20:01:07 Dirac allocation balance = 0 Dirac allocation balance = 0 Job submitted job_id='65dff3751e7c62f82290185d'-: 2024/02/28 20:01:09 completed: 2024/02/28 20:01:11 Dirac allocation balance = 0 Dirac allocation balance = 0 Job submitted job_id='65dff3791e7c62f82290185e'-: 2024/02/28 20:01:13 completed: 2024/02/28 20:01:15 Dirac allocation balance = 0 Dirac allocation balance = 0 Job submitted job_id='65dff37e1e7c62f82290185f'-: 2024/02/28 20:01:18 running: 2024/02/28 20:01:19 completed: 2024/02/28 20:01:21 Dirac allocation balance = 0 Dirac allocation balance = 0 Job submitted job_id='65dff3831e7c62f822901860'-: 2024/02/28 20:01:23 completed: 2024/02/28 20:01:25 Dirac allocation balance = 0 Dirac allocation balance = 0 Job submitted job_id='65dff3881e7c62f822901861'-: 2024/02/28 20:01:28 running: 2024/02/28 20:01:30 completed: 2024/02/28 20:01:31 Dirac allocation balance = 0 Solution | Frequency (0.0, 400.0, 0.0) **** (400, 0, 0) ******

Out [ ]:

<Figure size 640x480 with 1 Axes>