从零实现VQE：量子态+Hamiltonian+优化器

一、项目初始化

# 创建项目结构mkdir-p quantum_lab/{core,experiments,ui}cdquantum_lab# 创建requirements.txtcat>requirements.txt<<'EOF' streamlit>=1.28.0 numpy>=1.24.0 scipy>=1.11.0 matplotlib>=3.7.0 pandas>=2.0.0 plotly>=5.17.0 seaborn>=0.12.0 tqdm>=4.65.0 joblib>=1.3.0 EOF# 创建主文件touch__init__.pytouchmain.py

二、核心后端实现

2.1 量子态模拟器基类

# core/backend.pyfromabcimportABC,abstractmethodfromtypingimportList,Tuple,Dict,AnyimportnumpyasnpclassQuantumBackend(ABC):"""量子计算后端抽象基类"""@abstractmethoddefinitialize(self,n_qubits:int)->None:"""初始化量子寄存器"""pass@abstractmethoddefapply_gate(self,gate_name:str,targets:List[int],controls:List[int]=None,params:Dict=None)->None:"""应用量子门"""pass@abstractmethoddefget_statevector(self)->np.ndarray:"""获取当前量子态"""pass@abstractmethoddefmeasure(self,shots:int=1024)->Dict[str,int]:"""测量量子态"""pass@abstractmethoddefreset(self)->None:"""重置到|0>态"""pass@abstractmethoddefexpectation_value(self,pauli_string:str)->float:"""计算Pauli算符的期望值"""pass

2.2 Statevector模拟器实现

# core/statevector_backend.pyimportnumpyasnpfromtypingimportList,Dict,Anyimportitertoolsfrom.backendimportQuantumBackendclassStatevectorBackend(QuantumBackend):"""纯Python实现的量子态矢量模拟器"""# Pauli矩阵定义PAULI_MATRICES={'I':np.array([[1,0],[0,1]],dtype=complex),'X':np.array([[0,1],[1,0]],dtype=complex),'Y':np.array([[0,-1j],[1j,0]],dtype=complex),'Z':np.array([[1,0],[0,-1]],dtype=complex),}# 基础量子门矩阵GATE_MATRICES={'H':np.array([[1,1],[1,-1]],dtype=complex)/np.sqrt(2),'X':np.array([[0,1],[1,0]],dtype=complex),'Y':np.array([[0,-1j],[1j,0]],dtype=complex),'Z':np.array([[1,0],[0,-1]],dtype=complex),'S':np.array([[1,0],[0,1j]],dtype=complex),'T':np.array([[1,0],[0,np.exp(1j*np.pi/4)]],dtype=complex),'CNOT':None,# 特殊处理'SWAP':None,# 特殊处理}def__init__(self,n_qubits:int):self.n_qubits=n_qubits self.state=np.zeros(2**n_qubits,dtype=complex)self.state[0]=1.0# |00...0⟩self.gate_history=[]definitialize(self,n_qubits:int=None)->None:"""初始化量子寄存器"""ifn_qubits:self.n_qubits=n_qubits self.state=np.zeros(2**self.n_qubits,dtype=complex)self.state[0]=1.0self.gate_history=[]def_tensor_product(self,matrices:List[np.ndarray])->np.ndarray:"""计算张量积"""result=matrices[0]formatinmatrices[1:]:result=np.kron(result,mat)returnresultdef_apply_single_qubit_gate(self,gate:np.ndarray,target:int)->None:"""应用单量子比特门"""# 构造完整的变换矩阵matrices=[self.PAULI_MATRICES['I']]*self.n_qubits matrices[target]=gate full_matrix=self._tensor_product(matrices)# 应用变换self.state=full_matrix @ self.statedef_apply_controlled_gate(self,gate:np.ndarray,control:int,target:int)->None:"""应用受控门（CNOT, CZ等）"""# 构建受控门的矩阵表示dim=2**self.n_qubits matrix=np.eye(dim,dtype=complex)# 遍历所有基态foriinrange(dim):# 将索引转换为二进制字符串bits=format(i,f'0{self.n_qubits}b')# 检查控制位是否为1ifbits[control]=='1':# 计算目标位翻转后的索引target_bit=bits[target]new_target_bit='0'iftarget_bit=='1'else'1'# 构造新的二进制字符串new_bits=list(bits)new_bits[target]=new_target_bit new_bits_str=''.join(new_bits)# 计算新索引j=int(new_bits_str,2)# 对于CNOT，就是简单的交换振幅matrix[i,i]=0matrix[i,j]=1# 应用变换self.state=matrix @ self.statedef_rotation_gate(self,axis:str,angle:float)->np.ndarray:"""生成旋转门矩阵"""ifaxis=='x':returnnp.array([[np.cos(angle/2),-1j*np.sin(angle/2)],[-1j*np.sin(angle/2),np.cos(angle/2)]],dtype=complex)elifaxis=='y':returnnp.array([[np.cos(angle/2),-np.sin(angle/2)],[np.sin(angle/2),np.cos(angle/2)]],dtype=complex)elifaxis=='z':returnnp.array([[np.exp(-1j*angle/2),0],[0,np.exp(1j*angle/2)]],dtype=complex)else:raiseValueError(f"Unknown rotation axis:{axis}")defapply_gate(self,gate_name:str,targets:List[int],controls:List[int]=None,params:Dict=None)->None:"""应用量子门"""controls=controlsor[]params=paramsor{}# 记录门操作self.gate_history.append({'gate':gate_name,'targets':targets,'controls':controls,'params':params})# 处理旋转门ifgate_name.startswith('R'):axis=gate_name[1].lower()angle=params.get('theta',0.0)gate_matrix=self._rotation_gate(axis,angle)fortargetintargets:self._apply_single_qubit_gate(gate_matrix,target)# 处理标准单量子比特门elifgate_nameinself.GATE_MATRICESandself.GATE_MATRICES[gate_name]isnotNone:gate_matrix=self.GATE_MATRICES[gate_name]fortargetintargets:self._apply_single_qubit_gate(gate_matrix,target)# 处理CNOT门elifgate_name=='CNOT':iflen(controls)!=1orlen(targets)!=1:raiseValueError("CNOT gate requires exactly one control and one target")self._apply_controlled_gate(self.GATE_MATRICES['X'],controls[0],targets[0])# 处理CZ门elifgate_name=='CZ':iflen(controls)!=1orlen(targets)!=1:raiseValueError("CZ gate requires exactly one control and one target")self._apply_controlled_gate(self.PAULI_MATRICES['Z'],controls[0],targets[0])else:raiseValueError(f"Unsupported gate:{gate_name}")defget_statevector(self)->np.ndarray:"""获取当前量子态"""returnself.state.copy()defget_probabilities(self)->np.ndarray:"""获取测量概率分布"""returnnp.abs(self.state)**2defmeasure(self,shots:int=1024)->Dict[str,int]:"""测量量子态"""probs=self.get_probabilities()counts={}# 概率采样samples=np.random.choice(len(probs),size=shots,p=probs)# 统计结果forsampleinsamples:bitstring=format(sample,f'0{self.n_qubits}b')counts[bitstring]=counts.get(bitstring,0)+1returncountsdefexpectation_value(self,pauli_string:str)->float:"""计算Pauli算符的期望值"""# 解析Pauli字符串，如 "X0 Z1"terms=pauli_string.split()coeff=1.0iflen(terms)==0:return1.0# 检查是否有系数ifterms[0].replace('.','').replace('-','').isdigit():coeff=float(terms[0])terms=terms[1:]# 构建Pauli算符的矩阵表示matrices=[self.PAULI_MATRICES['I']]*self.n_qubitsforterminterms:pauli=term[0]qubit=int(term[1:])matrices[qubit]=self.PAULI_MATRICES[pauli]# 计算张量积pauli_matrix=self._tensor_product(matrices)# 计算期望值expectation=np.real(np.vdot(self.state,pauli_matrix @ self.state))returncoeff*expectationdefreset(self)->None:"""重置到|0>态"""self.initialize(self.n_qubits)def__str__(self)->str:"""字符串表示"""returnf"StatevectorBackend(n_qubits={self.n_qubits})"

2.3 可视化工具

# core/visualization.pyimportnumpyasnpimportmatplotlib.pyplotaspltfrommatplotlib.patchesimportFancyArrowPatchfrommpl_toolkits.mplot3dimportproj3dimportplotly.graph_objectsasgoclassQuantumVisualizer:"""量子态可视化工具"""@staticmethoddefplot_statevector(statevector:np.ndarray,title:str="Quantum State"):"""绘制量子态振幅和相位"""fig,(ax1,ax2)=plt.subplots(1,2,figsize=(12,4))n_qubits=int(np.log2(len(statevector)))basis_states=[format(i,f'0{n_qubits}b')foriinrange(len(statevector))]# 振幅图amplitudes=np.abs(statevector)ax1.bar(basis_states,amplitudes**2,alpha=0.7,color='blue')ax1.set_xlabel('Basis State')ax1.set_ylabel('Probability')ax1.set_title(f'{title}- Probabilities')ax1.tick_params(axis='x',rotation=45)# 相位图phases=np.angle(statevector)ax2.bar(basis_states,phases,alpha=0.7,color='red')ax2.set_xlabel('Basis State')ax2.set_ylabel('Phase (radians)')ax2.set_title(f'{title}- Phases')ax2.tick_params(axis='x',rotation=45)plt.tight_layout()returnfig@staticmethoddefplot_bloch_vector(theta:float,phi:float,title:str="Bloch Sphere"):"""在Bloch球上绘制量子态"""# 计算笛卡尔坐标x=np.sin(theta)*np.cos(phi)y=np.sin(theta)*np.sin(phi)z=np.cos(theta)# 创建Bloch球fig=plt.figure(figsize=(8,8))ax=fig.add_subplot(111,projection='3d')# 绘制球体u=np.linspace(0,2*np.pi,100)v=np.linspace(0,np.pi,100)xs=np.outer(np.cos(u),np.sin(v))ys=np.outer(np.sin(u),np.sin(v))zs=np.outer(np.ones(np.size(u)),np.cos(v))ax.plot_surface(xs,ys,zs,color='b',alpha=0.1)# 绘制坐标轴ax.quiver(0,0,0,1.5,0,0,color='r',arrow_length_ratio=0.1,label='X')ax.quiver(0,0,0,0,1.5,0,color='g',arrow_length_ratio=0.1,label='Y')ax.quiver(0