Combining Qiskit Runtime and MATLAB with new Quantum Support Package

Published in

Qiskit

11 min readAug 23, 2023

by Syed Farhan Ahmad, Hamed Mohammadbagherpoor, Francisco Martin Fernandez and Nate-Earnest Noble

MATLAB is a programming and numeric computing platform used by millions of engineers and scientists to analyze data, develop algorithms, and create models. And now, thanks to the new MATLAB Quantum Support Package, MATLAB users can use Qiskit Runtime to build, simulate, and run quantum algorithms on real quantum hardware — all without leaving the MATLAB platform.

The new integration is an exciting example of how Qiskit Runtime provides developers with seamless access to the full functionality of Qiskit, regardless of the programming language they use. In this blog, we highlight the key features of IBM’s Qiskit Runtime service and show how users can integrate it with the MATLAB Quantum Support Package.

What is Qiskit Runtime?

Qiskit Runtime is a quantum computing service and programming model that allows users to optimize workloads and efficiently execute them on quantum systems at scale. The Qiskit programming model was recently extended to allow building algorithms with primitives such as those implemented by IBM’s Qiskit Runtime.

Primitives provide a simplified interface for algorithms and applications. In particular, Qiskit Runtime now has two more abstract primitives: the Estimator and Sampler. They perform foundational quantum computing tasks and act as an entry point to IBM’s Qiskit Runtime service.

IBM Quantum provides accelerated execution of Primitive queries within a runtime. Qiskit Runtime also supports built-in error suppression and error mitigation features that can be configured by the user.

For reference, here is an example of Python code that submits an OpenQASM 3.0 string generated from a Bell State circuit to Qiskit Runtime using REST APIs:

from qiskit import QuantumCircuit, qasm3

# Creating a Quantum Circuit for Bell State
qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)
qc.measure_all()
qc.draw('mpl')

# Converting the qiskit quantum circuit to an OpenQASM 3.0 string
qasm_string = qasm3.dumps(qc)

# Submitting jobs to IBM Quantum backends using the Qiskit Runtime API
url = 'https://runtime-us-east.quantum-computing.ibm.com/jobs'
headers = {'Content-Type': 'application/json',
            'Authorization': 'Bearer '+'<IBM Quantum API Token>'}
job_input = {
    'program_id': 'sampler',
    "backend": 'ibm_lagos',
    "hub": "ibm-q",
    "group": "open",
    "project": "main",
    "params": {
       "circuits": qasm_string, 
       "parameters": [
              [],
          ],
        "circuit_indices": [
            0,
        ],
        "parameter_values": [
            [],   
        ]
    }
}

A Quantum Development Ecosystem

Establishing a quantum development ecosystem that is independent of any specific programming language can foster collaboration among industries, academia and governments along several application verticals. Integrations like the new MATLAB Quantum Support Package help lower the barrier to entry for those looking to start experimenting with quantum computing. In part, that’s because they enable access to quantum hardware in the cloud from a wider array of devices and applications — e.g., Android, industrial IoT appliances, smart healthcare devices, etc.

The Qiskit Runtime service bridges the gap between several programming languages and quantum computing frameworks because it is language-agnostic. It can accept OpenQASM 3.0 strings when submitting jobs to the quantum backend via REST API calls. As long as a programming language supports making REST API calls — which most languages do — and the associated quantum framework supports the generation of OpenQASM 3.0 strings from quantum circuits, users should be able to submit quantum jobs to IBM Quantum compute resources and retrieve results in a language-agnostic manner.

OpenQASM 3.0

OpenQASM is an imperative programming language designed for near-term quantum computing algorithms and applications. Quantum programs are described using the extended quantum circuit model with support for classical feed-forward flow control based on measurement outcomes. It forms a bridge between several quantum programming languages, giving us a universal set of instructions that can run on near-term quantum hardware.

OpenQASM 3.0 extends the quantum assembly language to include timing, pulse control and gate modifiers. There are two different timescales of quantum-classical interactions when it comes to the quantum devices we have today: real-time classical computations that must occur within the coherence time of the qubits; and near-time computations that have more flexible, less stringent timing. Dynamic circuits are based on real-time interactions where classical results provide crucial feedback for quantum operations as seen in quantum error correction routines.

The following sections will focus on how IBM’s Qiskit Runtime can seamlessly integrate with MATLAB with the help of OpenQASM and REST APIs, and no Python involved whatsoever.

MATLAB & the Quantum Support Package

The MATLAB Quantum Support Package lets users build, simulate, and run quantum algorithms. The support package enables you to prototype algorithms to accelerate solving problems such as optimization, scenario simulation, and artificial intelligence (AI) and machine learning tasks, and even invites you to try your hand at tackling currently unsolvable problems in chemistry and material simulations.

With the MATLAB Quantum Support Package, you can:

Build circuits to implement quantum algorithms using a variety of built-in and customizable composite gates
Verify implementation of algorithms using simulations on your local computer or a remote simulator accessed via cloud services
Execute Variational Quantum Algorithms (VQAs) with the help of variational quantum circuits

Qiskit Runtime & The MATLAB Quantum Support Package

With the advent of the Qiskit Runtime cloud service, access to IBM’s quantum hardware is no longer exclusive to Python. The Qiskit Runtime API sends quantum jobs to IBM Quantum compute resources as QASM strings. The MATLAB Quantum Support Package can generate QASM Strings from the quantum circuits that are natively created in MATLAB.

Below, we demonstrate how to encode a classical combinatorial optimization problem into a quantum optimization problem and solve it using the Variational Quantum Eigensolver (VQE) algorithm.

The Test Case: max-cut

Optimization problems appear frequently throughout applied mathematics and computer science, and require finding some optimal solution by minimizing or maximizing a cost function. Max-cut or Maximum Cut is an NP-Hard problem with no polynomial-time algorithm that can solve it.

For a graph, a maximum cut is a cut whose size is at least the size of any other cut. That is, it is a partition of the graph’s vertices into two complementary sets, such that the number of edges between the two sets is as large as possible. Finding such a cut is known as the max-cut problem. In a weighted max-cut problem, each edge is associated with its weight, where the objective is to maximize the total weight of edges between the two complementary sets [2].

Lets create a weighted graph in MATLAB using graph(s, t, weights) and plot it.

%% Generate a graph
s = [1 1 2 3 3 4];
t = [2 5 3 4 5 5];
weights = [1 1 1 1 1 1];
G = graph(s,t,weights);

%% Plot the graph
figure;
h = plot(G);
highlight(h,1:numnodes(G),'MarkerSize',20)

The weighted graph generated by MATLAB can be seen below:

Formal definition of max-cut problem

Consider an n-node undirected graph G = (V, E) where |V| = n with edge. A cut is defined as a partition of the original set V into two subsets. The cost function to be optimized is in this case the sum of weights of edges connecting points in the two different subsets, crossing the cut. The cost to maximize the global profit function is given by:

We can formulate and solve this problem classically using the MATLAB Optimization Toolbox. We begin by formulating the optimization problem and creating the objective function based on the graph’s information.

%% Solve the Maxcut problem Classically
function [sol,fval] =  classical_optimizer(G)
    %%%%% Number of Nodes in the graph
    N = numnodes(G);
    %%%%% Adjacency matrix
    W = full(adjacency(G,'weighted'));
    
    %%% Define binary variables 
    x = optimvar("x",N,1,Type="integer",LowerBound=0,UpperBound=1);
    prob = optimproblem;
    %%% Create the objective function based on the graph info
    for i=1:length(x)
        for j=1:length(x)
            T(i,j) = W(i,j)*x(i)*(1-x(j));
        end
    end
    
    OBJ = sum(sum(T));
    %%% Solve the problem classically
    prob.Objective = -OBJ;
    rng default % For reproducibility
    options = optimoptions('ga','Display','off');
    [sol,fval] = solve(prob,'Options',options); 

end

We then solve the problem classically, like so:

%% Solve the MaxCut problem Classically
[sol,fval] = classical_optimizer(G);

And we can obtain the results as seen below:

disp('Optimal Value:')
disp(sol.x)

disp('MaxCuts:')
disp(fval)

Optimal bit string: 01001

Optimal max-cut value: 5

Solving max-cut Problem in MATLAB using IBM Quantum Compute Resources

We can also solve the max-cut problem quantumly using the Variational Quantum Eigensolver (VQE) algorithm, which involves mapping the problem to an Ising Hamiltonian. We can do this with the following assignment:

where Zi is the Pauli Z operator that has eigenvalues +1 or -1. Doing this, we find that

In other terms, the weighted max-cut problem is equivalent to minimizing the Ising Hamiltonian

Initialize Runtime service and select backend:

Before we begin solving the max-cut problem in MATLAB, let’s connect to the IBM Quantum Compute Resources using an IBM Quantum API Token that can be obtained from the IBM Quantum Platform. We then initiate a Session and enable the Estimator primitive.

apiToken= "Insert your IBM Quantum API Token here";

service = QiskitRuntimeService(apiToken);
service.program_id = "estimator";
service.Start_session = true;
backend="ibm_sherbrooke"; 

%% Initiate the Session and Enable the Estimator Primitive
session = Session(service, backend);  
sampler = Estimator(session=session);

Problem Hamiltonian and Variational Quantum circuit (Ansatz):

To utilize the VQE algorithm for a max-cut problem we require a Pauli Hamiltonian that encodes the cost in a manner such that the minimum expectation value of the operator corresponds to the maximum number of edges between the nodes in two different groups.

For this simple example, the operator is a linear combination of terms with Z operators on nodes connected by an edge.

%%% Convert the Maxcut problem into an Ising Hamiltonian
[hamiltonian.Pauli_Term,hamiltonian.Coeffs,Offset_value, Offset_string] = Maxcut.ToIsing (G);

The output would be:

hamiltonian.Pauli_Term
"IIIZZ"    "IIZZI"    "IZZII"    "ZIIIZ"    "ZIZII"    "ZZIII"    "IIIII"

hamiltonian.Coeffs
 "0.5"    "0.5"    "0.5"    "0.5"    "0.5"    "0.5"    "-3"

Once the operator is constructed, we use MATLAB Quantum Support Package to create the ansatz circuit. TwoLocal ansaztes are created by specifying the number of qubits, repetitions/depth, number of parameters and the entanglement type.

function circuit = Twolocal(circuitinfo,parameters)
    Gates = [];
    for i =1:circuitinfo.reps
        if circuitinfo.entanglement=="linear"
            gates = [
                 ryGate([1:circuitinfo.number_qubits], parameters((i-1)*circuitinfo.number_qubits+1:i*circuitinfo.number_qubits))
                 cxGate(1:circuitinfo.number_qubits-1, 2:circuitinfo.number_qubits)
                    ];
            Gates = [Gates,gates'];
        end
    end
    final_rot = [ryGate([1:circuitinfo.number_qubits], parameters(circuitinfo.reps*circuitinfo.number_qubits+1:(circuitinfo.reps+1)*circuitinfo.number_qubits))]';
    Gates  = [Gates, final_rot];
    circuit = quantumCircuit(Gates);
end

We begin by choosing a TwoLocal variational circuit as your ansatz. In this case, a linear entanglement with 4 reps are used.

circuit.reps=4;
circuit.entanglement = "linear";
circuit.number_qubits = numnodes(G);
circuit.num_parameters = (circuit.reps+1)*circuit.number_qubits;

% Creating the ansatz based on circuit definition and associated angles
ansatz = Twolocal(circuit, angles);
plot(ansatz)

Variational circuit with optimized parameters (angles)

Define the cost function by using Estimator:

As with an iterative optimization procedure, we now need to define our cost function over which to minimize. We define the cost function as follows, computing the expectation value of our Hamiltonian with respect to the parameterized ansatz circuit using the Qiskit Runtime Estimator primitive:

%% Define the cost function to calculate the expectation value of the derived Hamiltonian
function [energy] = cost_function(parameters,arg)    

    % Construct the variational circuit 
    ansatz = Twolocal(arg.circuit, parameters);

    estimator = arg.estimator;
    job       = estimator.run(ansatz,arg.estimator.options.service,arg.hamiltonian);

    %%%% Retrieve the results back
    results = Job.retrieveResults(job.id,arg.estimator.options.service.Access_API);
    energy  = results.values;
end

We need to set the arguments such as maximum iteration, initial values for the circuit parameters, etc. for the MATLAB optimizer.

Minimize the cost function

Any classical optimizer can be used to minimize the cost function. On a real quantum system, an optimizer designed for non-smooth cost function landscapes usually does better. Here we use the surrogate global optimizer.

Because we are iteratively executing many calls to Runtime, we should make use of a Session in order to execute all calls within a single block. Moreover, for the max-cut problem, the solution is encoded in the output distribution of the ansatz circuit bound with the optimal parameters from the minimization. Therefore, we will need a Sampler primitive as well, and will instantiate it with the same Session after the optimizer converges.

%% Arguments for the optimizer 
arg.hamiltonian = hamiltonian;
arg.circuit     = circuit;
arg.estimator   = estimator;

%% Define the cost function
cost_func = @(theta) cost_function(theta,arg);

%% Set the parameters for the MATLAB surrogateopt global optimizer
x0 = -5*ones(circuit.num_parameters,1);
max_iter = 80;

lower_bound = repmat(-2*pi,circuit.num_parameters,1);
upper_bound = repmat( 2*pi,circuit.num_parameters,1);

options = optimoptions("surrogateopt",...
    "MaxFunctionEvaluations",max_iter, ...
    "PlotFcn","optimplotfval",...
    "InitialPoints",x0);

rng default %% For reproducibility

[angles,minEnergy] = surrogateopt(cost_func,lower_bound,upper_bound,options);

Find the solution to the problem using the optimized angles:

The optimized angles are then used to construct the final ansatz circuit and run on the actual hardware which gives us the solution to the original problem.

The job is executed on IBM Quantum compute resources using the Sampler primitive by passing in the ansatz to sampler.run(). The results are retrieved using Job.retrieveResults() which takes in the job ID.

ansatz = Twolocal(circuit, angles);

sampler = Sampler (session=session);
job = sampler.run(ansatz,sampler.options.service);
results = Job.retrieveResults(job.id,sampler.options.service.Access_API);

From the retrieved results, we then identify the solution which is the bitstring that gives the highest probability amongst all other bitstrings. Results are extracted from results.quasi_dists.

%%% Extract the Bitstring
string_data = string(fieldnames(results.quasi_dists));
bitstring_data = replace(string_data,"x","");
%%% Extract the probabilitites
probabilities = cell2mat(struct2cell(results.quasi_dists));

Finally, the results are plotted by invoking theplot_results() function which takes in the Graph, bitstring_data and probabilities obtained above.

% Plot the results and color the graph using the received bit-string

Maxcut.plot_results(G,bitstring_data,probabilities);

The function values, as seen in the Optimization Plot function, reach an optimal value of -4.9375 after 80 iterations.

The returned distribution from Qiskit Runtime Sampler Primitive shows the highest probability for the combination 01101 (recall that the 0th qubit is farthest right) as matches in the solution obtained classically for the max-cut problem.

Distribution obtained after running the final circuit with the optimized parameters (angles)

The final graph is shown below with the maximum cuts.

Conclusion

This blog post demonstrates the integration of IBM Qiskit Runtime with the MATLAB Quantum Support Package. An example of the max-cut problem was solved using the MATLAB Optimization Toolbox and Quantum Support Package, and the circuits were executed on IBM Quantum backends using the Qiskit Runtime Sampler primitive.

The code for the integrations can be found here [3]. We hope you’ll take a look and experiment with it in your work!