import time
import numpy as np
from matplotlib import pyplot as plt

import paddle
import paddle_quantum
from paddle_quantum.ansatz import Circuit
from paddle_quantum.hamiltonian import Hamiltonian
from paddle_quantum.loss import ExpecVal
from paddle_quantum.state import State
from paddle_quantum.qinfo import pauli_str_to_matrix, schmidt_decompose
import warnings

warnings.filterwarnings("ignore")

N = 10      # Number of qubits
SEED = 16   # Fix a random seed
ITR = 100   # Set the number of learning iterations
LR = 0.1    # Set the learning rate
D = 3       # Set the depth for QNN

def get_ground_state_info(H):

    # Calculate the eigenvalues and eigenvectors of H
    vals, vecs = paddle.linalg.eigh(H)
    # Retrieve the ground state
    ground_state = paddle_quantum.State(vecs[:, 0])
    # Retrieve the ground state energy
    ground_state_energy = vals.tolist()[0]
    print(f'The ground state energy is {ground_state_energy:.5f} Ha.')
    # Run Schmidt decomposition on the ground state.
    l, _, _ = schmidt_decompose(ground_state)
    print(f'Schmidt rank of the ground state is {l.size}.')

    return ground_state_energy

# Fix a random seed
np.random.seed(SEED)

# Hard code a random a Hamiltonian
coefs = [-0.8886258, 0.453882]
pauli_str = ['x0,z1,z2,z4,x5,y6,y7,x8,x9', 'y0,x1,x2,x3,y4,x5,z6,z7,y8,x9']
pauli_str_A = ['x0,z1,z2,z4', 'y0,x1,x2,x3,y4']     # pauli substring for system A
pauli_str_B = ['x0,y1,y2,x3,x4', 'x0,z1,z2,y3,x4']  # pauli substring for system B

# Convert relavent object into Tensor form
H = Hamiltonian(list(zip(coefs, pauli_str)))
H_mtr = paddle.to_tensor(H.construct_h_matrix())
coefs = paddle.to_tensor(coefs)
H_A = [pauli_str_to_matrix([[1., pstr]], n=N//2) for pstr in pauli_str_A]
H_A = paddle.to_tensor(np.stack(H_A))
H_B = [pauli_str_to_matrix([[1., pstr]], n=N-N//2) for pstr in pauli_str_B]
H_B = paddle.to_tensor(np.stack(H_B))

# calculate the ground state information
ground_state_energy = get_ground_state_info(H_mtr)

The ground state energy is -0.99783 Ha.
Schmidt rank of the ground state is 32.

# Construct parameterized circuit
def U_cir(N, D):
    
    cir = Circuit(N)  # Initialize an N-qubit-width circuit
    cir.complex_entangled_layer('full', N, D)   # Add quantum gates
    return cir  

# Apply a parameterized circuit on the conputational bases
# and return a tensor of shape [2**N, num_states]
def output_states(num_states, N, cir):
    # Create num_states-many computational bases
    basis = paddle.eye(2**N, num_states)
    
    # Acquire a matrix of the parameterized circuit
    U = cir.unitary_matrix()
    
    # Apply the parameterized circuit on these bases
    vec = U @ basis                         
    
    return vec

# Construct the distributed model
class DistributedVQE(paddle.nn.Layer):
    def __init__(self, N, D, S):
        super().__init__()
        paddle.seed(SEED)

        self.S = S  # Define constant S
        self.N = N
        self.cir = [U_cir(N//2, D), U_cir(N - N//2, D)]
        
    # The core logic of distributed VQE
    def forward(self):
        # Obtain U|k> and V|k> for subsystems A and B respectively 
        vec_A = output_states(self.S, self.N//2, self.cir[0])
        vec_B = output_states(self.S, self.N - self.N//2, self.cir[1])
        
        # Calculate tensor E_A, E_B, which have elements E_{ijt}^A and E_{ijt}^B, as per defined in above
        E_A = vec_A.conj().t() @ H_A @ vec_A
        E_B = vec_B.conj().t() @ H_B @ vec_B
        M = (coefs.reshape([-1, 1, 1]) * E_A * E_B).sum(0)

        # Find the minimal eigenvalue of M
        eigval = paddle.linalg.eigvalsh(M)
        loss = eigval[0]
        
        return loss

def train(model):
    start_time = time.time()    # To calculate the running time of this function
    params = sum([cir.parameters() for cir in model.cir], [])
    # We will use Adam, a gradient-based optimizer to optimize theta and phi
    opt = paddle.optimizer.Adam(learning_rate=LR, parameters=params)
    summary_loss = []           # Save loss history

    # Optimization iteration
    for itr in range(ITR):
        # Forward propagation to calculates the loss function
        loss = model()
        # Backward propagation to optimize the loss function
        loss.backward()
        opt.minimize(loss)
        opt.clear_grad()
        # Update optimization result
        summary_loss.append(loss.numpy())
        # Print itermediary result
        if (itr+1) % 20 == 0:
            print(f"iter: {itr+1}, loss: {loss.tolist()[0]: .4f} Ha")

    print(f'Ground truth is  {ground_state_energy:.4f} Ha')
    print(f'Training took {time.time() - start_time:.2f}s')
    
    plt.plot(list(range(ITR)), summary_loss, color='r', label='loss')
    plt.hlines(y=ground_state_energy, xmin=0, xmax=ITR, linestyle=':',  label='ground truth')
    plt.legend()
    plt.title(f'Loss for {type(model).__name__} on a {N}-qubit Hamiltonian')
    plt.show()

# Note that we manually set S = 4 as the Hamiltonian we just created interacts weakly across the subsystems.
# (See the Conclusion section for further description)
vqe = DistributedVQE(N, D, S=4)
train(vqe)

iter: 20, loss: -0.9244 Ha
iter: 40, loss: -0.9906 Ha
iter: 60, loss: -0.9968 Ha
iter: 80, loss: -0.9977 Ha
iter: 100, loss: -0.9978 Ha
Ground truth is  -0.9978 Ha
Training took 7.24s

class StandardVQE(paddle.nn.Layer):
    def __init__(self, N, D, S):
        super().__init__()
        paddle.seed(SEED)
        self.cir = U_cir(N, D)
        self.loss_fcn = ExpecVal(H)
      
    def forward(self):
        output_state = self.cir()
        loss = self.loss_fcn(output_state)
        return loss.cast('float64').flatten()

svqe = StandardVQE(N, D, S=1)
train(svqe)  # Train the standard VQE

iter: 20, loss: -0.8365 Ha
iter: 40, loss: -0.9852 Ha
iter: 60, loss: -0.9958 Ha
iter: 80, loss: -0.9975 Ha
iter: 100, loss: -0.9978 Ha
Ground truth is  -0.9978 Ha
Training took 1278.13s

Distributed Variational Quantum Eigensolver Based on Schmidt Decomposition¶

Overview¶

Schmidt Decomposition¶

Distributed VQE Based on Schmidt Decomposition¶

Paddle Quantum implementation¶

Conclusion¶

References¶