<jupyter_start><jupyter_text>Estimation of Quantum State Properties Based on the Classical Shadow  Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved.   Overview In [The Classical Shadow of Unknown Quantum States](./ClassicalShadow_Intro_EN.ipynb), we introduced some theoretical knowledge of the classical shadow and showed how to construct the classical shadow of a quantum state $\rho$. According to the theoretical derivation of the classical shadow, it is very suitable for estimating the linear properties of quantum states. At present, its basic applications are: quantum fidelity estimation, entanglement verification, the estimation of local observables' expectation values, and the estimation of global observables' expectation values [1]. Among them, the estimation of observables' expectation values widely appears in current quantum algorithms, such as [Variational Quantum Eigensolver](./VQE_EN.ipynb)  (VQE) for estimating the ground state energy of complex molecular Hamiltonian. Next, we will focus on the estimation of observables' expectation values algorithms based on the classical shadow and show how to use the shadow function in Paddle Quantum to do so.  Estimation of Observables' Expectation Values  Problem description In the field of quantum chemistry, one of the core tasks is to solve for ground state energy of the Hamiltonian $\hat{H}$ on a closed physical system on a quantum scale and its corresponding ground state. The main method is to prepare a parameterized trial wave function $|\Psi(\theta)\rangle$. Then, the parameter $\theta$ is being continuously adjusted and optimized to minimize the expected value $\langle\Psi(\theta)|\hat{H}| \Psi(\theta)\rangle$ using classical optimization algorithms (such as gradient descent method). The principle of this scheme is based on Rayleigh-Ritz variational principle,$$E_{0}=\min _{\theta}\langle\Psi(\theta)|\hat{H}| \Psi(\theta)\rangle \tag{1}$$where $E_{0}$ represents the ground state energy of the system. Numerically, the problem can be understood as solving for the minimum eigenvalue $\lambda_{\min }$ of a discretized Hamiltonian $\hat{H}$ (Hermitian matrix) and its corresponding eigenvector $\left|\Psi_{0}\right\rangle$. Where the classical shadow comes into play is to calculate the $\langle\Psi(\theta)|\hat{H}| \Psi(\theta)\rangle = \operatorname{tr}(\hat{H}\rho )$ part in every optimization iteration where $ \rho = | \Psi(\theta)\rangle\langle\Psi(\theta)|$.The problem is then transformed into: for a quantum state $\rho$ of $n$ qubits and an observable (Hamiltonian) $\hat{H}$ that can be written as a linear combination of a set of Pauli operators $\{I, X, Y, Z\}^{\otimes n}$,$$\hat{H}=\sum_{Q \in\{I, X, Y, Z\} ^{\otimes n}} \alpha_{Q} Q \quad \text{where} \quad \alpha_{Q} \in \mathbb{R}, \tag{2}$$how to estimate the observable's expectation value $\operatorname{tr}(\hat{H}\rho )$ using the classical shadow?The most intuitive method is to use each term of the Hamiltonian as a measurement base, and the corresponding Pauli measurements are made for the quantum state $\rho$. A certain number of repetitions are made for each of the measurements, and then the measurement results are being processed to obtain the estimation value. Here we refer to this method as the item-by-item measurement method.Readers can see that when both $n$ and the number of terms of the Hamiltonian $\hat{H}$ are small, we can get $\operatorname{tr}(\hat{H}\rho )$ through the item-by-item measurement method. However, when $n$ is large and the number of terms of $\hat{H}$ increases, the cost of the this method will increase significantly. The classical-shadow-based method that will be introduced can obtain the same precision estimation of $\operatorname{tr}(\hat{H}\rho )$ with less cost.  The improved algorithm based on the classical shadow When constructing the classical shadow, a critical step is to uniformly and randomly sample the unitary transformation from the fixed set. In [The Classical Shadow of Unknown Quantum States](./ClassicalShadow_Intro_EN.ipynb), we showed the case when the selected set is the Clifford group. When the selected set is a Clifford group on single qubit, the process of sampling and measurement is equivalent to making Pauli measurements on the quantum states. The classical shadow algorithm using random Pauli measurements (CS) is provided in Paddle Quantum. Briefly, in the CS algorithm, we repeatedly choose a Pauli basis for each qubit uniformly at random to measure the quantum state $\rho$ and estimate the observable's expectation value based on the measurement results, and the reader may refer to [1-2] to learn the details of the principle. Further, when the Pauli measurement base is no longer selected uniformly and randomly, improved algorithms are proposed [2-3]. Relevant algorithm functions are also provided in Paddle Quantum: Locally-biased classical shadows (LBCS) [2], Adaptive Pauli shadows (APS) [3]. Readers can refer to [1-3] to learn these algorithms in detail.  Paddle Quantum Implementation In Paddle Quantum, we provide the shadow function, which mainly includes two parts to use the above three algorithms to estimate the expectation value of the observable and obtain the classical shadow data of unknown quantum states. Next, we will show how to implement finding the ground state energy estimation of $H_{2}$ and $LiH$ based on the shadow function in Paddle Quantum.<jupyter_code># Import all the dependencies
import numpy as np
from numpy import pi as PI
import paddle
from paddle_quantum.circuit import UAnsatz
from paddle_quantum.VQE.chemistrysub import H2_generator
from paddle_quantum.utils import Hamiltonian<jupyter_output><empty_output><jupyter_text>Estimate the ground state energy of hydrogen molecule  ($H_{ 2}$) Import the Hamiltonian of $H_2$ with 4 qubits (for details, please refer to [Variational Quantum Eigensolver](./VQE_CN.ipynb) tutorial to obtain $H_2$ molecular Hamiltonian).<jupyter_code># Set up the Hamiltonian of hydrogen molecule
H2_pauli_str, H2_qubit = H2_generator()
# Construct a Hamiltonian class instance using the H2_pauli_str
H2_hamiltonian = Hamiltonian(H2_pauli_str)
print('H2 hamiltonian = ', H2_hamiltonian)<jupyter_output>H2 hamiltonian =  -0.04207897647782277 I0
0.17771287465139946 Z0
0.1777128746513994 Z1
-0.2427428051314046 Z2
-0.24274280513140462 Z3
0.17059738328801055 Z0, Z1
0.04475014401535163 Y0, X1, X2, Y3
-0.04475014401535163 Y0, Y1, X2, X3
-0.04475014401535163 X0, X1, Y2, Y3
0.04475014401535163 X0, Y1, Y2, X3
0.12293305056183797 Z0, Z2
0.1676831945771896 Z0, Z3
0.1676831945771896 Z1, Z2
0.12293305056183797 Z1, Z3
0.1762764080431959 Z2, Z3<jupyter_text>To show how to estimate the ground state energy using the classical-shadow-based algorithm, we first get the estimated ground state of $H_{2}$ using the VQE algorithm in Paddle Quantum.<jupyter_code>def U_theta(theta, hamiltonian, N, D):
    """
    Quantum Neural Network
    """
    
    # Initialize the quantum neural network according to the number of qubits N
    cir = UAnsatz(N)
    
    # Built-in {R_y + CNOT} circuit template
    cir.real_entangled_layer(theta[:D], D)
    
    # Lay R_y gates in the last row
    for i in range(N):
        cir.ry(theta=theta[D][i][0], which_qubit=i)
        
    # The quantum neural network acts on the default initial state |0000>
    cir.run_state_vector()
    
    # Calculate the expected value of a given Hamiltonian
    expectation_val = cir.expecval(hamiltonian)

    return expectation_val, cir

class StateNet(paddle.nn.Layer):
    """
    Construct the model net
    """

    def __init__(self, shape, dtype="float64"):
        super(StateNet, self).__init__()
        
        # Assign the theta parameter list to be the trainable parameter list of the circuit
        self.theta = self.create_parameter(shape=shape, 
                                           default_initializer=paddle.nn.initializer.Uniform(low=0.0, high=2*PI),
                                           dtype=dtype, is_bias=False)
        
    # Define loss function and forward propagation mechanism
    def forward(self, hamiltonian, N, D):
        
        # Calculate the loss function/expected value
        loss, cir = U_theta(self.theta, hamiltonian, N, D)

        return loss, cir<jupyter_output><empty_output><jupyter_text>After building the quantum neural network of the VQE algorithm and defining the loss function, we can estimate the ground state energy of $H_{2}$ and the quantum circuit corresponding to the ground state by training the quantum neural network.<jupyter_code>ITR = 80  # Set the number of optimization iterations
LR = 0.4   # Set the learning rate
D = 2      # Set the depth of the repetitive calculation module in QNN
N = H2_hamiltonian.n_qubits 

# Determine the parameter dimension of the network
net = StateNet(shape=[D + 1, N, 1])

# Generally speaking, we use Adam optimizer to obtain relatively good convergence,
# You can change it to SGD or RMS prop.
opt = paddle.optimizer.Adam(learning_rate=LR, parameters=net.parameters())

# Record optimization results
summary_iter, summary_loss = [], []

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

    # Forward propagation to calculate loss function
    loss, cir = net(H2_hamiltonian, N, D)

    # Use back propagation to minimize the loss function
    loss.backward()
    opt.minimize(loss)
    opt.clear_grad()

    # Record optimization results
    summary_loss.append(loss.numpy())
    summary_iter.append(itr)

    # Print result
    if itr % 20 == 0:
        print("iter:", itr, "loss:", "%.4f" % loss.numpy())
        print("iter:", itr, "Ground state energy:", "%.4f Ha" 
                                            % loss.numpy())
    if itr == ITR:
        print("\nThe trained circuit:") 
        print(cir)<jupyter_output>iter: 20 loss: -1.0865
iter: 20 Ground state energy: -1.0865 Ha
iter: 40 loss: -1.1284
iter: 40 Ground state energy: -1.1284 Ha
iter: 60 loss: -1.1355
iter: 60 Ground state energy: -1.1355 Ha
iter: 80 loss: -1.1360
iter: 80 Ground state energy: -1.1360 Ha

The trained circuit:
--Ry(-1.57)----*--------------x----Ry(4.713)----*--------------x----Ry(3.129)--
               |              |                 |              |               
--Ry(5.091)----x----*---------|----Ry(1.768)----x----*---------|----Ry(3.711)--
                    |         |                      |         |               
--Ry(0.997)---------x----*----|----Ry(4.696)---------x----*----|----Ry(3.180)--
                         |    |                           |    |               
--Ry(4.753)--------------x----*----Ry(4.701)--------------x----*----Ry(-0.37)--<jupyter_text>Introduction to shadow function At this point, we've obtained the quantum circuit for generating the ground state of $H_{2}$. We can run the `shadow_trace` function on this circuit to obtain the estimated ground state energy. In the `shadow_trace` function, our inputs are the Hamiltonian to be estimated, the number of samples, and the method selected. You can choose the sampling mode you want by specifying the parameter `method`.  Among them, CS has a broader application range and the fastest speed, but its estimation accuracy may be slightly poor; LBCS has a higher accuracy, but it runs slower as the number of terms of Hamiltonian grows; APS also has higher accuracy, but it runs slower when the number of qubits is large.<jupyter_code># The actual energy value corresponding to the estimated ground state
H2_energy = cir.expecval(H2_hamiltonian).numpy()[0]

# Sampling times
sample = 1500 
# Three algorithms are used to estimate the expectation value of observable.
H2_energy_CS = cir.shadow_trace(H2_hamiltonian, sample, method="CS")
H2_energy_LBCS = cir.shadow_trace(H2_hamiltonian, sample, method="LBCS")
H2_energy_APS = cir.shadow_trace(H2_hamiltonian, sample, method="APS")

print('H2 ground state energy = ', H2_energy)
print('H2 ground state energy CS= ', H2_energy_CS)
print('H2 ground state energy LBCS= ', H2_energy_LBCS)
print('H2 ground state energy APS= ', H2_energy_APS)<jupyter_output>H2 ground state energy =  -1.1360331063822373
H2 ground state energy CS=  -1.1652093330967717
H2 ground state energy LBCS=  -1.1213982275348622
H2 ground state energy APS=  -1.137720214050516<jupyter_text>Now let's use the item-by-item measurement method to estimate the ground state energy of $H_{2}$.<jupyter_code># Use the item-by-item measurement method to estimate the ground state energy
H2_energy_traditional = cir.expecval(H2_hamiltonian, shots=100).numpy()[0]
print('H2 ground state energy traditional = ', H2_energy_traditional)<jupyter_output>H2 ground state energy traditional =  -1.1367622727687419<jupyter_text>We can see that under 1500 samples, the estimated ground state energy by these three algorithms are very close to the energy of the estimated ground state by VQE. The item-by-item measurement method is to make 100 measurements for each item of the Hamiltonian, and there are 15 items of the Hamiltonian of $H_{2}$, which is equivalent to 1500 measurements in total. The difference between the obtained results and the energy of the estimated ground state by VQE is also tiny. In this small-scale situation, the classical-shadow-based algorithms do not show significant advantages over the item-by-item measurement method. But in large-scale qubits scenarios, this algorithm requires only constant-level growth of the number of the Hamiltonian terms. In contrast, the item-by-item measurement method or some other methods require polynomial-level or even exponential-level growth in the number of samples to get the same accuracy[1]. In fact, it is pointed out in [2] that for CS algorithm and LBCS algorithm, the average error of estimation $\epsilon$, variance $\operatorname {var}(\nu)$ and number of samples $S$ are related as follows:$$S = O(\epsilon^{-2} \operatorname{var}(\nu) ),  \tag{3}$$where the variance $\operatorname{var} (\nu) $ is independent of the number of samples, and the variance is related to the number of terms of the Hamiltonian. Therefore, knowing our expected precision, we can calculate the required number of samples. Similarly, our average error can also be defined according to the number of samples: $$\epsilon = \sqrt{\frac{\operatorname{var}}{S}}.  \tag{4}$$We can be seen that in the above experiments, the errors obtained by the CS algorithm and LBCS algorithm are within the accuracy of theoretical estimation in reference [3]. At the same time, Paddle Quantum provides a sampling function ` shadow_sample`, which supports pre-sampling of unknown quantum states, which is convenient for readers to explore other applications of the classical shadow. The specific usage is as follows:<jupyter_code>from paddle_quantum.shadow import shadow_sample

# Run the circuit in the vector form
H2_rho = np.array(cir.run_state_vector())
# Get the data of classical shadow and output it in the form of a list
H2_sample_data_CS = shadow_sample(H2_rho, H2_qubit, sample_shots=10, mode='state_vector', 
                                  hamiltonian=H2_hamiltonian, method='CS')
print(H2_sample_data_CS)<jupyter_output>[('yxxz', '0010'), ('xyzx', '0101'), ('yxyx', '1010'), ('yxyy', '1101'), ('yyyx', '0101'), ('xyxz', '0010'), ('zyxz', '1010'), ('xzyz', '0110'), ('zzzx', '1101'), ('xyxx', '0101')]<jupyter_text>Estimate the ground state energy of lithium hydride  ($LiH$) Next, we consider the ground state energy of $LiH$. First, we load from a pre-computed file to generate the molecular Pauli Hamiltonian of $LiH$ with 12 qubits.<jupyter_code>with open('./LiH_hamiltonian.txt', 'r') as lih_file:
    unprocessed_pauli_str = lih_file.read()
    LiH_pauli_str = [term.split(maxsplit=1) for term in unprocessed_pauli_str.split('\n')]
    LiH_pauli_str = [[float(term[0]), term[1]] for term in LiH_pauli_str]
    LiH_hamiltonian = Hamiltonian(LiH_pauli_str)<jupyter_output><empty_output><jupyter_text>Then, we also use the VQE algorithm to obtain the ground state energy and the ground state of $LiH $ and then use the classical shadow algorithm to estimate the ground state energy. Due to the large size of the $LiH$ molecular Hamiltonian, it will take a long time to train the VQE circuit. Hence we provide a set of pre-trained parameters, using which the users could test the classical-shadow-based algorithms on $LiH$ Hamiltonian directly.<jupyter_code># Load the pre-trained parameters
pretrained_parameters = paddle.load('LiH_VQE_parameters.pdtensor')
N = LiH_hamiltonian.n_qubits
# Run the VQE circuit with pre-trained parameters
energy, cir = U_theta(pretrained_parameters, LiH_hamiltonian, N, D)
print('The pre-trained VQE gets a ground state energy of: %.4f ' % energy.numpy())<jupyter_output>The pre-trained VQE gets a ground state energy of: -7.7720<jupyter_text>Once we have the circuit corresponding to the $LiH$ ground state, we can directly use the `shadow_trace` function to  perform random measurements. Also, since this molecular Hamiltonian has 631 terms, we specify `sample = 1262` for the function `shadow_trace` and `shots = 2` for the function `expecval` in order to ensure that the number of measurements is the same for both types of methods.Since ground state of $LiH$ has 12 qubits, there are $3^{12}$ possible combinations of measurements when doing different Pauli measurements on the ground state of $LiH$. So it is too random to just perform 1262 samples to get the valuation. Thus, we run each of the above four methods 20 times. Then take the mean of these 20 samples of data as the estimation for each algorithm, and calculate the sample variance for a simple comparison of the algorithms.（It may take at least 1 hour to run the following code blocks.）<jupyter_code>import time

begin = time.time()
estimator_CS = []
estimator_LBCS = []
estimator_APS = []
estimator_traditional = []

# The actual energy value corresponding to the estimated ground state
LiH_energy = cir.expecval(LiH_hamiltonian).numpy()[0]

# Number of repetition times
n = 20 

for i in range(n):
    LiH_energy_CS = cir.shadow_trace(LiH_hamiltonian, 1262, method="CS")
    LiH_energy_LBCS = cir.shadow_trace(LiH_hamiltonian, 1262, method="LBCS")
    LiH_energy_APS = cir.shadow_trace(LiH_hamiltonian, 1262, method="APS")
    LiH_energy_traditional = cir.expecval(LiH_hamiltonian, shots=2).numpy()[0]

    estimator_CS.append(LiH_energy_CS) 
    estimator_LBCS.append(LiH_energy_LBCS) 
    estimator_APS.append(LiH_energy_APS) 
    estimator_traditional.append(LiH_energy_traditional) 

ave_LiH_energy_CS = np.mean(estimator_CS)
ave_LiH_energy_LBCS = np.mean(estimator_LBCS)
ave_LiH_energy_APS = np.mean(estimator_APS)
ave_LiH_energy_traditional = np.mean(estimator_traditional)
end = time.time() 

print("LiH ground state energy = ", LiH_energy)
print("ave LiH ground state energy CS = ", ave_LiH_energy_CS)
print("ave LiH ground state energy LBCS = ", ave_LiH_energy_LBCS)
print("ave LiH ground state energy APS = ", ave_LiH_energy_APS)
print('ave LiH ground state energy traditional = ', ave_LiH_energy_traditional)
print('time = ', end-begin)<jupyter_output>LiH ground state energy =  -7.771980394176657
ave LiH ground state energy CS =  -7.835791570579005
ave LiH ground state energy LBCS =  -7.7622296662623445
ave LiH ground state energy APS =  -7.762836542787509
ave LiH ground state energy traditional =  -7.8964746269601465
time =  4206.216086864471<jupyter_text>From the results, the mean values obtained by the classical-shadow-based algorithms are closer to the energy of the estimated ground state of $LiH$ by VQE than the item-by-item measurement method, and the errors of the algorithms are all within the accuracy of theoretical estimation in reference [3]. So, what about the sample variance of each algorithm?<jupyter_code># Calculate the sample variance
variance_CS = []
variance_LBCS = []
variance_APS = []
variance_traditional = []

for i in range(n):
    variance_CS.append((estimator_CS[i] - ave_LiH_energy_CS) ** 2)
    variance_LBCS.append((estimator_LBCS[i] - ave_LiH_energy_LBCS) ** 2)
    variance_APS.append((estimator_APS[i] - ave_LiH_energy_APS) ** 2)
    variance_traditional.append((estimator_traditional[i] - ave_LiH_energy_traditional) ** 2)

var_CS = sum(variance_CS)/(n-1)
var_LBCS = sum(variance_LBCS)/(n-1)
var_APS = sum(variance_APS)/(n-1)
var_traditional = sum(variance_traditional)/(n-1)

print('LiH variance CS = ', var_CS)
print('LiH variance LBCS = ', var_LBCS)
print('LiH variance APS = ', var_APS)
print('LiH variance traditional = ', var_traditional)<jupyter_output>LiH variance CS =  0.034596840359755784
LiH variance LBCS =  0.016602696085670984
LiH variance APS =  0.0016603026356630662
LiH variance traditional =  0.13200055652163223
{'fraction_non_alphanumeric': 0.0907207554442089, 'fraction_numerical': 0.04249373675081904, 'mean_word_length': 4.241666666666666, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '367346', 'license': 'Apache-2.0', 'n_tokens_mistral': 6544, 'n_tokens_neox': 5861, 'n_words': 2414, 'path': 'tutorial/quantum_simulation/ClassicalShadow_Application_EN.ipynb', 'repo_name': 'imppresser/Quantum-1', 'size': '20756'}