# Import related packages
import paddle
from paddle_quantum.gate import Oracle
from paddle_quantum.ansatz import Circuit
from paddle_quantum.state import State, zero_state
from paddle_quantum.loss import StateFidelity
import numpy as np
from numpy import pi, log
from scipy import optimize
import warnings
warnings.filterwarnings('ignore')

# Define the unitary operator A
n_qubits = 1 # Single qubit
amp_operator = Circuit(n_qubits) # Construct the single-qubit unitary operator A
set_param = pi/3 # The parameter for the unitary operator A
amp_param = paddle.to_tensor([set_param, 0, 0], dtype='float32') # Transfer the parameter to tensor
amp_operator.u3(param=amp_param) # Input the parameter to a U3 gate to construct A

initial_state = zero_state(n_qubits) # Define the initial state as zero state
chi_0_state =amp_operator(initial_state) # The effect of applying A on zero state
quantum_amp =chi_0_state.measure()['1'] # Output the pre-set quantum amplitude for estimation
print(f"The quantum amplitude correspond to the unitary operator A is {quantum_amp}")

The quantum amplitude correspond to the unitary operator A is 0.25

# Construct the reflection operator R_chi
identity = paddle.to_tensor([[1, 0],[0, 1]], dtype="complex64") # Define the identity matrix
density_matrix_0 = chi_0_state.ket @ chi_0_state.bra # amp_0_state in the form of density matrix
r_chi = identity - paddle.to_tensor([2], dtype="complex64") * density_matrix_0 # Construct the reflection operator respect to state chi

# Construct the reflection operator R_good
flip = Circuit(n_qubits)
flip.x() 
one_state = flip(initial_state) # Construct the "one state" which is defined as good state
density_matrix_good = one_state.ket @ one_state.bra # Good state in the form of density matrix
r_good = identity - paddle.to_tensor([2], dtype="complex64") * density_matrix_good # Construct the reflection operator respect to the good state

# Construct the amplification operator Q
num_amplification = 1 # The total times of apply Q
Q_operator = Oracle(paddle.to_tensor([-1], dtype="complex64") * r_chi @ r_good, qubits_idx=0, num_qubits=1, depth=num_amplification) # Use Oracle to construct the amplification operator

# Classical Monte-Carlo method
N_sample_mc = 20000 # The total amount of sampling
N_good = 0 # Initialize the times of measuring good state
for i in range(N_sample_mc):
    result = chi_0_state.measure(shots=1) # Sample the chi state once
    N_good += result['1'] # If the measurement result is "1", N_good plus one
amp_mc = N_good/N_sample_mc # Estimate the quantum amplitude by the probability of mearsuring "1"
error_mc = abs(amp_mc - quantum_amp) # Absolute error of the estimation 
rate_mc = error_mc/quantum_amp # Relative error of the estimation
print(f"CMC Estimation result is a = {amp_mc}, absolute error is {error_mc}, relative error is {100*rate_mc}%")
print(f"CMC Query complexity is {N_sample_mc}")

CMC Estimation result is a = 0.24915, absolute error is 0.0008499999999999897, relative error is 0.33999999999999586%
CMC Query complexity is 20000

# Step 1 of MLAE

# Initialize parameters 
M = 25 # Maximum times of implementing Q
m = np.arange(0, M, 1) #LIS
N_k = 25 # Total amount of sampling
h_k = np.zeros(M) # Initialize the data space to store the times of measuring good state.

# Sampling process
for k in range(M):
    Q_operator_MLAE = Oracle(paddle.to_tensor([-1], dtype="complex64") * r_chi @ r_good, qubits_idx=0, num_qubits=1, depth=k)
    for i in range(N_k):
        rotated_state = Q_operator_MLAE(chi_0_state) # Implement amplification operator on initial state
        result = rotated_state.measure(shots = 1) # Measure and record the number of measuring good state
        h_k[k] += result['1']  # Store the number of good state for different "k"

# Step 2 of MLAE 
# Definition of the likelihood function (in logarithmic form)
params = ()
def likelihood(z, *params):
    theta = z
    f = 0
    for i in range(M):
        f = f + log(np.sin((2 * m[i] + 1) * theta) ** (2 * h_k[i]) * np.cos((2 * m[i] + 1) * theta) ** (2 * (N_k - h_k[i])))
    return (-f) # the following algorithm will find the minimum, minus sign corresponds the result to maximum.

# Step 3 of MLAE
# Use Brute-force search algorithm to find the maximum value of likelihood function
rranges = (0, pi/2, pi/200) # Range of grid for Brute-force algorithm
resbrute = optimize.brute(likelihood, (rranges,), args=params, full_output=True, finish=optimize.fmin) # Call Brute-force algorithm
theta_a = resbrute[0][0] # Theta corresponds to the minimum of -f
amp_MLAE = np.sin(theta_a) ** 2 # Quantum amplitude estimation from MLAE
error_MLAE = abs(amp_MLAE - quantum_amp) # Absolute estimation error
rate_MLAE = error_MLAE/quantum_amp # Relative estimation error
print(f"The maximum value for likelihood function is {-resbrute[1]} when theta = {resbrute[0][0]} rad")
print(f"MLAE result is {amp_MLAE}, absolute error is {error_MLAE}, relative error is {100 * rate_MLAE}%")

The maximum value for likelihood function is -243.068172605532 when theta = 0.5233743030763384 rad
MLAE result is 0.249805626294017, absolute error is 0.00019437370598299197, relative error is 0.07774948239319679%

# Quantum query complexity for MLAE 
N_q_MLAE = 0
for i in range(M):
    N_q_MLAE += N_k * (2 * i + 1)
print(f"The quantum query complexity for MLAE is {N_q_MLAE}")

The quantum query complexity for MLAE is 15625

# Define related functions for VQAE

def construct_Q(input_param: float, n_qubits: int, num_amplification):
    r""" Construct amplification operator Q.
    
    Args:
        input_param: input parameters for U3 rotation gate
        n_qubits: number of qubits in circuit
        num_amplification: number of implementing amplification operator
    
    Returns:
        Oracle: amplification operator with input parameter and number of implementation
    
    """
    # Define unitary operator A
    amp_operator = Circuit(n_qubits) 
    amp_param = paddle.to_tensor([input_param, 0, 0], dtype='float32') # Preset parameter for unitary operator
    amp_operator.u3(param=amp_param)

    initial_state = zero_state(n_qubits) # Define the initial state as zero state
    chi_0_state =amp_operator(initial_state) # Apply unitary operator A on zero state 

    # Construct reflection operator R_chi
    identity = paddle.to_tensor([[1, 0],[0, 1]], dtype="complex64") # Define identity operator
    density_matrix_0 = chi_0_state.ket @ chi_0_state.bra # Density matrix for amp_0_state
    r_chi = identity - paddle.to_tensor([2], dtype="complex64") * density_matrix_0

    # Construct reflection operator R_good
    flip = Circuit(n_qubits)
    flip.x()
    one_state = flip(initial_state) # Construct "one state" which is defined as good state
    density_matrix_good = one_state.ket @ one_state.bra # Density matrix for good state
    r_good = identity - paddle.to_tensor([2], dtype="complex64") * density_matrix_good
    
    # Return the amplification operator Q
    return Oracle(paddle.to_tensor(paddle.to_tensor([-1], dtype="complex64") * r_chi @ r_good), qubits_idx=0, num_qubits=n_qubits, depth=num_amplification)


def loss_fcn(input_state: State, period: int, start_state: State, target_param: float) -> paddle.Tensor:
    r""" Calculate the loss function for neural network.
    
    Args:
        input_state: input state
        period: period of variational process 
        start_state: current variational quantum state, for constructing target state
        target_param: parameter of target state
    
    Returns:
        loss: the value of loss function which is defined as the infidelity between the input state and target state
    
    """

    Q_var = construct_Q(target_param, int(1), period)
    Target_state = Q_var(start_state) 
    loss = 1 - StateFidelity(Target_state)(input_state)
    return loss


def cir_constructor(depth: int) -> Circuit:
    r""" Construct variational quantum circuit
    
    Args:
        depth: the depth of quantum circuit
    
    Returns:
        variational quantum circuit
    
    Note:
        For single qubit, one layer is defined as a combination of Ry and Rz rotation gate.
        The rotation angles are the parameters of quantum circuit.
    
    """
    cir = Circuit(1)
    
    for _ in range(depth):
        cir.ry()
        cir.rz()
    return cir

# Initialization of parameters (Notice: VQAE is based on the framework of MLAE)
M = 25 # Maximum times of implementing Q
m = np.arange(0, M, 1) # Linear increment sequence
N_k_vqae = 50 # Total amount of sampling

# Superparameters for variational approximation
theta_size = 6 # Amount of variational parameters (depends on the depth of circuit)
ITR = 500  # Number of iterations
LR = 0.001  # Learning rate 
SEED = 1  # Random number seed
paddle.seed(SEED)

<paddle.fluid.core_avx.Generator at 0x1fb5ba83d30>

# Sampling and variational approximation process
start_state = chi_0_state #start_state is the variational quantum state, equals to chi_0_state in the beginning
cycle = 5 # Variational period
h_k_vqae = np.zeros(M) # Initialize the data space to store the times of measuring good state.
for k in range(M):
    i = np.floor(k / cycle)
    j = k % cycle
    Q_operator_VQAE = construct_Q(set_param, int(1), j)
    for sample in range(N_k_vqae):
        rotated_state = Q_operator_VQAE(start_state) # Implement amplification operator on initial state
        result = rotated_state.measure(shots = 1) # Measure and record the number of measuring good state
        h_k_vqae[k] += result['1']  # Store the number of good state for different "k"

    if j == cycle - 1: # Condition to trigger variational approximation process
        # Create data space to store loss function
        loss_list = []

        # Define the variational circuit
        cir = cir_constructor(3)

        # Use Adam optimizer
        opt = paddle.optimizer.Adam(learning_rate=LR, parameters=cir.parameters())
        
        # Optimization iterations
        for itr in range(ITR):
            # Forward optimize the loss function
            loss = loss_fcn(cir(chi_0_state), cycle, start_state, set_param)

            # Backward optimize the loss function
            loss.backward()
            opt.minimize(loss)
            opt.clear_grad()

            # Record the loss function (learning curve)
            loss_list.append(loss.item())
        
        print("=================================")
        print(f'The {int(i + 1)} th variational approximation')
        print(f'The minimum of loss function  {loss_list[-1]}')
        print(f"Optimal variational parameters {cir.param.numpy()}")
        
        Q_test = construct_Q(set_param, int(1), cycle)
        test_state = Q_test(start_state) # Update the test state for fidelity calculation
        
        start_state = cir(chi_0_state)  # Update the variational quantum state 
        start_state.data.stop_gradient = True 
        print(f"The state fidelity between variational quantum state and target state is {StateFidelity(test_state)(start_state).numpy()[0]}")

=================================
The 1 th variational approximation
The minimum of loss function  0.011207759380340576
Optimal variational parameters [-0.39857474  4.4845576   3.8463612   0.5077383   2.1137552   5.696977  ]
The state fidelity between variational quantum state and target state is 0.9888777136802673
=================================
The 2 th variational approximation
The minimum of loss function  0.02215433120727539
Optimal variational parameters [3.538579  3.6868277 3.778493  0.8545991 3.3457727 5.151529 ]
The state fidelity between variational quantum state and target state is 0.9780022501945496
=================================
The 3 th variational approximation
The minimum of loss function  0.12489765882492065
Optimal variational parameters [2.9717562 4.5942483 3.8417864 1.0915955 2.9404604 4.45632  ]
The state fidelity between variational quantum state and target state is 0.8756692409515381
=================================
The 4 th variational approximation
The minimum of loss function  0.0724669098854065
Optimal variational parameters [5.293609  5.0283127 1.7944657 3.575252  2.10594   4.743199 ]
The state fidelity between variational quantum state and target state is 0.9278923273086548
=================================
The 5 th variational approximation
The minimum of loss function  0.007037222385406494
Optimal variational parameters [ 2.9615376   6.2987323   2.9560285  -0.05320466  1.5106332   3.6751769 ]
The state fidelity between variational quantum state and target state is 0.9930309057235718

# Definition of the likelihood function (in logarithmic form)
params = ()
def likelihood_vqae(z, *params):
    theta = z
    f = 0
    for i in range(M):
        f = f + log(np.sin((2 * m[i] + 1) * theta) ** (2 * h_k_vqae[i]) * np.cos((2 * m[i] + 1) * theta) ** (2 * (N_k_vqae - h_k_vqae[i])))
    return (-f)# the following algorithm will find the minimum, minus sign corresponds the result to maximum.

# Step 3 of MLAE
# Use Brute-force search algorithm to find the maximum value of likelihood function
rranges = (0, pi/2, pi/200) # Range of grid for Brute-force algorithm
resbrute = optimize.brute(likelihood_vqae, (rranges,), args=params, full_output=True, finish=optimize.fmin)  # Call Brute-force algorithm
theta_a = resbrute[0][0] # Theta corresponds to the minimum of -f
amp_VQAE = np.sin(theta_a) ** 2 # Quantum amplitude estimation from VQAE
error_VQAE = abs(amp_VQAE - quantum_amp) # Absolute estimation error
rate_VQAE = error_VQAE/quantum_amp # Relative estimation error
print("=================================")
print(f"The maximum value for likelihood function is {-resbrute[1]} when theta = {resbrute[0][0]} rad")
print(f"VQAE result is {amp_VQAE}, absolute error is {error_VQAE}, relative error is {100 * rate_VQAE}%")

=================================
The maximum value for likelihood function is -728.8282960799222 when theta = 0.5310066244864633 rad
VQAE result is 0.2564425882363815, absolute error is 0.006442588236381497, relative error is 2.5770352945525987%

# Calculate the quantum query complexity of VQAE
from math import floor
M_vqae = M
k_cycle = cycle
n_p = theta_size
n_s = ITR
n_f = n_s
N_var1 = 2 * n_p * n_s * n_f
N_var = N_var1 * (2 * k_cycle + 2) * floor(M_vqae/k_cycle)
N_samp = N_k_vqae * k_cycle * (k_cycle + 2) * floor(M_vqae/k_cycle) + N_k_vqae * (M_vqae % k_cycle) * (M_vqae % k_cycle + 2)
print(f"Sampling complexity is {N_samp}")
print(f"Variational complexity is {N_var}")
print(f"Quantum query complexity for VQAE is {N_var + N_samp}")

Sampling complexity is 8750
Variational complexity is 180000000
Quantum query complexity for VQAE is 180008750

Variational quantum amplitude estimation¶

Backgrounds¶

Brief introduction to QAE¶

Classical Monte-Carlo algorithm (CMC)¶

Maximum likelihood amplitude estimation（MLAE）¶

MLAE by Paddle Quantum¶

Comparison between MLAE and CMC¶

Variational quantum amplitude estimation¶

Pseudo code for VQAE¶

Variational approximation process¶

VQAE implementation by Paddle Quantum¶

Summary¶

References¶