Note

This page was generated from tutorials/operators/01_operator_flow.ipynb.

Operator Flow

Introduction

Qiskit provides classes representing states and operators and sums, tensor products, and compositions thereof. These algebraic constructs allow us to build expressions representing operators.

We introduce expressions by building them from Pauli operators. In subsequent sections we explore in more detail operators and states, how they are represented, and what we can do with them. In the last section we construct a state, evolve it with a Hamiltonian, and compute expectation values of an observable.

Pauli operators, sums, compositions, and tensor products

The most important base operators are the Pauli operators. The Pauli operators are represented like this.

from qiskit.opflow import I, X, Y, Z
print(I, X, Y, Z)

I X Y Z

These operators may also carry a coefficient.

print(1.5 * I)
print(2.5 * X)

1.5 * I
2.5 * X

These coefficients allow the operators to be used as terms in a sum.

print(X + 2.0 * Y)

1.0 * X
+ 2.0 * Y

Tensor products are denoted with a caret, like this.

print(X^Y^Z)

XYZ

Composition is denoted by the @ symbol.

print(X @ Y @ Z)

iI

In the preceding two examples, the tensor product and composition of Pauli operators were immediately reduced to the equivalent (possibly multi-qubit) Pauli operator. If we tensor or compose more complicated objects, the result is objects representing the unevaluated operations. That is, algebraic expressions.

For example, composing two sums gives

print((X + Y) @ (Y + Z))

1j * Z
+ -1j * Y
+ 1.0 * I
+ 1j * X

And tensoring two sums gives

print((X + Y) ^ (Y + Z))

1.0 * XY
+ 1.0 * XZ
+ 1.0 * YY
+ 1.0 * YZ

Let’s take a closer look at the types introduced above. First the Pauli operators.

(I, X)

(PauliOp(Pauli('I'), coeff=1.0), PauliOp(Pauli('X'), coeff=1.0))

Each Pauli operator is an instance of PauliOp, which wraps an instance of qiskit.quantum_info.Pauli, and adds a coefficient coeff. In general, a PauliOp represents a weighted tensor product of Pauli operators.

2.0 * X^Y^Z

PauliOp(Pauli('XYZ'), coeff=2.0)

For the encoding of the Pauli operators as pairs of Boolean values, see the documentation for qiskit.quantum_info.Pauli.

All of the objects representing operators, whether as “primitive”s such as PauliOp, or algebraic expressions carry a coefficient

print(1.1 * ((1.2 * X)^(Y + (1.3 * Z))))

1.2 * (
  1.1 * XY
  + 1.4300000000000002 * XZ
)

In the following we take a broader and deeper look at Qiskit’s operators, states, and the building blocks of quantum algorithms.

Part I: State Functions and Measurements

Quantum states are represented by subclasses of the class StateFn. There are four representations of quantum states: DictStateFn is a sparse representation in the computational basis, backed by a dict. VectorStateFn is a dense representation in the computational basis backed by a numpy array. CircuitStateFn is backed by a circuit and represents the state obtained by executing the circuit on the all-zero computational-basis state. OperatorStateFn represents mixed states via a density matrix. (As we will see later, OperatorStateFn is also used to represent observables.)

Several StateFn instances are provided for convenience. For example Zero, One, Plus, Minus.

from qiskit.opflow import (StateFn, Zero, One, Plus, Minus, H,
                           DictStateFn, VectorStateFn, CircuitStateFn, OperatorStateFn)

Zero and One represent the quantum states \(|0\rangle\) and \(|1\rangle\). They are represented via DictStateFn.

print(Zero, One)

DictStateFn({'0': 1}) DictStateFn({'1': 1})

Plus and Minus, representing states \((|0\rangle + |1\rangle)/\sqrt{2}\) and \((|0\rangle - |1\rangle)/\sqrt{2}\) are represented via circuits. H is a synonym for Plus.

print(Plus, Minus)

CircuitStateFn(
     ┌───┐
q_0: ┤ H ├
     └───┘
) CircuitStateFn(
     ┌───┐┌───┐
q_0: ┤ X ├┤ H ├
     └───┘└───┘
)

Indexing into quantum states is done with the eval method. These examples return the coefficients of the 0 and 1 basis states. (Below, we will see that the eval method is used for other computations, as well.)

print(Zero.eval('0'))
print(Zero.eval('1'))
print(One.eval('1'))
print(Plus.eval('0'))
print(Minus.eval('1'))

1.0
0.0
1.0
(0.7071067811865475+0j)
(-0.7071067811865475+8.7e-17j)

The dual vector of a quantum state, that is the bra corresponding to a ket is obtained via the adjoint method. The StateFn carries a flag is_measurement, which is False if the object is a ket and True if it is a bra.

Here, we construct \(\langle 1 |\).

One.adjoint()

DictStateFn({'1': 1}, coeff=1.0, is_measurement=True)

For convenience, one may obtain the dual vector with a tilde, like this

~One

DictStateFn({'1': 1}, coeff=1.0, is_measurement=True)

Algebraic operations and predicates

Many algebraic operations and predicates between StateFns are supported, including:

• + - addition

• - - subtraction, negation (scalar multiplication by -1)

• * - scalar multiplication

• / - scalar division

• @ - composition

• ^ - tensor product or tensor power (tensor with self n times)

• ** - composition power (compose with self n times)

• == - equality

• ~ - adjoint, alternating between a State Function and Measurement

Be very aware that these operators obey the Python rules for operator precedence, which might not be what you expect mathematically. For example, I^X + X^I will actually be parsed as I ^ (X + X) ^ I == 2 * (I^X^I) because Python evaluates + before ^. In these cases, you can use the methods (.tensor(), etc) or parentheses.

StateFns carry a coefficient. This allows us multiply states by a scalar, and so to construct sums.

Here, we construct \((2 + 3i)|0\rangle\).

(2.0 + 3.0j) * Zero

DictStateFn({'0': 1}, coeff=(2+3j), is_measurement=False)

Here, we see that adding two DictStateFns returns an object of the same type. We construct \(|0\rangle + |1\rangle\).

print(Zero + One)

DictStateFn({'0': 1.0, '1': 1.0})

Note that you must normalize states by hand. For example, to construct \((|0\rangle + |1\rangle)/\sqrt{2}\), we write

import math

v_zero_one = (Zero + One) / math.sqrt(2)
print(v_zero_one)

DictStateFn({'0': 1.0, '1': 1.0}) * 0.7071067811865475

In other cases, the result is a symbolic representation of a sum. For example, here is a representation of \(|+\rangle + |-\rangle\).

print(Plus + Minus)

SummedOp([
  CircuitStateFn(
       ┌───┐
  q_0: ┤ H ├
       └───┘
  ),
  CircuitStateFn(
       ┌───┐┌───┐
  q_0: ┤ X ├┤ H ├
       └───┘└───┘
  )
])

The composition operator is used to perform an inner product, which by default is held in an unevaluated form. Here is a representation of \(\langle 1 | 1 \rangle\).

print(~One @ One)

ComposedOp([
  DictMeasurement({'1': 1}),
  DictStateFn({'1': 1})
])

Note that the is_measurement flag causes the (bra) state ~One to be printed DictMeasurement.

Symbolic expressions may be evaluated with the eval method.

(~One @ One).eval()

1.0

(~v_zero_one @ v_zero_one).eval()

0.9999999999999998

Here is \(\langle - | 1 \rangle = \langle (\langle 0| - \langle 1|)/\sqrt{2} | 1\rangle\).

(~Minus @ One).eval()

(-0.7071067811865475-8.7e-17j)

The composition operator @ is equivalent to calling the compose method.

print((~One).compose(One))

ComposedOp([
  DictMeasurement({'1': 1}),
  DictStateFn({'1': 1})
])

Inner products may also be computed using the eval method directly, without constructing a ComposedOp.

(~One).eval(One)

1.0

Symbolic tensor products are constructed as follows. Here is \(|0\rangle \otimes |+\rangle\).

print(Zero^Plus)

TensoredOp([
  DictStateFn({'0': 1}),
  CircuitStateFn(
       ┌───┐
  q_0: ┤ H ├
       └───┘
  )
])

This may be represented as a simple (not compound) CircuitStateFn.

print((Zero^Plus).to_circuit_op())

CircuitStateFn(
     ┌───┐
q_0: ┤ H ├
     └───┘
q_1: ─────

)

Tensor powers are constructed using the caret ^ as follows. Here are \(600 (|11111\rangle + |00000\rangle)\), and \(|10\rangle^{\otimes 3}\).

print(600 * ((One^5) + (Zero^5)))
print((One^Zero)^3)

DictStateFn({'11111': 1.0, '00000': 1.0}) * 600.0
DictStateFn({'101010': 1})

The method to_matrix_op converts to VectorStateFn.

print(((Plus^Minus)^2).to_matrix_op())
print(((Plus^One)^2).to_circuit_op())
print(((Plus^One)^2).to_matrix_op().sample())

VectorStateFn(Statevector([ 0.25-6.1e-17j, -0.25+6.1e-17j,  0.25-6.1e-17j,
             -0.25+6.1e-17j, -0.25+6.1e-17j,  0.25-6.1e-17j,
             -0.25+6.1e-17j,  0.25-6.1e-17j,  0.25-6.1e-17j,
             -0.25+6.1e-17j,  0.25-6.1e-17j, -0.25+6.1e-17j,
             -0.25+6.1e-17j,  0.25-6.1e-17j, -0.25+6.1e-17j,
              0.25-6.1e-17j],
            dims=(2, 2, 2, 2)))
CircuitStateFn(
     ┌───┐
q_0: ┤ X ├
     ├───┤
q_1: ┤ H ├
     ├───┤
q_2: ┤ X ├
     ├───┤
q_3: ┤ H ├
     └───┘
)
{'0101': 0.265625, '1101': 0.2568359375, '1111': 0.2490234375, '0111': 0.228515625}

Constructing a StateFn is easy. The StateFn class also serves as a factory, and can take any applicable primitive in its constructor and return the correct StateFn subclass. Currently the following primitives can be passed into the constructor, listed alongside the StateFn subclass they produce:

• str (equal to some basis bitstring) -> DictStateFn

• dict -> DictStateFn

• Qiskit Result object -> DictStateFn

• list -> VectorStateFn

• np.ndarray -> VectorStateFn

• Statevector -> VectorStateFn

• QuantumCircuit -> CircuitStateFn

• Instruction -> CircuitStateFn

• OperatorBase -> OperatorStateFn

print(StateFn({'0':1}))
print(StateFn({'0':1}) == Zero)

print(StateFn([0,1,1,0]))

from qiskit.circuit.library import RealAmplitudes
print(StateFn(RealAmplitudes(2)))

DictStateFn({'0': 1})
True
VectorStateFn(Statevector([0.+0.j, 1.+0.j, 1.+0.j, 0.+0.j],
            dims=(2, 2)))
CircuitStateFn(
     ┌──────────────────────────────────────────────────────────┐
q_0: ┤0                                                         ├
     │  RealAmplitudes(θ[0],θ[1],θ[2],θ[3],θ[4],θ[5],θ[6],θ[7]) │
q_1: ┤1                                                         ├
     └──────────────────────────────────────────────────────────┘
)

/home/computertreker/git/qiskit/qiskit/.tox/docs/lib/python3.7/site-packages/sympy/core/expr.py:3951: SymPyDeprecationWarning:

expr_free_symbols method has been deprecated since SymPy 1.9. See
https://github.com/sympy/sympy/issues/21494 for more info.

  deprecated_since_version="1.9").warn()

Part II: PrimitiveOps

The basic Operators are subclasses of PrimitiveOp. Just like StateFn, PrimitiveOp is also a factory for creating the correct type of PrimitiveOp for a given primitive. Currently, the following primitives can be passed into the constructor, listed alongside the PrimitiveOp subclass they produce:

• Terra’s Pauli -> PauliOp

• Instruction -> CircuitOp

• QuantumCircuit -> CircuitOp

• 2d List -> MatrixOp

• np.ndarray -> MatrixOp

• spmatrix -> MatrixOp

• Terra’s quantum_info.Operator -> MatrixOp

from qiskit.opflow import X, Y, Z, I, CX, T, H, S, PrimitiveOp

Matrix elements

The eval method returns a column from an operator. For example, the Pauli \(X\) operator is represented by a PauliOp. Asking for a column returns an instance of the sparse representation, a DictStateFn.

X

PauliOp(Pauli('X'), coeff=1.0)

print(X.eval('0'))

DictStateFn({'1': (1+0j)})

It follows that indexing into an operator, that is obtaining a matrix element, is performed with two calls to the eval method.

We have \(X = \left(\begin{matrix} 0 & 1 \\ 1 & 0 \end{matrix} \right)\). And the matrix element \(\left\{X \right\}_{0,1}\) is

X.eval('0').eval('1')

(1+0j)

Here is an example using the two qubit operator CX, the controlled X, which is represented by a circuit.

print(CX)
print(CX.to_matrix().real) # The imaginary part vanishes.

q_0: ──■──
     ┌─┴─┐
q_1: ┤ X ├
     └───┘
[[1. 0. 0. 0.]
 [0. 0. 0. 1.]
 [0. 0. 1. 0.]
 [0. 1. 0. 0.]]

CX.eval('01')  # 01 is the one in decimal. We get the first column.

VectorStateFn(Statevector([0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j],
            dims=(2, 2)), coeff=1.0, is_measurement=False)

CX.eval('01').eval('11')  # This returns element with (zero-based) index (1, 3)

(1+0j)

Applying an operator to a state vector

Applying an operator to a state vector may be done with the compose method (equivalently, @ operator). Here is a representation of \(X | 1 \rangle = |0\rangle\).

print(X @ One)

ComposedOp([
  X,
  DictStateFn({'1': 1})
])

A simpler representation, the DictStateFn representation of \(|0\rangle\), is obtained with eval.

(X @ One).eval()

DictStateFn({'0': (1+0j)}, coeff=(1+0j), is_measurement=False)

The intermediate ComposedOp step may be avoided by using eval directly.

X.eval(One)

DictStateFn({'0': (1+0j)}, coeff=(1+0j), is_measurement=False)

Composition and tensor products of operators are effected with @ and ^. Here are some examples.

print(((~One^2) @ (CX.eval('01'))).eval())

print(((H^5) @ ((CX^2)^I) @ (I^(CX^2)))**2)
print((((H^5) @ ((CX^2)^I) @ (I^(CX^2)))**2) @ (Minus^5))
print(((H^I^I)@(X^I^I)@Zero))

(1+0j)
                    ┌─────────────┐┌───┐               ┌─────────────┐┌───┐
q_0: ───────■───────┤ circuit-205 ├┤ H ├───────■───────┤ circuit-205 ├┤ H ├
          ┌─┴─┐     └─────────────┘├───┤     ┌─┴─┐     └─────────────┘├───┤
q_1: ─────┤ X ├────────────■───────┤ H ├─────┤ X ├────────────■───────┤ H ├
          └───┘          ┌─┴─┐     ├───┤     └───┘          ┌─┴─┐     ├───┤
q_2: ───────■────────────┤ X ├─────┤ H ├───────■────────────┤ X ├─────┤ H ├
          ┌─┴─┐          └───┘     ├───┤     ┌─┴─┐          └───┘     ├───┤
q_3: ─────┤ X ├────────────■───────┤ H ├─────┤ X ├────────────■───────┤ H ├
     ┌────┴───┴────┐     ┌─┴─┐     ├───┤┌────┴───┴────┐     ┌─┴─┐     ├───┤
q_4: ┤ circuit-218 ├─────┤ X ├─────┤ H ├┤ circuit-218 ├─────┤ X ├─────┤ H ├
     └─────────────┘     └───┘     └───┘└─────────────┘     └───┘     └───┘
CircuitStateFn(
     ┌───┐┌───┐     ┌───┐          ┌───┐
q_0: ┤ X ├┤ H ├──■──┤ H ├───────■──┤ H ├─────
     ├───┤├───┤┌─┴─┐└───┘┌───┐┌─┴─┐└───┘┌───┐
q_1: ┤ X ├┤ H ├┤ X ├──■──┤ H ├┤ X ├──■──┤ H ├
     ├───┤├───┤└───┘┌─┴─┐├───┤└───┘┌─┴─┐├───┤
q_2: ┤ X ├┤ H ├──■──┤ X ├┤ H ├──■──┤ X ├┤ H ├
     ├───┤├───┤┌─┴─┐└───┘├───┤┌─┴─┐└───┘├───┤
q_3: ┤ X ├┤ H ├┤ X ├──■──┤ H ├┤ X ├──■──┤ H ├
     ├───┤├───┤└───┘┌─┴─┐├───┤└───┘┌─┴─┐├───┤
q_4: ┤ X ├┤ H ├─────┤ X ├┤ H ├─────┤ X ├┤ H ├
     └───┘└───┘     └───┘└───┘     └───┘└───┘
)
CircuitStateFn(
     ┌─────────────┐
q_0: ┤0            ├─────
     │             │
q_1: ┤1 Pauli(XII) ├─────
     │             │┌───┐
q_2: ┤2            ├┤ H ├
     └─────────────┘└───┘
)

print(~One @ Minus)

ComposedOp([
  DictMeasurement({'1': 1}),
  CircuitStateFn(
       ┌───┐┌───┐
  q_0: ┤ X ├┤ H ├
       └───┘└───┘
  )
])

Part III ListOp and subclasses

ListOp

ListOp is a container for effectively vectorizing operations over a list of operators and states.

from qiskit.opflow import ListOp

print((~ListOp([One, Zero]) @ ListOp([One, Zero])))

ComposedOp([
  ListOp([
    DictMeasurement({'1': 1}),
    DictMeasurement({'0': 1})
  ]),
  ListOp([
    DictStateFn({'1': 1}),
    DictStateFn({'0': 1})
  ])
])

For example, the composition above is distributed over the lists (ListOp) using the simplification method reduce.

print((~ListOp([One, Zero]) @ ListOp([One, Zero])).reduce())

ListOp([
  ListOp([
    ComposedOp([
      DictMeasurement({'1': 1}),
      DictStateFn({'1': 1})
    ]),
    ComposedOp([
      DictMeasurement({'1': 1}),
      DictStateFn({'0': 1})
    ])
  ]),
  ListOp([
    ComposedOp([
      DictMeasurement({'0': 1}),
      DictStateFn({'1': 1})
    ]),
    ComposedOp([
      DictMeasurement({'0': 1}),
      DictStateFn({'0': 1})
    ])
  ])
])

ListOps: SummedOp, ComposedOp, TensoredOp

ListOp, introduced above, is useful for vectorizing operations. But, it also serves as the superclass for list-like composite classes. If you’ve already played around with the above, you’ll notice that you can easily perform operations between OperatorBases which we may not know how to perform efficiently in general (or simply haven’t implemented an efficient procedure for yet), such as addition between CircuitOps. In those cases, you may receive a ListOp result (or subclass thereof) from your operation representing the lazy execution of the operation. For example, if you attempt to add together a DictStateFn and a CircuitStateFn, you’ll receive a SummedOp representing the sum of the two. This composite State function still has a working eval (but may need to perform a non-scalable computation under the hood, such as converting both to vectors).

These composite OperatorBases are how we construct increasingly complex and rich computation out of PrimitiveOp and StateFn building blocks.

Every ListOp has four properties: * oplist - The list of OperatorBases which may represent terms, factors, etc. * combo_fn - The function taking a list of complex numbers to an output value which defines how to combine the outputs of the oplist items. For broadcasting simplicity, this function is defined over NumPy arrays. * coeff - A coefficient multiplying the primitive. Note that coeff can be int, float, complex or a free Parameter object (from qiskit.circuit in Terra) to be bound later using my_op.bind_parameters. * abelian - Indicates whether the Operators in oplist are known to mutually commute (usually set after being converted by the AbelianGrouper converter).

Note that ListOp supports typical iteration overloads, so you can use indexing like my_op[4] to access the OperatorBases in oplist.

OperatorStateFn

We mentioned above that OperatorStateFn represents a density operator. But, if the is_measurement flag is True, then OperatorStateFn represents an observable. The expectation value of this observable can then be constructed via ComposedOp. Or, directly, using eval. Recall that the is_measurement flag (property) is set via the adjoint method.

Here we construct the observable corresponding to the Pauli \(Z\) operator. Note that when printing, it is called OperatorMeasurement.

print(StateFn(Z).adjoint())
StateFn(Z).adjoint()

OperatorMeasurement(Z)

OperatorStateFn(PauliOp(Pauli('Z'), coeff=1.0), coeff=1.0, is_measurement=True)

Here, we compute \(\langle 0 | Z | 0 \rangle\), \(\langle 1 | Z | 1 \rangle\), and \(\langle + | Z | + \rangle\), where \(|+\rangle = (|0\rangle + |1\rangle)/\sqrt{2}\).

print(StateFn(Z).adjoint().eval(Zero))
print(StateFn(Z).adjoint().eval(One))
print(StateFn(Z).adjoint().eval(Plus))

(1+0j)
(-1+0j)
0j

Part IV: Converters

Converters are classes that manipulate operators and states and perform building blocks of algorithms. Examples include changing the basis of operators and Trotterization. Converters traverse an expression and perform a particular manipulation or replacement, defined by the converter’s convert() method, of the Operators within. Typically, if a converter encounters an OperatorBase in the recursion which is irrelevant to its conversion purpose, that OperatorBase is left unchanged.

import numpy as np
from qiskit.opflow import I, X, Y, Z, H, CX, Zero, ListOp, PauliExpectation, PauliTrotterEvolution, CircuitSampler, MatrixEvolution, Suzuki
from qiskit.circuit import Parameter
from qiskit import Aer

Evolutions, exp_i(), and the EvolvedOp

Every PrimitiveOp and ListOp has an .exp_i() function such that H.exp_i() corresponds to \(e^{-iH}\). In practice, only a few of these Operators have an efficiently computable exponentiation (such as MatrixOp and the PauliOps with only one non-identity single-qubit Pauli), so we need to return a placeholder, or symbolic representation, (similar to how SummedOp is a placeholder when we can’t perform addition). This placeholder is called EvolvedOp, and it holds the OperatorBase to be exponentiated in its .primitive property.

Qiskit operators fully support parameterization, so we can use a Parameter for our evolution time here. Notice that there’s no “evolution time” argument in any function. The Operator flow exponentiates whatever operator we tell it to, and if we choose to multiply the operator by an evolution time, \(e^{-iHt}\), that will be reflected in our exponentiation parameters.

Weighted sum of Pauli operators

A Hamiltonian expressed as a linear combination of multi-qubit Pauli operators may be constructed like this.

two_qubit_H2 =  (-1.0523732 * I^I) + \
                (0.39793742 * I^Z) + \
                (-0.3979374 * Z^I) + \
                (-0.0112801 * Z^Z) + \
                (0.18093119 * X^X)

Note that two_qubit_H2 is represented as a SummedOp whose terms are PauliOps.

print(two_qubit_H2)

-1.0523732 * II
+ 0.39793742 * IZ
- 0.3979374 * ZI
- 0.0112801 * ZZ
+ 0.18093119 * XX

Next, we multiply the Hamiltonian by a Parameter. This Parameter is stored in the coeff property of the SummedOp. Calling exp_i() on the result wraps it in EvolvedOp, representing exponentiation.

evo_time = Parameter('θ')
evolution_op = (evo_time*two_qubit_H2).exp_i()
print(evolution_op) # Note, EvolvedOps print as exponentiations
print(repr(evolution_op))

e^(-i*1.0*θ * (
  -1.0523732 * II
  + 0.39793742 * IZ
  - 0.3979374 * ZI
  - 0.0112801 * ZZ
  + 0.18093119 * XX
))
EvolvedOp(PauliSumOp(SparsePauliOp([[False, False, False, False],
               [False, False,  True, False],
               [False, False, False,  True],
               [False, False,  True,  True],
               [ True,  True, False, False]],
              coeffs=[-1.0523732 +0.j,  0.39793742+0.j, -0.3979374 +0.j, -0.0112801 +0.j,
  0.18093119+0.j]), coeff=1.0*θ), coeff=1.0)

We construct h2_measurement, which represents two_qubit_H2 as an observable.

h2_measurement = StateFn(two_qubit_H2).adjoint()
print(h2_measurement)

OperatorMeasurement(-1.0523732 * II
+ 0.39793742 * IZ
- 0.3979374 * ZI
- 0.0112801 * ZZ
+ 0.18093119 * XX)

We construct a Bell state \(|\Phi_+\rangle\) via \(\text{CX} (H\otimes I) |00\rangle\).

bell = CX @ (I ^ H) @ Zero
print(bell)

CircuitStateFn(
     ┌───┐
q_0: ┤ H ├──■──
     └───┘┌─┴─┐
q_1: ─────┤ X ├
          └───┘
)

Here is the expression \(H e^{-iHt} |\Phi_+\rangle\).

evo_and_meas = h2_measurement @ evolution_op @ bell
print(evo_and_meas)

ComposedOp([
  OperatorMeasurement(-1.0523732 * II
  + 0.39793742 * IZ
  - 0.3979374 * ZI
  - 0.0112801 * ZZ
  + 0.18093119 * XX),
  e^(-i*1.0*θ * (
    -1.0523732 * II
    + 0.39793742 * IZ
    - 0.3979374 * ZI
    - 0.0112801 * ZZ
    + 0.18093119 * XX
  )),
  CircuitStateFn(
       ┌───┐
  q_0: ┤ H ├──■──
       └───┘┌─┴─┐
  q_1: ─────┤ X ├
            └───┘
  )
])

Typically, we want to approximate \(e^{-iHt}\) using two-qubit gates. We achieve this with the convert method of PauliTrotterEvolution, which traverses expressions applying trotterization to all EvolvedOps encountered. Although we use PauliTrotterEvolution here, there are other possibilities, such as MatrixEvolution, which performs the exponentiation exactly.

trotterized_op = PauliTrotterEvolution(trotter_mode=Suzuki(order=2, reps=1)).convert(evo_and_meas)
# We can also set trotter_mode='suzuki' or leave it empty to default to first order Trotterization.
print(trotterized_op)

ComposedOp([
  OperatorMeasurement(-1.0523732 * II
  + 0.39793742 * IZ
  - 0.3979374 * ZI
  - 0.0112801 * ZZ
  + 0.18093119 * XX),
  CircuitStateFn(
  global phase: 1.0523732*θ
       ┌───┐     ┌───┐┌───┐┌──────────────────┐┌───┐┌───┐┌───┐»
  q_0: ┤ H ├──■──┤ H ├┤ X ├┤ Rz(0.18093119*θ) ├┤ X ├┤ H ├┤ X ├»
       └───┘┌─┴─┐├───┤└─┬─┘└──────────────────┘└─┬─┘├───┤└─┬─┘»
  q_1: ─────┤ X ├┤ H ├──■────────────────────────■──┤ H ├──■──»
            └───┘└───┘                              └───┘     »
  «     ┌──────────────────┐┌───┐┌──────────────────┐┌──────────────────┐┌───┐»
  «q_0: ┤ Rz(-0.0112801*θ) ├┤ X ├┤ Rz(0.39793742*θ) ├┤ Rz(0.39793742*θ) ├┤ X ├»
  «     └──────────────────┘└─┬─┘├──────────────────┤├──────────────────┤└─┬─┘»
  «q_1: ──────────────────────■──┤ Rz(-0.3979374*θ) ├┤ Rz(-0.3979374*θ) ├──■──»
  «                              └──────────────────┘└──────────────────┘     »
  «     ┌──────────────────┐┌───┐┌───┐┌───┐┌──────────────────┐┌───┐┌───┐
  «q_0: ┤ Rz(-0.0112801*θ) ├┤ X ├┤ H ├┤ X ├┤ Rz(0.18093119*θ) ├┤ X ├┤ H ├
  «     └──────────────────┘└─┬─┘├───┤└─┬─┘└──────────────────┘└─┬─┘├───┤
  «q_1: ──────────────────────■──┤ H ├──■────────────────────────■──┤ H ├
  «                              └───┘                              └───┘
  )
])

trotterized_op contains a Parameter. The bind_parameters method traverses the expression binding values to parameter names as specified via a dict. In this case, there is only one parameter.

bound = trotterized_op.bind_parameters({evo_time: .5})

bound is a ComposedOp. The second factor is the circuit. Let’s draw it to verify that the binding has taken place.

bound[1].to_circuit().draw()

global phase: 0.52619
     ┌───┐     ┌───┐┌───┐┌────────────────────────┐┌───┐┌───┐┌───┐»
q_0: ┤ H ├──■──┤ H ├┤ X ├┤ Rz(0.0904655950000000) ├┤ X ├┤ H ├┤ X ├»
     └───┘┌─┴─┐├───┤└─┬─┘└────────────────────────┘└─┬─┘├───┤└─┬─┘»
q_1: ─────┤ X ├┤ H ├──■──────────────────────────────■──┤ H ├──■──»
          └───┘└───┘                                    └───┘     »
«     ┌──────────────────────────┐┌───┐┌───────────────────────┐ »
«q_0: ┤ Rz(-0.00564005000000000) ├┤ X ├┤ Rz(0.198968710000000) ├─»
«     └──────────────────────────┘└─┬─┘├───────────────────────┴┐»
«q_1: ──────────────────────────────■──┤ Rz(-0.198968700000000) ├»
«                                      └────────────────────────┘»
«     ┌───────────────────────┐ ┌───┐┌──────────────────────────┐┌───┐┌───┐»
«q_0: ┤ Rz(0.198968710000000) ├─┤ X ├┤ Rz(-0.00564005000000000) ├┤ X ├┤ H ├»
«     ├───────────────────────┴┐└─┬─┘└──────────────────────────┘└─┬─┘├───┤»
«q_1: ┤ Rz(-0.198968700000000) ├──■────────────────────────────────■──┤ H ├»
«     └────────────────────────┘                                      └───┘»
«     ┌───┐┌────────────────────────┐┌───┐┌───┐
«q_0: ┤ X ├┤ Rz(0.0904655950000000) ├┤ X ├┤ H ├
«     └─┬─┘└────────────────────────┘└─┬─┘├───┤
«q_1: ──■──────────────────────────────■──┤ H ├
«                                         └───┘

Expectations

Expectations are converters that enable the computation of expectation values of observables. They traverse an Operator tree, replacing OperatorStateFns (observables) with equivalent instructions which are more amenable to computation on quantum or classical hardware. For example, if we want to measure the expectation value of an Operator o expressed as a sum of Paulis with respect to some state function, but can only access diagonal measurements on quantum hardware, we can create an observable ~StateFn(o) and use a PauliExpectation to convert it to a diagonal measurement and circuit pre-rotations to append to the state.

Another interesting Expectation is the AerPauliExpectation, which converts the observable into a CircuitStateFn containing a special expectation snapshot instruction which Aer can execute natively with high performance.

# Note that XX was the only non-diagonal measurement in our H2 Observable
print(PauliExpectation(group_paulis=False).convert(h2_measurement))

SummedOp([
  ComposedOp([
    OperatorMeasurement(-1.0523732 * II),
    II
  ]),
  ComposedOp([
    OperatorMeasurement(0.39793742 * IZ),
    II
  ]),
  ComposedOp([
    OperatorMeasurement(-0.3979374 * ZI),
    II
  ]),
  ComposedOp([
    OperatorMeasurement(-0.0112801 * ZZ),
    II
  ]),
  ComposedOp([
    OperatorMeasurement(0.18093119 * ZZ),
         ┌───┐
    q_0: ┤ H ├
         ├───┤
    q_1: ┤ H ├
         └───┘
  ])
])

By default group_paulis=True, which will use the AbelianGrouper to convert the SummedOp into groups of mutually qubit-wise commuting Paulis. This reduces circuit execution overhead, as each group can share the same circuit execution.

print(PauliExpectation().convert(h2_measurement))

SummedOp([
  ComposedOp([
    OperatorMeasurement(0.18093119 * ZZ
    - 1.0523732 * II),
         ┌───┐
    q_0: ┤ H ├
         ├───┤
    q_1: ┤ H ├
         └───┘
  ]),
  ComposedOp([
    OperatorMeasurement(0.39793742 * IZ
    - 0.3979374 * ZI
    - 0.0112801 * ZZ),
    II
  ])
])

Note that converters act recursively, that is, they traverse an expression applying their action only where possible. So we can just convert our full evolution and measurement expression. We could have equivalently composed the converted h2_measurement with our evolution CircuitStateFn. We proceed by applying the conversion on the entire expression.

diagonalized_meas_op = PauliExpectation().convert(trotterized_op)
print(diagonalized_meas_op)

SummedOp([
  ComposedOp([
    OperatorMeasurement(0.18093119 * ZZ
    - 1.0523732 * II),
    CircuitStateFn(
    global phase: 1.0523732*θ
         ┌───┐     ┌───┐┌───┐┌──────────────────┐┌───┐┌───┐┌───┐»
    q_0: ┤ H ├──■──┤ H ├┤ X ├┤ Rz(0.18093119*θ) ├┤ X ├┤ H ├┤ X ├»
         └───┘┌─┴─┐├───┤└─┬─┘└──────────────────┘└─┬─┘├───┤└─┬─┘»
    q_1: ─────┤ X ├┤ H ├──■────────────────────────■──┤ H ├──■──»
              └───┘└───┘                              └───┘     »
    «     ┌──────────────────┐┌───┐┌──────────────────┐┌──────────────────┐┌───┐»
    «q_0: ┤ Rz(-0.0112801*θ) ├┤ X ├┤ Rz(0.39793742*θ) ├┤ Rz(0.39793742*θ) ├┤ X ├»
    «     └──────────────────┘└─┬─┘├──────────────────┤├──────────────────┤└─┬─┘»
    «q_1: ──────────────────────■──┤ Rz(-0.3979374*θ) ├┤ Rz(-0.3979374*θ) ├──■──»
    «                              └──────────────────┘└──────────────────┘     »
    «     ┌──────────────────┐┌───┐┌───┐┌───┐┌──────────────────┐┌───┐┌───┐┌───┐
    «q_0: ┤ Rz(-0.0112801*θ) ├┤ X ├┤ H ├┤ X ├┤ Rz(0.18093119*θ) ├┤ X ├┤ H ├┤ H ├
    «     └──────────────────┘└─┬─┘├───┤└─┬─┘└──────────────────┘└─┬─┘├───┤├───┤
    «q_1: ──────────────────────■──┤ H ├──■────────────────────────■──┤ H ├┤ H ├
    «                              └───┘                              └───┘└───┘
    )
  ]),
  ComposedOp([
    OperatorMeasurement(0.39793742 * IZ
    - 0.3979374 * ZI
    - 0.0112801 * ZZ),
    CircuitStateFn(
    global phase: 1.0523732*θ
         ┌───┐     ┌───┐┌───┐┌──────────────────┐┌───┐┌───┐┌───┐»
    q_0: ┤ H ├──■──┤ H ├┤ X ├┤ Rz(0.18093119*θ) ├┤ X ├┤ H ├┤ X ├»
         └───┘┌─┴─┐├───┤└─┬─┘└──────────────────┘└─┬─┘├───┤└─┬─┘»
    q_1: ─────┤ X ├┤ H ├──■────────────────────────■──┤ H ├──■──»
              └───┘└───┘                              └───┘     »
    «     ┌──────────────────┐┌───┐┌──────────────────┐┌──────────────────┐┌───┐»
    «q_0: ┤ Rz(-0.0112801*θ) ├┤ X ├┤ Rz(0.39793742*θ) ├┤ Rz(0.39793742*θ) ├┤ X ├»
    «     └──────────────────┘└─┬─┘├──────────────────┤├──────────────────┤└─┬─┘»
    «q_1: ──────────────────────■──┤ Rz(-0.3979374*θ) ├┤ Rz(-0.3979374*θ) ├──■──»
    «                              └──────────────────┘└──────────────────┘     »
    «     ┌──────────────────┐┌───┐┌───┐┌───┐┌──────────────────┐┌───┐┌───┐
    «q_0: ┤ Rz(-0.0112801*θ) ├┤ X ├┤ H ├┤ X ├┤ Rz(0.18093119*θ) ├┤ X ├┤ H ├
    «     └──────────────────┘└─┬─┘├───┤└─┬─┘└──────────────────┘└─┬─┘├───┤
    «q_1: ──────────────────────■──┤ H ├──■────────────────────────■──┤ H ├
    «                              └───┘                              └───┘
    )
  ])
])

Now we bind multiple parameter values into a ListOp, followed by eval to evaluate the entire expression. We could have used eval earlier if we bound earlier, but it would not be efficient. Here, eval will convert our CircuitStateFns to VectorStateFns through simulation internally.

evo_time_points = list(range(8))
h2_trotter_expectations = diagonalized_meas_op.bind_parameters({evo_time: evo_time_points})

Here are the expectation values \(\langle \Phi_+| e^{iHt} H e^{-iHt} |\Phi_+\rangle\) corresponding to the different values of the parameter.

h2_trotter_expectations.eval()

array([-0.88272211+0.0e+00j, -0.88272211+0.0e+00j, -0.88272211+0.0e+00j,
       -0.88272211+0.0e+00j, -0.88272211+0.0e+00j, -0.88272211+2.8e-17j,
       -0.88272211+5.6e-17j, -0.88272211+0.0e+00j])

Executing CircuitStateFns with the CircuitSampler

The CircuitSampler traverses an Operator and converts any CircuitStateFns into approximations of the resulting state function by a DictStateFn or VectorStateFn using a quantum backend. Note that in order to approximate the value of the CircuitStateFn, it must 1) send the state function through a depolarizing channel, which will destroy all phase information and 2) replace the sampled frequencies with square roots of the frequency, rather than the raw probability of sampling (which would be the equivalent of sampling the square of the state function, per the Born rule.)

sampler = CircuitSampler(backend=Aer.get_backend('aer_simulator'))
# sampler.quantum_instance.run_config.shots = 1000
sampled_trotter_exp_op = sampler.convert(h2_trotter_expectations)
sampled_trotter_energies = sampled_trotter_exp_op.eval()
print('Sampled Trotterized energies:\n {}'.format(np.real(sampled_trotter_energies)))

Sampled Trotterized energies:
 [-0.88272211 -0.88272211 -0.88272211 -0.88272211 -0.88272211 -0.88272211
 -0.88272211 -0.88272211]

Note again that the circuits are replaced by dicts with square roots of the circuit sampling probabilities. Take a look at one sub-expression before and after the conversion:

print('Before:\n')
print(h2_trotter_expectations.reduce()[0][0])
print('\nAfter:\n')
print(sampled_trotter_exp_op[0][0])

Before:

ComposedOp([
  OperatorMeasurement(0.18093119 * ZZ
  - 1.0523732 * II),
  CircuitStateFn(
       ┌───┐     ┌───┐┌───┐┌───────┐┌───┐┌───┐┌───┐┌───────┐┌───┐┌───────┐»
  q_0: ┤ H ├──■──┤ H ├┤ X ├┤ Rz(0) ├┤ X ├┤ H ├┤ X ├┤ Rz(0) ├┤ X ├┤ Rz(0) ├»
       └───┘┌─┴─┐├───┤└─┬─┘└───────┘└─┬─┘├───┤└─┬─┘└───────┘└─┬─┘├───────┤»
  q_1: ─────┤ X ├┤ H ├──■─────────────■──┤ H ├──■─────────────■──┤ Rz(0) ├»
            └───┘└───┘                   └───┘                   └───────┘»
  «     ┌───────┐┌───┐┌───────┐┌───┐┌───┐┌───┐┌───────┐┌───┐┌───┐┌───┐
  «q_0: ┤ Rz(0) ├┤ X ├┤ Rz(0) ├┤ X ├┤ H ├┤ X ├┤ Rz(0) ├┤ X ├┤ H ├┤ H ├
  «     ├───────┤└─┬─┘└───────┘└─┬─┘├───┤└─┬─┘└───────┘└─┬─┘├───┤├───┤
  «q_1: ┤ Rz(0) ├──■─────────────■──┤ H ├──■─────────────■──┤ H ├┤ H ├
  «     └───────┘                   └───┘                   └───┘└───┘
  )
])

After:

ComposedOp([
  OperatorMeasurement(0.18093119 * ZZ
  - 1.0523732 * II),
  DictStateFn({'11': 0.7036456405748563, '00': 0.7105510625563796})
])

/home/computertreker/git/qiskit/qiskit/.tox/docs/lib/python3.7/site-packages/sympy/core/expr.py:3951: SymPyDeprecationWarning:

expr_free_symbols method has been deprecated since SymPy 1.9. See
https://github.com/sympy/sympy/issues/21494 for more info.

  deprecated_since_version="1.9").warn()

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

/home/computertreker/git/qiskit/qiskit/.tox/docs/lib/python3.7/site-packages/qiskit/aqua/__init__.py:86: DeprecationWarning: The package qiskit.aqua is deprecated. It was moved/refactored to qiskit-terra For more information see <https://github.com/Qiskit/qiskit-aqua/blob/main/README.md#migration-guide>
  warn_package('aqua', 'qiskit-terra')

Version Information

Qiskit Software	Version
qiskit-terra	0.18.2
qiskit-aer	0.8.2
qiskit-ignis	0.6.0
qiskit-ibmq-provider	0.16.0
qiskit-aqua	0.9.5
qiskit	0.29.1
qiskit-nature	0.2.2
qiskit-finance	0.3.0
qiskit-optimization	0.2.3
qiskit-machine-learning	0.2.1
System information
Python	3.7.12 (default, Nov 22 2021, 14:57:10) [GCC 11.1.0]
OS	Linux
CPUs	32
Memory (Gb)	125.71650314331055
Tue Jan 04 11:16:44 2022 EST

This code is a part of Qiskit

© Copyright IBM 2017, 2022.

This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.

Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.