Variational algorithms

Recently, variational quantum algorithms are actively studied, where optimal values of parameters in parametric quantum circuits are searched. In this section, we first describe estimation of gradient of operator expectations, and then see how to construct one of the variational algorithms, variational quantum eigensolver (VQE), using the gradient.

Prerequisite

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

!pip install "quri-parts[qulacs]"

Gradient of operator expectation

Variational algorithms often involves minimization of a certain cost function, which is defined as an expectation value of an operator for a certain parametric quantum state: \(f(\boldsymbol{\theta}) = \langle O\rangle_\boldsymbol{\theta} = \langle \psi(\boldsymbol{\theta})|O|\psi(\boldsymbol{\theta})\rangle\), where \(O\) is the operator and \(\psi(\boldsymbol{\theta})\) is a parametric state with parameters \(\boldsymbol{\theta} = \theta_0, \ldots, \theta_{m-1}\). In such minimization, gradient of the cost function is often used, defined as:

\begin{equation} \nabla_\boldsymbol{\theta} f(\boldsymbol{\theta}) = \left( \dfrac{\partial \langle O\rangle_\boldsymbol{\theta}}{\partial \theta_0}, \ldots, \dfrac{\partial \langle O\rangle_\boldsymbol{\theta}}{\partial \theta_{m-1}} \right) \end{equation}

Below we will see two methods to estimate such gradient: numerical gradient and parameter shift rule.

Let’s prepare a target operator and a parametric state here:

from quri_parts.core.operator import Operator, pauli_label

op = Operator({
    pauli_label("X0 Y1"): 0.5 + 0.5j,
    pauli_label("Z0 X1"): 0.2,
})

from math import pi
from quri_parts.circuit import LinearMappedUnboundParametricQuantumCircuit, CONST

param_circuit = LinearMappedUnboundParametricQuantumCircuit(2)
param_circuit.add_H_gate(0)
param_circuit.add_CNOT_gate(0, 1)

theta, phi = param_circuit.add_parameters("theta", "phi")
param_circuit.add_ParametricRX_gate(0, {theta: 1/2, phi: 1/3, CONST: pi/2})
param_circuit.add_ParametricRZ_gate(1, {theta: 1/3, phi: -1/2, CONST: -pi/2})

from quri_parts.core.state import ParametricCircuitQuantumState

param_state = ParametricCircuitQuantumState(2, param_circuit)

Numerical gradient

QURI Parts defines quri_parts.core.estimator.GradientEstimator interface to estimate the gradient. A simple gradient estimator using numerical differentiation is provided:

from quri_parts.core.estimator.gradient import create_numerical_gradient_estimator
from quri_parts.qulacs.estimator import create_qulacs_vector_concurrent_parametric_estimator

qulacs_concurrent_parametric_estimator = create_qulacs_vector_concurrent_parametric_estimator()
gradient_estimator = create_numerical_gradient_estimator(
    qulacs_concurrent_parametric_estimator,
    delta=1e-4,
)

gradient = gradient_estimator(op, param_state, [0.2, 0.3])
print("Estimated gradient:", gradient.values)

Estimated gradient: [(0.0004866565750383245-0.013872819366045341j), (0.042165661391369014+0.020809229047680233j)]

When creating a gradient estimator, two arguments are given. The first one is a ConcurrentParametricQuantumEstimator, which is used to estimate the expectation values at slightly shifted parameter values. The second one delta is the step size for numerical differentiation. The gradient estimator is invoked with the operator, the parametric state and the parameter values at which the gradient is evaluated, and returns the estimated gradient values.

Gradient evaluation by parameter shift rule

Parameter shift rule is a method to evaluate derivatives of an expectation value of an operator for a state generated by a parametric circuit [1]_. QURI Parts provides a way to estimate the gradient using parameter shift rule. Here, it is assumed that any parametric gate in a parametric circuit is defined as \(\exp(-i\theta P/2)\), where \(P\) is a Pauli product. All parametric gates defined in QURI Parts satisfy this condition. We also assume that all the gate parameters depend on the circuit parameters (at most) linearly.

[1]

Mitarai, K. and Negoro, M. and Kitagawa, M. and Fujii, K., Phys. Rev. A 98, 032309 (2018). arXiv:1803.00745.

from quri_parts.core.estimator.gradient import create_parameter_shift_gradient_estimator
from quri_parts.qulacs.estimator import create_qulacs_vector_concurrent_parametric_estimator

qulacs_concurrent_parametric_estimator = create_qulacs_vector_concurrent_parametric_estimator()
gradient_estimator = create_parameter_shift_gradient_estimator(
    qulacs_concurrent_parametric_estimator,
)

gradient = gradient_estimator(op, param_state, [0.2, 0.3])
print("Estimated gradient:", gradient.values)

Estimated gradient: [(0.0004866565766964738-0.013872819366718317j), (0.04216566140053679+0.020809229050077496j)]

If you would like to know how the gradient evaluation by parameter shift rule works in detail, read the following explanation. You can skip it and go to the next topic Variational quantum eigensolver.

Explanation of how gradient evaluation by parameter shift rule works

When evaluating the gradient with parameter shift rule, parameters of each parametric gates need to be shifted independently, even if they depend on the same circuit parameters. It is also necessary to compute derivative of each gate parameter with respect to the circuit parameters so that we can use chain rule of differentiation. Therefore we need the followings:

• The parametric circuit where each gate parameters are treated as independent (UnboundParametricQuantumCircuit in QURI Parts).

• Parameter shifts for each gate parameters for each circuit parameters.

• Differential coefficents corresponding to each parameter shifts.

The following function computes and returns the above:

from quri_parts.circuit.parameter_shift import ShiftedParameters
from quri_parts.core.state import ParametricCircuitQuantumState

def get_raw_param_state_and_shifted_parameters(state, params):
    param_mapping = state.parametric_circuit.param_mapping
    raw_circuit = state.parametric_circuit.primitive_circuit()
    parameter_shift = ShiftedParameters(param_mapping)
    derivatives = parameter_shift.get_derivatives()
    shifted_parameters = [
        d.get_shifted_parameters_and_coef(params) for d in derivatives
    ]

    raw_param_state = ParametricCircuitQuantumState(state.qubit_count, raw_circuit)

    return raw_param_state, shifted_parameters

# Example
raw_state, shifted_params_and_coefs = get_raw_param_state_and_shifted_parameters(
    param_state, [0.2, 0.3]
)

for i, params_and_coefs in enumerate(shifted_params_and_coefs):
    print(f"Parameter shifts for circuit parameter {i}:")
    for p, c in params_and_coefs:
        print(f"  gate params: {p}, coefficient: {c}")

Parameter shifts for circuit parameter 0:
  gate params: (3.3415926535897933, -1.6541296601282298), coefficient: 0.25
  gate params: (1.7707963267948965, -3.224925986923126), coefficient: -0.16666666666666666
  gate params: (1.7707963267948965, -0.08333333333333326), coefficient: 0.16666666666666666
  gate params: (0.19999999999999996, -1.6541296601282298), coefficient: -0.25
Parameter shifts for circuit parameter 1:
  gate params: (3.3415926535897933, -1.6541296601282298), coefficient: 0.16666666666666666
  gate params: (1.7707963267948965, -3.224925986923126), coefficient: 0.25
  gate params: (1.7707963267948965, -0.08333333333333326), coefficient: -0.25
  gate params: (0.19999999999999996, -1.6541296601282298), coefficient: -0.16666666666666666

We then obtain the gradient by 1) estimating the expectation value of the operator for each shifted gate parameters, and 2) sum up them with the corresponding coefficients multiplied. This can be done as follows:

from quri_parts.qulacs.estimator import create_qulacs_vector_concurrent_parametric_estimator

def get_parameter_shift_gradient(op, raw_state, shifted_params_and_coefs):
    # Collect gate parameters to be evaluated
    gate_params = set()
    for params_and_coefs in shifted_params_and_coefs:
        for p, _ in params_and_coefs:
            gate_params.add(p)
    gate_params_list = list(gate_params)

    # Prepare a parametric estimator
    estimator = create_qulacs_vector_concurrent_parametric_estimator()

    # Estimate the expectation values
    estimates = estimator(op, raw_state, gate_params_list)
    estimates_dict = dict(zip(gate_params_list, estimates))

    # Sum up the expectation values with the coefficients multiplied
    gradient = []
    for params_and_coefs in shifted_params_and_coefs:
        g = 0.0
        for p, c in params_and_coefs:
            g += estimates_dict[p].value * c
        gradient.append(g)

    return gradient

# Example
gradient = get_parameter_shift_gradient(op, raw_state, shifted_params_and_coefs)
print("Estimated gradient:", gradient)

Estimated gradient: [(0.00048665657669647033-0.01387281936671833j), (0.04216566140053679+0.020809229050077496j)]

These functions can be combined to create a GradientEstimator as follows:

from collections.abc import Sequence
from dataclasses import dataclass

# This is a return type of GradientEstimator
@dataclass
class _Estimates:
    values: Sequence[complex]
    error_matrix = None

def parameter_shift_gradient_estimator(op, state, params):
    raw_state, shifted_params_and_coefs = get_raw_param_state_and_shifted_parameters(state, params)
    gradient = get_parameter_shift_gradient(op, raw_state, shifted_params_and_coefs)
    return _Estimates(gradient)

# Example
gradient = parameter_shift_gradient_estimator(op, param_state, [0.2, 0.3])
print("Estimated gradient:", gradient.values)

Estimated gradient: [(0.00048665657669647033-0.01387281936671833j), (0.04216566140053679+0.020809229050077496j)]

Variational quantum eigensolver

Variational quantum eigensolver (VQE) is a method to optimize an expectation value of an operator (e.g. energy of a molecule) over parametrized quantum states. There are two major components in VQE:

• Ansatz: A parametric quantum circuit which generates the parametrized quantum states subject to optimization

• Optimizer: A method to numerically optimize the expectation value of the operator

Ansatz

In context of VQE, ansatz refers to a parametric quantum circuit used for generating parametrized quantum states for which expectation values of the target operator is evaluated. You can define a (LinearMapped)UnboundParametricQuantumCircuit on your own, or use a well-known ansatz defined in quri_parts.algo.ansatz package. In this example we use a hardware-efficient ansatz [1]_:

[1]

Kandala, A., Mezzacapo, A., Temme, K. et al. Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets. Nature 549, 242–246 (2017).

from quri_parts.algo.ansatz import HardwareEfficient

hw_ansatz = HardwareEfficient(qubit_count=4, reps=3)

In order to evaluate the expectation value, the parametrized quantum state is necessary, which is obtained by applying the ansatz to a specific initial state. Here we use a computational basis state \(|0011\rangle\).

from quri_parts.core.state import ComputationalBasisState, ParametricCircuitQuantumState, apply_circuit

cb_state = ComputationalBasisState(4, bits=0b0011)
parametric_state = apply_circuit(hw_ansatz, cb_state)

Optimizer

An optimizer searches optimal parameters that minimize a given cost function. In context of VQE, the cost function is the expectation value of the target operator. Some optimizers use only the cost function itself, while others use gradient of the cost function for efficient optimization. You can use optimizers provided by libraries such as scipy.optimize, or ones provided in quri_parts.algo.optimizer package. In this example we use Adam [1]_, which uses the gradient.

[1]

Diederik P. Kingma, Jimmy Ba, Adam: A Method for Stochastic Optimization. arXiv:1412.6980 (2014)

from quri_parts.algo.optimizer import Adam

# You can pass optional parameters. See the reference for details
adam_optimizer = Adam()

Running VQE

We first define a target operator, whose expectation value is subject to the optimization:

from quri_parts.core.operator import Operator, pauli_label, PAULI_IDENTITY

# This is Jordan-Wigner transformed Hamiltonian of a hydrogen molecule
hamiltonian = Operator({
    PAULI_IDENTITY: 0.03775110394645542,
    pauli_label("Z0"): 0.18601648886230593,
    pauli_label("Z1"): 0.18601648886230593,
    pauli_label("Z2"): -0.2694169314163197,
    pauli_label("Z3"): -0.2694169314163197,
    pauli_label("Z0 Z1"): 0.172976101307451,
    pauli_label("Z0 Z2"): 0.12584136558006326,
    pauli_label("Z0 Z3"): 0.16992097848261506,
    pauli_label("Z1 Z2"): 0.16992097848261506,
    pauli_label("Z1 Z3"): 0.12584136558006326,
    pauli_label("Z2 Z3"): 0.17866777775953396,
    pauli_label("X0 X1 Y2 Y3"): -0.044079612902551774,
    pauli_label("X0 Y1 Y2 X3"): 0.044079612902551774,
    pauli_label("Y0 X1 X2 Y3"): 0.044079612902551774,
    pauli_label("Y0 Y1 X2 X3"): -0.044079612902551774,
})

Using this operator and the parametric state prepared above, we can define the cost function as a function of the circuit parameters:

from quri_parts.qulacs.estimator import create_qulacs_vector_parametric_estimator

estimator = create_qulacs_vector_parametric_estimator()

def cost_fn(param_values):
    estimate = estimator(hamiltonian, parametric_state, param_values)
    return estimate.value.real

We also define gradient of the cost function using numerical gradient:

import numpy as np
from quri_parts.core.estimator.gradient import create_numerical_gradient_estimator
from quri_parts.qulacs.estimator import create_qulacs_vector_concurrent_parametric_estimator

qulacs_concurrent_parametric_estimator = create_qulacs_vector_concurrent_parametric_estimator()
gradient_estimator = create_numerical_gradient_estimator(
    qulacs_concurrent_parametric_estimator,
    delta=1e-4,
)

def grad_fn(param_values):
    estimate = gradient_estimator(hamiltonian, parametric_state, param_values)
    return np.asarray([g.real for g in estimate.values])

Then we can run VQE with a QURI Parts optimizer:

from quri_parts.algo.optimizer import OptimizerStatus

def vqe(operator, init_params, cost_fn, grad_fn, optimizer):
    opt_state = optimizer.get_init_state(init_params)
    while True:
        opt_state = optimizer.step(opt_state, cost_fn, grad_fn)
        if opt_state.status == OptimizerStatus.FAILED:
            print("Optimizer failed")
            break
        if opt_state.status == OptimizerStatus.CONVERGED:
            print("Optimizer converged")
            break
    return opt_state

init_params = [0.1] * hw_ansatz.parameter_count
result = vqe(hamiltonian, init_params, cost_fn, grad_fn, adam_optimizer)
print("Optimized value:", result.cost)
print("Optimized parameter:", result.params)
print("Iterations:", result.niter)
print("Cost function calls:", result.funcalls)
print("Gradient function calls:", result.gradcalls)

Optimizer converged
Optimized value: -1.1119813406307852
Optimized parameter: [ 5.47178291e-02  8.40762138e-02  5.12253344e-02  8.19750289e-02
 -9.72099777e-03 -1.16141829e-01 -3.06727509e-03  9.66792838e-01
  1.27323903e-01  1.04790835e-01  1.27097746e-01  9.40512680e-02
 -1.60419268e-02  9.92326569e-01 -3.35897820e-02  9.91027220e-01
  6.44048147e-02  2.49943780e-04  6.43611653e-02 -5.72092116e-03
 -1.48640084e-02 -1.16555432e-01 -3.59503991e-02  9.79005522e-01
  1.67652637e-02 -2.35033521e-01  1.34115103e-02 -2.24492480e-01
 -2.91851961e-02  4.35033663e-01 -3.52284766e-03  4.24493109e-01]
Iterations: 24
Cost function calls: 25
Gradient function calls: 24

You can also run VQE with a SciPy optimizer:

from scipy.optimize import minimize

def vqe_scipy(operator, init_params, cost_fn, grad_fn, method, options):
    return minimize(cost_fn, init_params, jac=grad_fn, method=method, options=options)

init_params = [0.1] * hw_ansatz.parameter_count
bfgs_options = {
    "gtol": 1e-6,
}
result = vqe_scipy(hamiltonian, init_params, cost_fn, grad_fn, "BFGS", bfgs_options)
print(result.message)
print("Optimized value:", result.fun)
print("Optimized parameter:", result.x)
print("Iterations:", result.nit)
print("Cost function calls:", result.nfev)
print("Gradient function calls:", result.njev)

Optimization terminated successfully.
Optimized value: -1.129904784280104
Optimized parameter: [ 7.82483153e-04  4.29969941e-02  6.61286966e-01  2.10598472e-03
  3.12868377e-01 -4.24005711e-02 -1.39244497e+00 -1.19450656e-03
  3.36566899e-01  8.08495521e-05  6.57323966e-01 -2.87283611e-01
  6.78503238e-01  1.15585829e-01  2.19568571e+00 -2.24222650e-03
  1.57043690e+00 -6.09357373e-08 -4.37512017e-04  1.77385571e-01
  1.79573603e-01 -1.42632685e-01 -2.29170046e-01  1.36120809e-02
  1.23449374e+00  1.10763438e-01 -5.92928055e-04  1.10784204e-01
 -5.21553526e-01  8.92367689e-02 -6.31244124e-01  8.92156968e-02]
Iterations: 174
Cost function calls: 181
Gradient function calls: 181