Quantum Simulation Modeling (QuaSiMo)

Table of Contents ¶

• Overview
• Simulation Model
• Workflow
• Cost Function Evaluate
• Implement a new workflow

Overview ¶

The QuaSiMo library provides domain-specific tools for quantum simulation on quantum computers. It supports problems such as ground-state energy computations or time-dependent simulations.

The library comprises the main drivers, so-called workflows (QuantumSimulationWorkflow), which encapsulate the procedure (classical and quantum routines) to solve the quantum simulation problem. The input to QuaSiMo’s workflows is a QuantumSimulationModel , specifying all the parameters of the quantum simulation problem, e.g., the observable operator, the Hamiltonian operator, which may be different from the observable or is time-dependent, etc.

The QuaSiMo library provides a ModelFactory factory to facilitate QuantumSimulationModel creation for common use cases.

Built-in workflows can be retrieved from the registry by using getWorkflow helper function with the corresponding workflow name, such as vqe, qaoa, qite, etc.

Some workflow may require additional configurations, which should be provided when calling getWorkflow. Please refer to specific workflow sections for information about their supported configurations.

The retrieved workflow instance can then be used to solve the QuantumSimulationModel problem using the execute method, which returns the result information specific to that workflow, such as the ground-state energy for variational quantum eigensolver workflow or the time-series expectation values for time-dependent quantum simulation.

For advanced users or workflow developers, the QuaSiMo library also provides interfaces for state-preparation (ansatz) circuit generation and cost function (observable) evaluator.

The first can be used during workflow execution to construct the quantum circuits for evaluation, e.g., variational circuits of a particular structure or first-order Trotter circuits for Hamiltonian evolution.

The latter provides an abstraction for quantum backend execution and post-processing actions to compute the expectation value of an observable operator. For example, the cost function evaluator may add necessary gates to change the basis according to the observable operators, analyze the bitstring result to extract the expectation value. For more information, please refer to the Cost Function Evaluate section .

Quantum Simulation Model ¶

The input to all QuaSiMo workflows is a QuantumSimulationModel, which can be generated by the createModel method from the ModelFactory factory. There are a couple of createModel overloads.

At the minimum, the problem model contains information about the target observable whose expectation value is the optimization objective of the workflow. The observable is constructed as a QCOR operator . For example,

• In C++:

// Create the Deuteron Hamiltonian
auto H = 5.907 - 2.1433 * X(0) * X(1) - 2.143 * Y(0) * Y(1) + 0.21829 * Z(0) -
          6.125 * Z(1);
auto problemModel = QuaSiMo::ModelFactory::createModel(H);

• In Python:

# Define the deuteron hamiltonian 
H = -2.1433 * X(0) * X(1) - 2.1433 * \
    Y(0) * Y(1) + .21829 * Z(0) - 6.125 * Z(1) + 5.907

problemModel = QuaSiMo.ModelFactory.createModel(H)

If a specific state-preparation ansatz circuit is to be used, a QCOR kernel function can also be provided to the createModel method as follows,

# Define the deuteron hamiltonian 
H = -2.1433 * X(0) * X(1) - 2.1433 * \
    Y(0) * Y(1) + .21829 * Z(0) - 6.125 * Z(1) + 5.907

# Define the quantum kernel
@qjit
def ansatz(q : qreg, theta : float):
    X(q[0])
    Ry(q[1], theta)
    CX(q[1], q[0])

# Create the problem model, provide the state 
# prep circuit, Hamiltonian and note how many qubits
# and variational parameters 
num_params = 1
problemModel = QuaSiMo.ModelFactory.createModel(ansatz, H, num_params)

The above example is for Python. C++ API is completely identical and is shown in the VQE section below.

When the Hamiltonian operator is different from the observable operator, e.g. for systems with time-dependent Hamiltonian, users can provide this Hamiltonian as a function from time value to concrete operator value. Noted that a static Hamiltonian can be returned if the system is time-independent. Please refer to the time-dependent simulation section for examples.

Workflow ¶

The QuaSiMo library has built-in implementations for Variational Quantum Eigensolver algorithm (VQE), Quantum Approximate Optimization Algorithm (QAOA), Quantum Imaginary Time Evolution (QITE), and general time-dependent evolution workflows. Here, we show simple examples for common use cases and explain extra configuration parameters that each workflow may take.

Variational Quantum Eigensolver - VQE ¶

A Variational Quantum Eigensolver algorithm workflow instance can be retrieved from the QuaSiMo registry by calling getWorkflow("vqe") as shown in the below example.

#include "qcor_qsim.hpp"

// Define a fixed ansatz as a QCOR kernel
__qpu__ void ansatz(qreg q, double theta) {
  X(q[0]);
  auto exponent_op = X(0) * Y(1) - Y(0) * X(1);
  exp_i_theta(q, theta, exponent_op);
}

int main(int argc, char **argv) {
  // Create the Deuteron Hamiltonian
  auto H = 5.907 - 2.1433 * X(0) * X(1) - 2.143 * Y(0) * Y(1) + 0.21829 * Z(0) -
           6.125 * Z(1);
  const auto num_qubits = 2;
  const auto num_params = 1;
  auto problemModel =
      QuaSiMo::ModelFactory::createModel(ansatz, H, num_qubits, num_params);
  auto optimizer = createOptimizer("nlopt");
  // Instantiate a VQE workflow with the nlopt optimizer
  auto workflow = QuaSiMo::getWorkflow("vqe", {{"optimizer", optimizer}});

  // Result should contain the ground-state energy along with the optimal
  // parameters.
  auto result = workflow->execute(problemModel);

  const auto energy = result.get<double>("energy");
  std::cout << "Ground-state energy = " << energy << "\n";
  return 0;
}

from qcor import * 

# Define the deuteron hamiltonian 
H = -2.1433 * X(0) * X(1) - 2.1433 * \
    Y(0) * Y(1) + .21829 * Z(0) - 6.125 * Z(1) + 5.907

# Define the quantum kernel by providing a 
# python function that is annotated with qjit for 
# quantum just in time compilation
@qjit
def ansatz(q : qreg, theta : float):
    X(q[0])
    Ry(q[1], theta)
    CX(q[1], q[0])

# Create the problem model, provide the state 
# prep circuit, Hamiltonian and note how many qubits
# and variational parameters 
num_params = 1
problemModel = QuaSiMo.ModelFactory.createModel(ansatz, H, num_params)
      
# Create the NLOpt derivative free optimizer
optimizer = createOptimizer('nlopt')

# Create the VQE workflow
workflow = QuaSiMo.getWorkflow('vqe', {'optimizer': optimizer})

# Execute and print the result
result = workflow.execute(problemModel)
energy = result['energy']
print(energy)

As we may see, users can provide a specific optimizer instance to the VQE workflow to be used. By default, the nlopt(gradient-free) optimizer will be used if none is provided.

VQE workflow returns the following information:

Quantum Approximate Optimization Algorithm - QAOA ¶

The Quantum Approximate Optimization Algorithm (QAOA) workflow has the following configuration parameters:

A simple QAOA workflow execution is shown in the below example.

#include "qcor_qsim.hpp"

int main(int argc, char **argv) {
  // Create the Deuteron Hamiltonian
  auto H = 5.907 - 2.1433 * X(0) * X(1) - 2.143 * Y(0) * Y(1) + 0.21829 * Z(0) -
           6.125 * Z(1);
  auto problemModel = QuaSiMo::ModelFactory::createModel(H);
  auto optimizer = createOptimizer("nlopt");
  // Instantiate a QAOA workflow with the nlopt optimizer
  // "steps" = the (p) param in QAOA algorithm.
  auto workflow = QuaSiMo::getWorkflow("qaoa", {{"optimizer", optimizer}, {"steps", 8}});

  // Result should contain the ground-state energy along with the optimal
  // parameters.
  auto result = workflow->execute(problemModel);

  const auto energy = result.get<double>("energy");
  std::cout << "Ground-state energy = " << energy << "\n";
  return 0;
}

from qcor import * 

# Solving deuteron problem using QuaSiMo QAOA

# Define the deuteron hamiltonian 
H = -2.1433 * X(0) * X(1) - 2.1433 * \
    Y(0) * Y(1) + .21829 * Z(0) - 6.125 * Z(1) + 5.907

problemModel = QuaSiMo.ModelFactory.createModel(H)
      
# Create the NLOpt derivative free optimizer
optimizer = createOptimizer('nlopt')

# Create the QAOA workflow with p = 8 steps
workflow = QuaSiMo.getWorkflow('qaoa', {'optimizer': optimizer, 'steps': 8})

# Execute and print the result
result = workflow.execute(problemModel)
energy = result['energy']
print(energy)

QAOA workflow returns the following information:

Quantum Imaginary Time Evolution - QITE ¶

The Quantum Imaginary Time Evolution (QITE) workflow requires the following configuration parameters:

Below is a simple example of using the QITE workflow.

#include "qcor_qsim.hpp"
int main(int argc, char **argv) {
  auto ham = 0.7071067811865475 * X(0) + 0.7071067811865475 * Z(0);
  // Number of QITE time steps and step size
  const int nbSteps = 25;
  const double stepSize = 0.1;
  auto problemModel = QuaSiMo::ModelFactory::createModel(ham);
  auto workflow =
      QuaSiMo::getWorkflow("qite", {{"steps", nbSteps}, {"step-size", stepSize}});
  auto result = workflow->execute(problemModel);
  const auto energy = result.get<double>("energy");
  // Array of energy values at each QITE step
  const auto energyAtStep = result.get<std::vector<double>>("exp-vals");
  std::cout << "Ground state energy: " << energy << "\n";
  return 0;
}

from qcor import * 

# Target H = 1/sqrt(2)(X + Z)
H = 0.7071067811865475 * X(0) + 0.7071067811865475 * Z(0)

# The number of Trotter steps 
nbSteps = 25

# The Trotter step size
stepSize = 0.1

# Create the problem model
problemModel = QuaSiMo.ModelFactory.createModel(H)

# Create the QITE workflow
workflow = QuaSiMo.getWorkflow('qite', {'steps': nbSteps, 'step-size': stepSize})

# Execute and print the result
result = workflow.execute(problemModel)

# Final energy:
energy = result['energy']
print('Ground state energy =', energy)

# Energy values along QITE steps
td_energy_vals = result['exp-vals']
print(td_energy_vals)

QITE workflow returns the following information:

Time-dependent Simulation ¶

The time-dependent simulation (td-evolution) workflow requires the following configuration parameters:

To execute a time-dependent simulation workflow, the problem model should provide a Hamiltonian function
taking time value as the input and returning the Hamiltonian operator.

This is illustrated in the below example.

#include "qcor_qsim.hpp"

int main(int argc, char **argv) {
  // Time-dependent Hamiltonian
  QuaSiMo::TdObservable H = [](double t) {
    // Parameters:
    const double Jz = 2 * M_PI * 2.86265 * 1e-3;
    const double epsilon = Jz; 
    const double omega = 4.8 * 2 * M_PI * 1e-3;
    return -Jz * Z(0) * Z(1) - Jz * Z(1) * Z(2) -
           epsilon * std::cos(omega * t) * (X(0) + X(1) + X(2));
  };

  // Observable = average magnetization
  auto observable = (1.0 / 3.0) * (Z(0) + Z(1) + Z(2));

  // Example: build model and TD workflow 
  auto problemModel = QuaSiMo::ModelFactory::createModel(observable, H);
  // Trotter step = 3fs, number of steps = 100 -> end time = 300fs
  auto workflow = QuaSiMo::getWorkflow(
      "td-evolution", {{"dt", 3.0}, {"steps", 100}});

  // Result should contain the observable expectation value along Trotter steps.
  auto result = workflow->execute(problemModel);
  // Get the observable values (average magnetization)
  const auto obsVals = result.get<std::vector<double>>("exp-vals");

  return 0;
}

from qcor import *
import numpy as np

# Time-dependent Hamiltonian: 
# Returns the Pauli operators at a time point.
def td_hamiltonian(t):
  Jz = 2 * np.pi * 2.86265 * 1e-3
  epsilon = Jz
  omega = 4.8 * 2 * np.pi * 1e-3
  return -Jz * Z(0) * Z(1)  - Jz * Z(1) * Z(2) + (-epsilon * np.cos(omega * t)) * (X(0) + X(1) + X(2)) 

# Observable = average magnetization
observable = (1.0 / 3.0) * (Z(0) + Z(1) + Z(2))

# Build model and TD workflow 
problemModel = QuaSiMo.ModelFactory.createModel(observable, td_hamiltonian)
# TD workflow with hyper-parameters: 
# Trotter step = 3fs, number of steps = 100 -> end time = 300fs
workflow = QuaSiMo.getWorkflow(
      "td-evolution", {"dt": 3.0, "steps": 100})

# Result contains the observable expectation value along Trotter steps.
result = workflow.execute(problemModel)
obsVals = result["exp-vals"]
print(obsVals)

Time-dependent simulation workflow returns the following information:

Cost Function (Observable) Evaluate ¶

A common utility that is required across QuaSiMo workflows is the ability to observe the expectation value of an arbitrary operator at a specific quantum state (e.g., after the ansatz circuit or after a certain number of Trotter steps).

This capability is provided by the CostFunctionEvaluator interface in the QuaSiMo library.

There are two CostFunctionEvaluator implementations currently available: the default one, which is based on the partial tomography method, and the quantum phase estimation (qpe)-based method, which can be optionally error-mitigated via the protocol described in this paper .

All workflows accept an optional evaluator key in addition to their set of configuration parameters to specify the CostFunctionEvaluator (via a name string or an instantiated instance). If none provided, the partial-tomography (default) one will be used.

The CostFunctionEvaluator configuration parameters, if any, can be provided in the workflow parameter pack, which will be forwarded to the CostFunctionEvaluator appropriately.

The CostFunctionEvaluator will use the runtime QPU instance to execute the necessary sub-circuits to evaluate the operator expectation value.

Partial Tomography ¶

The default CostFunctionEvaluator is based on the partial tomography method, whereby change-of-basis gates are added to measure the expectation values of each product term in the operator.

Users don’t need to provide the evaluator key (and set it to default) when creating the workflow if they wish to use this method because all workflows assume this option.

It’s worth noting that error mitigation (if any) for this cost evaluator is performed at the QPU level via conventional QCOR error mitigation decorators.

Quantum Phase Estimation ¶

The qpe cost evaluator is based on the time-series quantum phase estimation method to estimate the expectation value. In particular, it makes use of the Prony signal processing method to extract the expectation value estimate from the time-series QPE signal.

The qpe cost evaluator has the following configuration parameters:

The qpe cost evaluator can be explicitly instantiated by using the getObjEvaluator helper function and providing the
operator to be evaluated, the name key (qpe), and any additional parameters.

For example,

• In C++:

auto observable = Z(0) + Z(1) + Z(2);
auto vqpeEvaluator = QuaSiMo::getObjEvaluator(observable, "qpe", {{"verified", true}});

• In Python:

observable = Z(0) + Z(1) + Z(2)
vqpeEvaluator = QuaSiMo.getObjEvaluator(observable, 'qpe', {'verified': True})

Implement a new workflow ¶

Please refer to this section for a tutorial on how to create a new QuaSiMo workflow plugin.