Playing with cuQuantum (5) — VQE Part 2

Purpose

In Playing with cuQuantum (4) — VQE, I performed VQE in a rather forced manner. Since I want to do something a bit more elegant this time, I'd like to use the cuQuantum API to calculate expectation values.

I also plan to run and compare it with the previous brute-force version.

Problem

I want to minimize the expectation value calculation for $\mathcal{H} = Z \otimes X$. Essentially, I want to approximate $\ket{\psi} = \ket{1} \otimes \ket{+} = \frac{1}{\sqrt{2}}(\ket{10} + \ket{11})$ or $\ket{\psi^\prime} = \ket{0} \otimes \ket{-} = \frac{1}{\sqrt{2}}(\ket{00} - \ket{01})$ using a PQC (Parameterized Quantum Circuit).

Brute-force Version

Import the necessary modules:

from qiskit import QuantumCircuit, Aer
from qiskit.circuit.library import TwoLocal
from qiskit.quantum_info import SparsePauliOp
from qiskit.quantum_info import Statevector
from qiskit.visualization import plot_histogram
from cuquantum import CircuitToEinsum, contract
import numpy as np
import cupy as cp
from scipy.optimize import minimize

Then, perform the forced implementation just like last time:

def create_ansatz(num_qubits):
    ansatz = TwoLocal(num_qubits, 'ry', 'cz')
    return ansatz, ansatz.inverse()

hamiltonian = SparsePauliOp.from_list([
    ('ZX', 1.0)
])
hamiltonian = cp.array(hamiltonian.to_matrix().reshape(2, 2, 2, 2))

num_qubits = 2
ansatz, ansatz_dagger = create_ansatz(num_qubits)

def expectation(theta, ansatz, ansatz_dagger, hamiltonian):
    assert len(theta) == 4 * ansatz.num_qubits
    ansatz = ansatz.bind_parameters(theta)
    ansatz_dagger = ansatz_dagger.bind_parameters(theta)

    converter = CircuitToEinsum(ansatz)
    expr, operands = converter.state_vector()
    vec = contract(expr, *operands)

    converter = CircuitToEinsum(ansatz_dagger)
    expr, operands = converter.amplitude('0' * ansatz.num_qubits)
    # 'a,b,cb,da,efcd...->' ==> 'cb,da,efcd...->ab'
    out, expr = expr[:ansatz.num_qubits*2].replace(',', ''), \
                expr[ansatz.num_qubits*2:]
    expr += out
    vec_dagger = contract(expr, *(operands[ansatz.num_qubits:]))

    val = contract('ab,cdba,dc->', vec, hamiltonian, vec_dagger)
    return float(val.real)

Experiment:

theta = np.random.random(4 * num_qubits)
args = (ansatz, ansatz_dagger, hamiltonian)

result = minimize(expectation, theta, args=args, method='Powell')
print(f'optimal_value: {result.fun}')
print(f'x: {result.x}')

optimal_value: -1.0000000000000007
x: [-2.16629773 0.87642794 0.48274178 0.58879033 0.17536848 0.69861945
0.40712598 0.19401039]

Result:

qc = TwoLocal(2, 'ry', 'cz').decompose()
qc = qc.bind_parameters(result.x)
qc.measure_all()
backend = Aer.get_backend('aer_simulator')
counts = backend.run(qc).result().get_counts()
plot_histogram(counts)

It might seem like a somewhat ambiguous result. However, I've included it here because it is actually quite interesting upon closer inspection. This roughly means that:

$$\begin{align*} \ket{\psi} = \sqrt{\frac{3}{5}} \frac{\ket{00} - \ket{01}}{\sqrt{2}} + \sqrt{\frac{2}{5}} \frac{\ket{10} + \ket{11}}{\sqrt{2}} = \sqrt{\frac{3}{5}} \ket{0-} + \sqrt{\frac{2}{5}} \ket{1+} \end{align*}$$

is being produced. Let's verify this later. To double-check, calculating it gives:

$$\newcommand \parenthetical[1] { \left( #1 \right) } \begin{align*} \braket{\psi | Z \!\otimes\! X | \psi} &= \parenthetical{ \sqrt{\frac{3}{5}} \bra{0-} + \sqrt{\frac{2}{5}} \bra{1+} } \! (Z \!\otimes\! X) \! \parenthetical{ \sqrt{\frac{3}{5}} \ket{0-} + \sqrt{\frac{2}{5}} \ket{1+} } \\ & = \parenthetical{ \sqrt{\frac{3}{5}} \bra{0-} + \sqrt{\frac{2}{5}} \bra{1+} } \! \parenthetical{ -\sqrt{\frac{3}{5}} \ket{0-} - \sqrt{\frac{2}{5}} \ket{1+} } = - \frac{3}{5} - \frac{2}{5} = -1 \end{align*}$$

That's how it works. That said, since phase information is lost during measurement, let's determine it using the state vector.

qc.remove_final_measurements()
sv = Statevector(qc)
sv.draw('latex')

$$\begin{align*} 0.5492966902 \ket{00} - 0.5492966858 \ket{01} + 0.4452787292 \ket{10} + 0.4452787335 \ket{11} \end{align*}$$

In short, since the eigenvalue -1 of the Hamiltonian is degenerate, a superposition of eigenvectors was obtained.

Actually, these proportions fluctuated and were unstable with each VQE run. In the following "cuQuantum API Version" section, I obtained the result I "initially expected to see." Since the results fluctuate there as well with each VQE execution, it just happened to turn out this way this time.

cuQuantum API Version

Expectation value calculation implementation:

If the Hamiltonian is a Pauli operator, the expectation method can be used, making the expectation value calculation significantly easier.

num_qubits = 2
ansatz, _ = create_ansatz(num_qubits)

def expectation_cutn(theta, ansatz, hamiltonian: str):
    assert len(theta) == 4 * ansatz.num_qubits
    ansatz = ansatz.bind_parameters(theta)

    converter = CircuitToEinsum(ansatz)
    expr, operands = converter.expectation(hamiltonian)
    val = contract(expr, *operands)

    return float(val.real)

Experiment:

theta = np.random.random(4 * num_qubits)
args = (ansatz, hamiltonian)

hamiltonian = 'XZ'

result = minimize(expectation_cutn, theta, args=args, method='Powell')
print(f'optimal_value: {result.fun}')
print(f'x: {result.x}')

optimal_value: -1.0000000000000004
x: [2.37477257 0.14751172 1.11527912 0.41963828 0.3659612 0.52473461
0.84868946 0.13640107]

Result:

qc = TwoLocal(2, 'ry', 'cz').decompose()
qc = qc.bind_parameters(result.x)
qc.measure_all()
backend = Aer.get_backend('aer_simulator')
counts = backend.run(qc).result().get_counts()
plot_histogram(counts)

This roughly means that:

$$\begin{align*} \ket{\psi} = \frac{1}{\sqrt{2}}(\ket{00} - \ket{01}) = \ket{0-} \end{align*}$$

is being produced. To double-check, the state vector is calculated as follows:

qc.remove_final_measurements()
sv = Statevector(qc)
sv.draw('latex')

$$\begin{align*} -0.7053123039 \ket{00} + 0.7053123164 \ket{01} - 0.0503442633 \ket{10} - 0.0503442692 \ket{11} \end{align*}$$

is obtained.

Summary

In both methods tried this time, the minimum eigenvalue of the Hamiltonian, -1, was obtained. However, the conclusion is that for Pauli operators, the expectation value calculation is much easier as it can be done in a single call to the cuQuantum API.

As for optimization methods, I'm not sure which one is better, but I felt that 'SLSQP' often resulted in a mixture of degenerate states. 'Powell' and 'COBYLA' seemed to produce a single eigenvector more clearly. This might not be essential and could just be a coincidence, though.

I also felt that caution is needed because the qubit ordering is reversed when defining a Hamiltonian with the Qiskit API's SparsePauliOp versus providing it to the cuQuantum API's expectation method.