Estimate expectation value of operators

One of the most basic procedures in quantum algorithms is measuring operators with a given quantum state. There are several ways to measure operators. We start with estimation of expectation value of operators since it is a basic operation found in many NISQ algorithms.

Prerequisite

QURI Parts modules used in this tutorial: quri-parts-core and quri-parts-qulacs. You can install them as follows:

!pip install "quri-parts[qulacs]"

Pauli label

A Pauli label (also known as Pauli string) is a product of Pauli matrices for single qubits. A Pauli label \(X_0 Y_2 Z_4\), a product of \(X\) on qubit 0, \(Y\) on qubit 2 and \(Z\) on qubit 4, can be defined as follows:

from quri_parts.core.operator import pauli_label
label = pauli_label("X0 Y2 Z4")
# You can put spaces between a Pauli name and a qubit index
label = pauli_label("X 0 Y 2 Z 4")
print(label)

X0 Y2 Z4

PAULI_IDENTITY represents a Pauli label with zero Pauli matrices:

from quri_parts.core.operator import PAULI_IDENTITY
print(PAULI_IDENTITY)

I

A Pauli label is immutable and hashable. It is a (frozen) set of pairs of a qubit index and a SinglePauli enum, and you can iterate over it:

for pair in label:
    print(pair)

(0, <SinglePauli.X: 1>)
(4, <SinglePauli.Z: 3>)
(2, <SinglePauli.Y: 2>)

for index, matrix in label:
    print(f"qubit index: {index}, Pauli matrix: {matrix}")

qubit index: 0, Pauli matrix: 1
qubit index: 4, Pauli matrix: 3
qubit index: 2, Pauli matrix: 2

SinglePauli is an IntEnum. You can either use SinglePauli.X, SinglePauli.Y and SinglePauli.Z, or just integers 1, 2, 3:

from quri_parts.core.operator import SinglePauli
print(SinglePauli.X == 1)
print(SinglePauli.Y == 2)
print(SinglePauli.Z == 3)

True
True
True

Operator

An operator is defined as a set of pairs of a Pauli label and its complex coefficient. For example, an operator \((0.5 + 0.5i) X_0 Y_1 + 0.2i Z_0 Z_2 + 0.3 + 0.4i\) can be defined as follows:

from quri_parts.core.operator import Operator
op = Operator({
    pauli_label("X0 Y1"): 0.5 + 0.5j,
    pauli_label("Z0 Z2"): 0.2j,
    PAULI_IDENTITY: 0.3 + 0.4j,
})
print(op)

(0.5+0.5j)*X0 Y1 + 0.2j*Z0 Z2 + (0.3+0.4j)*I

You can also construct an operator by adding single Pauli terms:

op = Operator()
op.add_term(pauli_label("X0 Y1"), 0.5 + 0.5j)
op.add_term(pauli_label("Z0 Z2"), 0.2j)
op.constant = 0.3 + 0.4j
print(op)
print(f"Number of terms: {op.n_terms}")

# When adding an existing term, its coefficient is summed
op.add_term(pauli_label("X0 Y1"), 0.5)
print(op)
print(f"Number of terms: {op.n_terms}")

# If the sum of coefficient is zero, the term gets removed
op.add_term(pauli_label("X0 Y1"), -1.0 - 0.5j)
print(op)
print(f"Number of terms: {op.n_terms}")

(0.5+0.5j)*X0 Y1 + 0.2j*Z0 Z2 + (0.3+0.4j)*I
Number of terms: 3
(1+0.5j)*X0 Y1 + 0.2j*Z0 Z2 + (0.3+0.4j)*I
Number of terms: 3
0.2j*Z0 Z2 + (0.3+0.4j)*I
Number of terms: 2

To get a hermitian conjugated operator:

conj = op.hermitian_conjugated()
print(conj)

-0.2j*Z0 Z2 + (0.3-0.4j)*I

An Operator is internally a dict, so you can do the followings:

p = pauli_label("Z0 Z2")
coef = op[p]
print(f"Coefficient of {p} = {coef}")

op[p] = 0.4
coef = op[p]
print(f"Coefficient of {p} = {coef}")

for label, coef in op.items():
    print(f"Coefficient of {label} = {coef}")

Coefficient of Z0 Z2 = 0.2j
Coefficient of Z0 Z2 = 0.4
Coefficient of Z0 Z2 = 0.4
Coefficient of I = (0.3+0.4j)

Quantum state

Several types of quantum states are available. Here we introduce the most basic states: computational basis states. Other types of states will be described later.

A computational basis state is a quantum state where each qubit is in 0 or 1 eigenstate. To construct a computational basis state for 5 qubits:

from quri_parts.core.state import ComputationalBasisState
state1 = ComputationalBasisState(5, bits=0b10100)
print(state1)

ComputationalBasisState(qubit_count=5, bits=0b10100, phase=0π/2)

Here bits=0b10100 means that qubit 0 is in \(|0\rangle\), qubit 1 is in \(|0\rangle\), qubit 2 is in \(|1\rangle\), qubit 3 is in \(|0\rangle\) and qubit 4 is in \(|1\rangle\). We use 0-based indices for qubits, and bits for qubits are ordered from the least significant bit to the most significant bit.

You can also create a superposition of two computational basis states. Note that the resulting state is not a computational basis state anymore. comp_basis_superposition() takes four arguments. First two are computational basis states to be superposed. The third argument \(\theta\) determines the weight for the superposition and the fourth argument \(\phi\) determines the phase factor for the superposition: coefficients for two states are given as \(\cos\theta\) and \(e^{i\phi}\sin\theta\).

import math
from quri_parts.core.state import comp_basis_superposition
state2 = ComputationalBasisState(5, bits=0b01011)
sp_state = comp_basis_superposition(state1, state2, math.pi/2, math.pi/4)
print(sp_state)

GeneralCircuitQuantumState(n_qubits=5, circuit=<quri_parts.circuit.circuit.ImmutableQuantumCircuit object at 0x7f78a3de2ee0>)

Estimator

To estimate the expectation value of an operator with a given state, you can use a QuantumEstimator. The QuantumEstimator itself (defined in quri_parts.core.estimator package) is an abstract interface and you need a concrete instance for actual estimation. Various implementations of QuantumEstimator interface can use various methods for estimation: direct calculation with state vector simulation, sampling simulation, sampling on real devices, or more advanced ones.

One of the most basic estimators is one that calculates the expectation value directly by state vector simulation with Qulacs:

from quri_parts.qulacs.estimator import create_qulacs_vector_estimator
# First create the estimator
estimator = create_qulacs_vector_estimator()
# Pass the operator and state to the estimator
estimate = estimator(op, sp_state)
# The return value contains the estimated value and the estimation error
print(f"Estimated expectation value: {estimate.value}")
# (The estimation error is zero for state vector calculation)
print(f"Estimation error: {estimate.error}")

Estimated expectation value: (-0.10000000000000003+0.4j)
Estimation error: 0.0

The return value of a QuantumEstimator contains the estimated expectation value (.value) and the estimation error (.error).

There is also a ConcurrentQuantumEstimator interface, which estimates multiple operators and multiple states at once. The interface accept either one of the followings:

• One operator, multiple states

• Multiple operators, one state

• The same number of operators and states

For the case of Qulacs, you can create a ConcurrentQuantumEstimator specifying a concurrent.futures.Executor (default is None, meaning no parallelization) and concurrency (default is 1). Note that since Qulacs itself has multithreading support, using ThreadPoolExecutor or ProcessPoolExecutor may not have performance improvement.

from concurrent.futures import ThreadPoolExecutor
from quri_parts.qulacs.estimator import create_qulacs_vector_concurrent_estimator
# Create an executor (optional)
executor = ThreadPoolExecutor(max_workers=4)
# Create the concurrent estimator
estimator = create_qulacs_vector_concurrent_estimator(executor, concurrency=4)

estimates = estimator(op, [state1, state2])
for i, est in enumerate(estimates):
    print(f"State {i}: value={est.value}, error={est.error}")

# You can pass a PauliLabel to an estimator too
estimates = estimator([label, op], [sp_state])
for i, est in enumerate(estimates):
    print(f"Operator {i}: value={est.value}, error={est.error}")

estimates = estimator([label, op], [state1, state2])
for i, est in enumerate(estimates):
    print(f"Operator {i}, state {i}: value={est.value}, error={est.error}")

State 0: value=(-1+0j), error=0.0
State 1: value=(1+0j), error=0.0
Operator 0: value=(1+0j), error=0.0
Operator 1: value=(-0.10000000000000003+0.4j), error=0.0
Operator 0, state 0: value=(1+0j), error=0.0
Operator 1, state 1: value=(-0.10000000000000003+0.4j), error=0.0