Subspace-search Variational Quantum Eigensolver¶

Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved.

Overview¶

  • In this tutorial, we will show how to train a quantum neural network (QNN) through Paddle Quantum to find the entire energy spectrum of a quantum system.

  • First, import the following packages.

In [1]:
import numpy
from numpy import pi as PI
import paddle 
import paddle_quantum
from paddle import matmul
from paddle_quantum.ansatz import Circuit
from paddle_quantum.qinfo import random_pauli_str_generator, pauli_str_to_matrix
from paddle_quantum.linalg import dagger

Background¶

  • Variational Quantum Eigensolver (VQE) [1-3] is one of the most promising applications for near-term quantum computing. One of the its powerful versions is SSVQE [4], which can be used to find the ground state and the excited state of a physical system's Hamiltonian. Mathematically, one can interpret it as solving the eigenvalues and eigenvectors of a Hermitian matrix. The set of eigenvalues of the Hamiltonian is called the energy spectrum.
  • Next, we will use a brief example to demonstrate how to solve this problem by training a QNN, that is, to solve the energy spectrum of a given Hamiltonian HH.

Hamiltonian¶

  • For a specific molecule that needs to be analyzed, we need its geometry, charge, and spin multiplicity to obtain the Hamiltonian (in Pauli products form) describing the system. Specifically, through our built-in quantum chemistry toolkit, fermionic-to-qubit mapping technology can be used to output the qubit Hamiltonian.
  • As a simple demonstration of SSVQE, we provide a random 2-qubit Hamiltonian.
In [2]:
N = 2  # Number of qubits
SEED = 14  # Fixed random seed
In [5]:
# Generate random Hamiltonian represented by Pauli string
numpy.random.seed(SEED)
hamiltonian = random_pauli_str_generator(N, terms=10)
print("Random Hamiltonian in Pauli string format = \n", hamiltonian)

# Generate matrix representation of Hamiltonian
complex_dtype = paddle_quantum.get_dtype()
H = pauli_str_to_matrix(hamiltonian, N).astype(complex_dtype)

Random Hamiltonian in Pauli string format = 
 [[0.9152074787317819, 'x1,y0'], [-0.2717604556798945, 'z0'], [0.3628495008719168, 'x0'], [-0.5050129214094752, 'x1'], [-0.6971554357833791, 'y0,x1'], [0.8651151857574237, 'x0,y1'], [0.7409989105435002, 'y0'], [-0.39981603921243236, 'y0'], [0.06862640764702, 'z0'], [-0.7647553733438246, 'y1']]

Building a quantum neural network¶

  • To implement SSVQE, we first need to design a QNN U(θ)U(θ) (parameterized quantum circuit). In this tutorial, we provide a predefined universal quantum circuit template suitable for 2 qubits. Theoretically, this template has enough expressibility to simulate arbitrary 2-qubit unitary operation [5]. The specific implementation requires 3 CNOTCNOT gates plus 15 single-qubit rotation gates ∈{Ry,Rz}∈{Ry,Rz}.

  • One can randomly initialize the QNN parameters →θθ→ containing 15 parameters.

In [6]:
def U_theta(num_qubits: int) -> Circuit:
    """
    U_theta
    """
    # Initialize the quantum neural network according to the number of qubits/network width
    cir = Circuit(num_qubits)
    
    # Call the built-in quantum neural network template
    cir.universal_two_qubits([0, 1])

    # Return the circuit of the quantum neural network
    return cir

Training model and loss function¶

  • After setting up the Hamiltonian and the quantum neural network architecture, we will further define the parameters to be trained, the loss function and optimization methods. For a detailed inspection of the theory of SSVQE, please refer to the original paper [4].

  • By acting the quantum neural network U(θ)U(θ) on a set of orthogonal initial states (one can take the computational basis {|00⟩,|01⟩,|10⟩,|11⟩}{|00⟩,|01⟩,|10⟩,|11⟩}), we will get the output states {|ψ1(θ)⟩,|ψ2(θ)⟩,|ψ3(θ)⟩,|ψ4(θ)⟩}{|ψ1(θ)⟩,|ψ2(θ)⟩,|ψ3(θ)⟩,|ψ4(θ)⟩}.

  • Further, the loss function in the SSVQE model generally consists of expectation value of each output quantum state |ψk(θ)⟩|ψk(θ)⟩ given the Hamiltonian HH. More specifically, it's the weighted summation of the energy expectation value. In this example, the default weight vector is →w=[4,3,2,1]w→=[4,3,2,1].

  • The loss function is defined as:

L(θ)=2n∑k=1wk∗⟨ψk(θ)|H|ψk(θ)⟩.(1)(1)L(θ)=∑k=12nwk∗⟨ψk(θ)|H|ψk(θ)⟩.
In [7]:
class Net(paddle.nn.Layer):
    def __init__(self, num_qubits: int):
        super(Net, self).__init__()
        
        # Build quantum neural network
        self.cir = U_theta(num_qubits)

    # Define loss function and forward propagation mechanism
    def forward(self, H, N):
        
        # Calculate the loss function
        U = self.cir.unitary_matrix()
        loss_struct = paddle.real(matmul(matmul(dagger(U), H), U))

        # Enter the computational basis to calculate the expected value 
        # which is equivalent to taking the diagonal element of U^dagger*H*U
        loss_components = []
        for i in range(len(loss_struct)):
            loss_components.append(loss_struct[i][i])
        
        # Weighted summation of loss function
        loss = 0
        for i in range(len(loss_components)):
            weight = 4 - i
            loss += weight * loss_components[i]
        
        return loss, loss_components, self.cir

Hyper-parameters¶

Before training the quantum neural network, we also need to set up several hyper-parameters, mainly the learning rate LR, the number of iterations ITR. Here we set the learning rate to be LR = 0.3 and the number of iterations ITR = 100. One can adjust these hyper-parameters accordingly and check how they influence the training performance.

In [8]:
ITR = 100  # Set the total number of iterations of training
LR = 0.3  # Set the learning rate

Training process¶

  • After setting all the parameters of SSVQE model, we need to convert all the data into Tensor in the PaddlePaddle, and then train the quantum neural network.
  • We use Adam Optimizer in training, and one can also call other optimizers provided in PaddlePaddle.
In [10]:
paddle.seed(SEED)

# We need to convert numpy.ndarray to Tensor supported in Paddle
hamiltonian = paddle.to_tensor(H)

# Determine the parameter dimension of the network
net = Net(N)

# We use Adam optimizer for better performance
# One can change it to SGD or RMSprop.
opt = paddle.optimizer.Adam(learning_rate=LR, parameters=net.parameters())

# Optimization loop
for itr in range(1, ITR + 1):

    # Forward propagation calculates the loss function and returns the estimated energy spectrum
    loss, loss_components, cir = net(hamiltonian, N)

    # Under the dynamic graph mechanism, use back propagation to minimize the loss function
    loss.backward()
    opt.minimize(loss)
    opt.clear_grad()

    # Print training results
    if itr% 10 == 0:
        print('iter:', itr,'loss:','%.4f'% loss.numpy()[0])

iter: 10 loss: -4.5668
iter: 20 loss: -5.3998
iter: 30 loss: -5.6210
iter: 40 loss: -5.8872
iter: 50 loss: -5.9246
iter: 60 loss: -5.9471
iter: 70 loss: -5.9739
iter: 80 loss: -5.9833
iter: 90 loss: -5.9846
iter: 100 loss: -5.9848

Benchmarking¶

We have now completed the training of the quantum neural network, and we will verify the results by comparing them with theoretical values.

  • The theoretical Hamiltonian eigenvalues are solved by the linear algebra package in NumPy;
  • We compare the energy of each energy level obtained by training QNN with the theoretical value.
  • It can be seen that the training output is very close to the exact value.
In [11]:
def output_ordinalvalue(num):
    r"""
    Convert to ordinal value

    Args:
        num (int): input number

    Return:
        (str): output ordinal value
    """
    if num == 1:
        return str(num) + "st"
    elif num == 2:
        return str(num) + "nd"
    elif num == 3:
        return str(num) + "rd"
    else:
        return str(num) + 'th'

for i in range(len(loss_components)):
    if i == 0:
        print('The estimated ground state energy is: ', loss_components[i].numpy())
        print('The theoretical ground state energy is: ', numpy.linalg.eigh(H)[0][i])
    else:
        print('The estimated {} excited state energy is: {}'.format(
            output_ordinalvalue(i), loss_components[i].numpy())
        )
        print('The theoretical {} excited state energy is: {}'.format(
            output_ordinalvalue(i), numpy.linalg.eigh(H)[0][i])
        )

The estimated ground state energy is:  [-2.1876235]
The theoretical ground state energy is:  -2.187902
The estimated 1st excited state energy is: [-0.13721023]
The theoretical 1st excited state energy is: -0.13704127073287964
The estimated 2nd excited state energy is: [0.85251486]
The theoretical 2nd excited state energy is: 0.8523274064064026
The estimated 3rd excited state energy is: [1.4723194]
The theoretical 3rd excited state energy is: 1.4726158380508423

References¶

[1] Peruzzo, A. et al. A variational eigenvalue solver on a photonic quantum processor. Nat. Commun. 5, 4213 (2014).

[2] McArdle, S., Endo, S., Aspuru-Guzik, A., Benjamin, S. C. & Yuan, X. Quantum computational chemistry. Rev. Mod. Phys. 92, 015003 (2020).

[3] Cao, Y. et al. Quantum chemistry in the age of quantum computing. Chem. Rev. 119, 10856–10915 (2019).

[4] Nakanishi, K. M., Mitarai, K. & Fujii, K. Subspace-search variational quantum eigensolver for excited states. Phys. Rev. Res. 1, 033062 (2019).

[5] Vatan, F. & Williams, C. Optimal quantum circuits for general two-qubit gates. Phys. Rev. A 69, 032315 (2004).