Variational quantum linear solver论文复现代码

复现工作思想：

 成果：采样1000，迭代10次，hardware-efficient ansatz,3qubit,3层ansatz,利用代码随机生成5种类型的ansatz(单旋转比特门顺序、数目不同)，每一种ansatz调参100次，代价函数成功收敛到0。

全局代价函数：

 全局代价函数计算：

电路搭建：（原文的hard overlap test电路不理解如何得到论文中的数值，我的实现利用的是swap test电路）

 

我后期进行思考的时候对PPT内容进行了修改。即使内积是含有虚部的，我们 利用swap test是计算不出虚部的，但是这里有计算实部和虚部，主要是输出态形式和swap test输出态形式不同：

 

复现代码：

（1）实现fixed ansatz，hardmard test && swap test 电路的搭建以及验证电路搭建是否正确

import qiskit
import math
from math import pi
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister,IBMQ,execute
from qiskit import Aer, transpile, assemble
from qiskit.providers.aer import QasmSimulator
from qiskit.visualization import circuit_drawer, plot_histogram
import random
import numpy as np
from scipy.optimize import minimize
from math import sqrt
# 一个单元格有多个输出
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = 'all'



 

# 定义函数fixed-ansatz,用于搭建门电路V，V|000>=|x(alpha)>
def apply_fixed_ansatz(circ, qubits, parameters,count_type):
    # count_type列表指定了每一层单旋转比特门的组成，列表parameters指定每个单比特旋转门的角度
     if((len(qubits)==3) or (len(qubits)==6)):
        for iz in range(0, 3):
            #随机数为1:旋转门x,随机数为2：旋转门y,随机数为3：旋转门z
            count_gate= count_type[0][iz]
            if count_gate ==1:
                circ.rx(parameters[0][iz], qubits[iz])
            if count_gate==2:
                circ.ry(parameters[0][iz], qubits[iz])
            if count_gate==3:
                circ.rz(parameters[0][iz], qubits[iz])
                
        
        circ.cz(qubits[0], qubits[1])
        circ.cz(qubits[2], qubits[0])

        for iz in range(0, 3):
            #生成3个随机数,随机数为1:旋转门x,随机数为2：旋转门y,随机数为3：旋转门z
            count_gate= count_type[1][iz]
            if count_gate ==1:
                circ.rx(parameters[1][iz], qubits[iz])
            if count_gate==2:
                circ.ry(parameters[1][iz], qubits[iz])
            if count_gate==3:
                circ.rz(parameters[1][iz], qubits[iz])        
        

        circ.cz(qubits[1], qubits[2])
        circ.cz(qubits[2], qubits[0])

        for iz in range(0, 3):
            #随机数为1:旋转门x,随机数为2：旋转门y,随机数为3：旋转门z
            count_gate= count_type[2][iz]
            if count_gate ==1:
                circ.rx(parameters[2][iz], qubits[iz])
            if count_gate==2:
                circ.ry(parameters[2][iz], qubits[iz])
            if count_gate==3:
                circ.rz(parameters[2][iz], qubits[iz])         
        
    
            
   

            
# Creates the Hadamard test,for calculating <psi|psi>
def had_test(gate_type, qubits, auxiliary_index, parameters,count_type):
    # 创建局部变量，使得每执行一次循环，就会生成一个新的电路
    circ = QuantumCircuit(4,4)
    circ.h(auxiliary_index)
    apply_fixed_ansatz(circ, qubits, parameters,count_type)

    for ie in range(0, len(gate_type[0])):
        if (gate_type[0][ie] == 1):
            circ.cz(auxiliary_index, qubits[ie])
#         if (gate_type[0][ie] == 0):
#             circ.cx(auxiliary_index, qubits[ie])
#             circ.cx(auxiliary_index, qubits[ie])


    for ie in range(0, len(gate_type[1])):
        if (gate_type[1][ie] == 1):
            circ.cz(auxiliary_index, qubits[ie])
#         if (gate_type[1][ie] == 0):
#             circ.cx(auxiliary_index, qubits[ie])
#             circ.cx(auxiliary_index, qubits[ie])
            
    
    circ.h(auxiliary_index)
    
    #测量辅助比特
    circ.measure(0,0)
    return circ

# #验证函数had_test()是否正确
# circ1 = had_test([[0,0,0], [0,0,0]],[1, 2, 3], 0, [[1,1,1],[1,1,1],[1,1,1]],[[1,3,3],[1,2,1],[2,3,1]])
# #电路可视化
# circ1.draw(output='mpl')
# simulator = Aer.get_backend('qasm_simulator')
# job = execute(circ1,simulator,shots=1000)
# result = job.result()
# # Returns counts
# counts = result.get_counts(circ1)
# # 统计辅助比特测得1的概率
# overall_sum_1=0
# if ('0001' in counts.keys()):
#     m_sum = float(counts["0001"])/1000
# else:
#     m_sum = 0
                
# overall_sum_1+=(1-2*m_sum)
# print('the final outcome of had_test: '+str(overall_sum_1))
# # 输出直方图
# plot_histogram(counts)



# Creates the swap test, for calculating |<b|psi>|^2

def swap_test(gate_type, qubits, auxiliary_index, parameters,count_type):
    circ = QuantumCircuit(7,7)
#     #制备|b>
#     circ.x(6)
#     circ.h(5)
    # swap test左侧H门
    circ.h(auxiliary_index)
    # 实现V门，用于制备|x(alpha)>
    apply_fixed_ansatz(circ,qubits, parameters,count_type)
    

    # 实现qubit1,2,3上的受控A
    for ie in range(0, len(gate_type[0])):
        if (gate_type[0][ie] == 1):
            circ.cz(auxiliary_index, qubits[ie])
#         if (gate_type[0][ie] == 0):
#             circ.cx(auxiliary_index, qubits[ie])
#             circ.cx(auxiliary_index, qubits[ie])
            
    # 实现qubit4,5,6上的受控A
    for ie in range(0, len(gate_type[1])):
        if (gate_type[1][ie] == 1):
            circ.cz(auxiliary_index, qubits[ie+3])
#         if (gate_type[1][ie] == 0):
#             circ.cx(auxiliary_index, qubits[ie+3])
#             circ.cx(auxiliary_index, qubits[ie+3])

    # 多qubit,CSWAP的实现
    circ.cswap(0,1,4)
    circ.cswap(0,2,5)
    circ.cswap(0,3,6)
    
 
    #swap test辅助比特右侧H门
    circ.h(auxiliary_index)
    circ.measure(0,0)
    return circ


# #验证函数swap_test()是否正确
# circ1 = swap_test([[0,1,0], [1,0,0]],[1, 2, 3,4,5,6], 0, [[1,1,1],[1,1,1],[1,1,1]],[[1,3,3],[1,2,1],[2,3,1]])
# #电路可视化
# circ1.draw(output='mpl')
# simulator = Aer.get_backend('qasm_simulator')
# job = execute(circ1,simulator,shots=1000)
# result = job.result()
# # Returns counts
# counts = result.get_counts(circ1)
# plot_histogram(counts)
# # 统计辅助比特测得1的概率
# overall_sum_2=0
# if ('0000001' in counts.keys()):
#     m_sum = float(counts["0000001"])/1000
# else:
#     m_sum = 0
                
# overall_sum_2+=(1-2*m_sum)
# print('the final outcome of swap_test: '+str(overall_sum_2))

（2）定义计算,的相关函数

#Implements the entire cost function on the quantum circuit (sampling, 1000 shots)


def calculate_had_test(parameters,count_type):

    global opt
    
    #初始化参数，用于存储<psi|psi>的计算结果
    overall_sum_1 = 0
    # 对传过来的列表x0进行切片
    parameters = [parameters[0:3], parameters[3:6], parameters[6:9]]
    count_type = [count_type[0:3], count_type[3:6], count_type[6:9]]


    for i in range(0, len(gate_set)):
        for j in range(0, len(gate_set)):
            backend = Aer.get_backend('aer_simulator')
            
            # 统计两个受控门的系数乘积  
            multiply = coefficient_set[i]*coefficient_set[j]
            
            # 接收had_test()函数返回的电路
            circ1 = had_test([gate_set[i], gate_set[j]], [1, 2, 3], 0, parameters,count_type)
            simulator = Aer.get_backend('qasm_simulator')
            job = execute(circ1,simulator,shots=100)
            result = job.result()
            # Returns counts
            outputstate = result.get_counts(circ1)
            
            # 统计辅助比特测得1的概率
            if ('0001' in outputstate.keys()):
                m_sum = float(outputstate["0001"])/100
            else:
                m_sum = 0
                
            # 内积
            overall_sum_1+=(abs(multiply))*(1-2*m_sum)
#     print('calculate_had_test:   '+str(overall_sum_1))
    return overall_sum_1





#计算内积的模方


def calculate_swap_test(parameters,count_type):

    global opt
    overall_sum_2 = 0
    parameters = [parameters[0:3], parameters[3:6], parameters[6:9]]
    count_type = [count_type[0:3], count_type[3:6], count_type[6:9]]


    for i in range(0, len(gate_set)):
        for j in range(0, len(gate_set)):
            multiply = coefficient_set[i]*coefficient_set[j]
            mult = 1     
            backend = Aer.get_backend('aer_simulator')
            circ = swap_test([gate_set[i], gate_set[j]],[1, 2, 3,4,5,6], 0, parameters,count_type)
            job = execute(circ, backend,shots=100)
            result = job.result()
            # Returns counts
            outputstate = result.get_counts(circ)
         
            #计算辅助比特测得1的概率

            if ('0000001' in outputstate.keys()):
                m_sum = float(outputstate["0000001"])/100
            else:
                m_sum = 0

            mult = mult*(1-2*m_sum)
            overall_sum_2+=(abs(multiply))*mult
#             print('calculate_swap_test:   '+str(overall_sum_2))
    return overall_sum_2

（3） 代价函数数值计算

def calculate_cost_function2(parameters,count_type):

    global opt

    # 获取计算cost所需数值
    overall_sum_1=calculate_had_test(parameters,count_type)
    overall_sum_2=calculate_swap_test(parameters,count_type)
    
    # 计算代价函数
    outcome= 1-float(overall_sum_2/overall_sum_1)


    # 确保代价函数的数值大于等于0
    final_outcome=abs(outcome)

    # 返回代价函数值
    return final_outcome 

（4）主函数定义，调参，改变ansatz组成以及返回代价函数的最小值、最小值所对应的参数以及ansatz组成

# 主函数
coefficient_set = [0.3, 0.5, 0.2]
gate_set = [[0, 0, 0], [0, 1, 0],[1,0,0]]

for i in range(0,1):
    #随机生成9个参数，用于搭建ansatz,用列表count_type来存储生成的随机数
    count_type = []
    for iz in range(0, 9):
        #生成3个随机数,随机数为1:旋转门x,随机数为2：旋转门y,随机数为3：旋转门z
        count_gate= random.randint(1, 3)
        count_type.append(count_gate)

                 
                  
    # 创建列表存储代价函数值,用于存储同一种ansatz下不同参数对应的代价函数值
  
    cost_function_value2 = []
    
    # 定义字典，存储随机生成的parameters,最小代价函数值和最小代价函数所对应的parameters构成了一对键值对
   
    her_parameters2 = {}
    
    for j in range(0,100):
        x0=[float(random.randint(0,3000))/1000 for i in range(0, 9)]

        
        out2 = minimize(calculate_cost_function2, x0,count_type, method="COBYLA", options={'maxiter':10})

        cost_function_value2.append(out2['fun'])
        # 记录代价函数值和参数之间的对应关系
        her_parameters2[out2['fun']]=out2['x']
    

      
        

        
    #当前ansatz下代价函数的最小值，将代价函数值从小到大排序，输出最小值
    cost_function_value2.sort()

    print('VQLS代价函数的收敛效果，归一化：')
    print(cost_function_value2)

    # 代价函数最小值
    print(cost_function_value2[0])

    # 输出代价函数最小时所对应的参数
    print(her_parameters2[cost_function_value2[0]])

    

 数值结果实例：

VQLS代价函数的收敛效果，归一化：
[0.0010570824524314126, 0.0018952854773752748, 0.0027380590203833766, 0.003416856492027387, 0.0038834951456307998, 0.00392376681614337, 0.005199306759099187, 0.011776251226692791, 0.025481019240769842, 0.035673187571921616, 0.037225792155993864, 0.03749147920927087, 0.047769360269360406, 0.04978354978354982, 0.0550791556728234, 0.058047493403694306, 0.06104972375690643, 0.06764168190127973, 0.06896551724137934, 0.09025174076057862, 0.09208972845336472, 0.09312977099236641, 0.09620872401141456, 0.10249307479224345, 0.10842645726366673, 0.11380267474135763, 0.1204761904761904, 0.12392179413456028, 0.12450787401574825, 0.126213592233009, 0.13102004187855232, 0.13670999551770535, 0.1495031616982836, 0.15706346885054134, 0.1675764192139737, 0.1771561771561776, 0.18119975262832422, 0.19633867276887895, 0.19959879638916767, 0.20328425821064555, 0.21076233183856496, 0.2112903225806454, 0.21566401816118053, 0.23907225691347023, 0.24129651860744306, 0.24899008446566284, 0.25098039215686263, 0.25314685314685337, 0.25764192139737996, 0.2589531680440771, 0.2602827763496144, 0.26580921757770626, 0.2848275862068965, 0.28923476005188076, 0.29241877256317683, 0.3096418732782369, 0.3221621621621621, 0.32710280373831846, 0.32868904369854335, 0.3407821229050281, 0.3444380069524914, 0.38206896551724123, 0.39854743531547887, 0.39929605005866253, 0.4001141552511416, 0.4024810475534113, 0.4057142857142857, 0.41701011614837014, 0.4305335446381404, 0.4507678410117434, 0.4554189085318985, 0.4849215922798552, 0.4908448442017349, 0.5017041581458759, 0.5276497695852533, 0.5356923076923078, 0.5357615894039736, 0.5577689243027888, 0.5602294455066923, 0.5635359116022098, 0.5876502279320348, 0.5915682967959528, 0.5962301587301587, 0.5972927241962774, 0.6029825548677548, 0.6033376123234918, 0.6115107913669066, 0.6231721034870643, 0.6740027510316369, 0.6787741203178208, 0.7085173501577287, 0.716549295774648, 0.7179814755669116, 0.7194875252865813, 0.7740569668976136, 0.8255905511811024, 0.8273195876288659, 0.828113063407181, 0.8533402922755742, 0.8668639053254439]
0.0010570824524314126
[3.441 2.623 3.138 3.402 2.655 0.376 1.67  1.736 3.053]

​

结果说明该ansatz下，100个不同的参数，代价函数收敛值最小为0.001，参数的选取为[3.441,2.623....].