<think>错误分析:
在`create_circuit`方法中,当添加双量子比特门(如CNOT或CZ)后,我们尝试添加双量子比特去极化噪声。错误信息指出:`cirq.depolarize`期望作用于1个量子比特,但我们传入了2个量子比特。
这是因为`cirq.depolarize`是一个单量子比特噪声通道。对于双量子比特,我们需要使用`cirq.depolarize`的双量子比特版本,即使用`cirq.two_qubit_depolarize`。
修改方案:
1. 对于单量子比特门后的噪声,使用`cirq.depolarize(p)`作用于单个量子比特。
2. 对于双量子比特门后的噪声,使用`cirq.two_qubit_depolarize(p)`作用于两个量子比特。
另外,注意`cirq.two_qubit_depolarize`的参数p是总错误率,其含义与单量子比特去极化噪声的参数p不同。在双量子比特去极化信道中,总错误率p被分配到15个Pauli错误上(除了单位算子)。因此,双量子比特去极化信道的定义与之前TensorCircuit中的定义(p/15)是兼容的,但这里我们直接使用给定的p值。
但是,我们需要注意:在TensorCircuit中,双量子比特去极化噪声的信道定义为:p/15 * (15种Pauli错误)。而Cirq的`two_qubit_depolarize`函数也是用相同的定义:以概率p应用一个随机的双量子比特Pauli错误(不包括恒等操作)。因此,我们可以直接使用相同的p值。
因此,修改噪声添加方式:
- 对于单比特门:添加`cirq.depolarize(self.single_qubit_depolarizing_p).on(qubit)`
- 对于双比特门:添加`cirq.two_qubit_depolarize(self.two_qubit_depolarizing_p).on(qubit1, qubit2)`
同时,比特翻转噪声是单量子比特噪声,我们使用`cirq.bit_flip(p).on(qubit)`。
另外,注意在双量子比特门后添加噪声时,我们使用`two_qubit_depolarize`,并且作用于两个量子比特。
修改`create_circuit`方法中的噪声部分:
原代码:
if gate[0] == "cx":
...
if self.noise:
circuit.append(cirq.depolarize(self.two_qubit_depolarizing_p).on(qubits[gate[1]], qubits[gate[2]]))
修改为:
if gate[0] == "cx":
...
if self.noise:
# 使用双量子比特去极化噪声
circuit.append(cirq.two_qubit_depolarize(self.two_qubit_depolarizing_p).on(qubits[gate[1]], qubits[gate[2]]))
同样,对CZ门做同样修改。
对于单量子比特旋转门,噪声添加方式不变(使用单量子比特去极化)。
另外,在电路最后添加比特翻转噪声时,我们使用单量子比特比特翻转噪声,所以保持不变。
同时,注意哈密顿量构建部分:原代码中构建的哈密顿量是一个`cirq.PauliSum`对象,但在分解时,我们使用了`for term in self.hamiltonian_cirq`,这实际上会遍历每个单独的PauliString项(如果Ham是由多个PauliString相加而成)。但是,Cirq的`PauliSum`对象实际上是一个单项的集合,我们可以直接使用`self.hamiltonian_cirq`的`expectation_from_state_vector`或通过密度矩阵计算期望值。
但是,在混合态计算中,我们需要将哈密顿量分解为单独的PauliString项。我们可以使用`self.hamiltonian_cirq`的展开形式。
我们可以这样做:在初始化时,将哈密顿量分解为单独的PauliString列表。因为`PauliSum`是一个可迭代对象,所以我们可以将其转换为列表。
修改`build_heisenberg_hamiltonian`方法,将构建的哈密顿量分解为单独的PauliString项,并存储为列表。
但是注意:构建哈密顿量时,我们使用了多个PauliString相加。Cirq的`PauliSum`对象会自动合并同类项,所以我们可以直接使用`self.hamiltonian_cirq`的展开(通过`.terms`属性)?
不过,在计算期望值时,我们需要对哈密顿量的每一项单独计算。我们可以使用`self.hamiltonian_cirq`的`__iter__`方法进行迭代,或者使用`cirq.PauliSum.collect()`方法展开为单项式列表。
实际上,我们可以使用`self.hamiltonian_cirq.to_z_basis_terms()`来获取每一项,但这里我们直接使用`self.hamiltonian_cirq`的迭代(它已经是一个单项式的和)。
但是,在构建哈密顿量时,我们使用的是`cirq.PauliSum`的加法,所以它内部是一个单项式的列表。因此,在计算期望值时,我们可以直接迭代`self.hamiltonian_cirq`,每次迭代得到一个项(`PauliString`)和它的系数(注意:在构建时,每个项都乘以了系数,所以每个PauliString项已经带有了系数)。
因此,在计算期望值时,我们不需要再单独乘以系数。
但是,在能量计算部分,我们目前对每个项的处理是:
混合态:energy += trace_val # 这里trace_val已经包含了该项的系数
纯态:energy += exp_val.real # exp_val是该项的期望值,已经包含系数
因此,不需要修改。
另外,注意在纯态计算中,我们使用了`term.expectation_from_state_vector`,这里term是一个`PauliString`,它已经包含了系数。所以计算正确。
我们还需要注意:在构建哈密顿量时,我们使用了`ham += ...`,这会产生一个`PauliSum`。但是,当我们在最后返回`ham`时,它是一个`PauliSum`对象。在计算期望值时,迭代这个`PauliSum`对象会得到每一项(即每个PauliString乘以系数)。
所以,哈密顿量构建部分不需要修改。
修改后的代码:
在`create_circuit`方法中,将双量子比特噪声改为`cirq.two_qubit_depolarize`。
同时,由于错误信息还指出`qubits[gate[1]]`和`qubits[gate[2]]`是两个量子比特,所以我们使用`on`方法作用于两个量子比特。
另外,注意`two_qubit_depolarize`是一个双量子比特门操作,所以它的`on`方法需要两个量子比特。
修改后,重新运行。
此外,还有一个潜在问题:在`compute_energy`方法中,混合态模拟时,我们使用`result.final_density_matrix`,然后计算`trace_val = np.trace(density_matrix @ op_matrix).real`。这里`op_matrix`是PauliString的矩阵表示,但注意PauliString可能作用于多个量子比特,且可能不是全空间(只作用在部分量子比特上)。我们应该使用整个希尔伯特空间上的矩阵表示。
但是,在构建PauliString时,我们使用的是具体的量子比特。在计算矩阵表示时,`op_matrix`是作用在指定量子比特上的矩阵(形状为(2^k, 2^k)),而密度矩阵的维度是(2^n, 2^n)。为了使乘法正确,我们需要将PauliString的矩阵扩展到整个系统空间(即其他量子比特上作用单位矩阵)。然而,Cirq的`PauliString.matrix()`方法返回的矩阵仅作用于它自己的量子比特,并不是整个系统。因此,我们需要将密度矩阵和PauliString都表示为整个系统空间的算子。
这是一个复杂的问题。另一种方法是使用`qsimcirq`提供的期望值计算函数,但qsimcirq目前没有直接提供计算混合态期望值的函数。我们可以考虑使用`cirq.DensityMatrixSimulator`来模拟,然后利用`cirq.expectation`函数计算期望值,但这样效率可能较低。
或者,我们可以在计算期望时,将密度矩阵和PauliString都表示为整个希尔伯特空间上的矩阵,然后直接计算Tr(ρ·H)。但是,当量子比特数较大时(例如12个量子比特),密度矩阵大小为4096x4096,整个哈密顿量有O(n)项(对于海森堡模型,有3n项),计算量很大(4096x4096矩阵乘法需要约10^9次运算,乘以3*12=36次,总运算量约为3.6e10,这在一台普通计算机上是可以接受的,但可能较慢)。
不过,我们注意到海森堡模型的哈密顿量是稀疏的,我们可以利用稀疏矩阵乘法。但密度矩阵是稠密的。因此,计算整个哈密顿量的期望值需要多次矩阵乘法,每次涉及一个4096x4096的矩阵。
另一种方法:我们可以使用`qsimcirq`的混合态模拟功能,但直接计算期望值需要自己实现。
考虑到性能,我们可以尝试将哈密顿量预先计算为一个完整的矩阵(当n=12时,矩阵大小为4096x4096,占用内存4096^2 * 8(复数)= 134217728字节≈134MB,这对于现代计算机来说是可以接受的)。然后,计算Tr(ρ·H)就变成一次矩阵乘法和取迹。
修改混合态能量计算方式:
在初始化时,预先计算整个哈密顿量矩阵(作为稀疏矩阵或稠密矩阵)。由于12量子比特的矩阵大小为4096x4096,我们可以存储为稠密矩阵(134MB)。然后,在混合态计算中,直接计算`np.trace(density_matrix @ h_matrix)`。
同样,在纯态计算中,我们也可以计算`<ψ|H|ψ>`,即`state_vector.conj() @ (h_matrix @ state_vector)`。
这样修改可以避免分解哈密顿量,并且计算更直接。
但是,我们需要注意:构建整个哈密顿量矩阵需要消耗一定的时间和内存。对于12量子比特,这是可行的,但对于更大的系统(例如30量子比特)则不可行。所以这里我们只针对12量子比特使用。
修改步骤:
1. 在`__init__`中,构建哈密顿量矩阵(稠密矩阵或稀疏矩阵)。由于我们已经在`build_heisenberg_hamiltonian`中构建了`PauliSum`,我们可以将其转换为矩阵。
2. 使用`cirq.unitary`函数将`PauliSum`转换为矩阵(注意:`PauliSum`实现了`_unitary_`方法?)。但是,`PauliSum`可能没有直接提供矩阵表示。我们可以使用`cirq.dense`模块中的`pauli_sum_to_dense`函数?或者,我们可以使用`cirq.unitary`只能用于单一门,不能用于整个哈密顿量。
3. 我们可以通过将哈密顿量作用在单位矩阵上,然后线性组合得到整个矩阵?但这样计算量很大。
4. 另一种方法:我们遍历哈密顿量中的每个PauliString,计算每个PauliString的矩阵(使用`cirq.unitary`或`PauliString.matrix()`方法),然后乘以系数并累加。注意,每个PauliString的矩阵是作用在它自己的量子比特集的子空间上,我们需要扩展到整个系统(即其他量子比特上为单位矩阵)。
5. 但是,由于我们已经有了所有量子比特的列表,我们可以为每个PauliString计算它在整个系统上的矩阵表示。我们可以使用`cirq.kron`来扩展,但这样很麻烦。Cirq提供了`cirq.pauli_string_expectation`,但没有直接提供矩阵扩展。
6. 一个简单的方法:使用`cirq.density_matrix_simulator`的`simulate_expectation_values`方法,但这样会重新模拟电路,而我们已经在混合态模拟中得到了密度矩阵。
7. 我们决定:在初始化时,预先计算整个哈密顿量的矩阵(稠密矩阵)。使用以下函数:
h_matrix = self.hamiltonian_cirq.matrix(qubit_order=qubits)
其中qubits是量子比特列表(顺序很重要)。
8. 但是,`cirq.PauliSum`有一个`matrix`方法,它要求指定`qubits_that_should_be_present`(所有量子比特)和`qubit_order`。
9. 因此,我们在`__init__`中构建哈密顿量矩阵:
qubits = cirq.LineQubit.range(self.n_qubit)
self.h_matrix = self.hamiltonian_cirq.matrix(qubit_order=qubits)
10. 在计算能量时:
纯态:energy = (np.vdot(state_vector, self.h_matrix @ state_vector)).real
混合态:energy = np.trace(density_matrix @ self.h_matrix).real
这样修改后,能量计算将变得简单直接,而且不需要遍历每一项。
但是,注意:哈密顿量矩阵的大小为`(2**n, 2**n)`,当n=12时,矩阵元素个数为2^24=16777216,每个元素是复数(16字节),总内存为16777216*16=268435456字节≈268MB。如果n=16,则矩阵大小为(65536,65536),元素个数2^32≈4e9,需要32GB内存(每个元素16字节)。所以对于12量子比特(268MB)是可行的。
因此,我们修改`__init__`,在构建哈密顿量后,计算其矩阵形式(稠密矩阵)并存储。
修改步骤:
在`__init__`方法中:
qubits = cirq.LineQubit.range(self.n_qubit)
self.h_matrix = self.hamiltonian_cirq.matrix(qubit_order=qubits)
然后,修改`compute_energy`方法:
if self.noise:
... # 模拟得到密度矩阵
energy = np.trace(density_matrix @ self.h_matrix).real
else:
... # 模拟得到态矢量
energy = (np.vdot(state_vector, self.h_matrix @ state_vector)).real
注意:在纯态计算中,我们也可以使用`np.vdot`计算,这样避免了显式地构建`<ψ|`,而是使用`np.vdot(state_vector, ...)`。
但是,由于`state_vector`是一维数组(列向量),我们可以将其视为列向量,则`<ψ|H|ψ>`等价于`state_vector.conj().dot(self.h_matrix.dot(state_vector))`,即`np.vdot(state_vector, self.h_matrix.dot(state_vector))`(因为`np.vdot`会自动取第一个参数的共轭)。
或者直接:`(state_vector.conj().T @ self.h_matrix @ state_vector).real`,但这样可能更慢。
使用`np.vdot`更高效:`np.vdot(a, b)`是`a.conj().dot(b)`。
因此,纯态能量计算:`energy = np.vdot(state_vector, self.h_matrix.dot(state_vector)).real`
混合态:`energy = np.trace(density_matrix @ self.h_matrix).real`
注意:`density_matrix`和`self.h_matrix`都是二维数组。
由于我们使用了整个哈密顿量的矩阵,所以避免了分解。这样修改后,代码更简洁,且对于12量子比特是可行的。
下面,我们整合这些修改:
1. 在`__init__`中存储量子比特列表并计算哈密顿量矩阵。
2. 修改噪声添加部分,使用`cirq.two_qubit_depolarize`。
3. 修改能量计算方式为直接使用哈密顿量矩阵。
由于我们使用了整个矩阵,所以不再需要原来计算期望值的循环。
修改后的代码:
注意:在`build_heisenberg_hamiltonian`方法中,我们使用了`qubits`变量,但在`__init__`中,我们需要再次使用量子比特列表来计算矩阵。因此,我们在`__init__`中构建量子比特列表,并传递给`build_heisenberg_hamiltonian`,或者在`build_heisenberg_hamiltonian`中构建量子比特列表,然后返回并存储。
修改`__init__`:
self.qubits = cirq.LineQubit.range(self.n_qubit)
self.hamiltonian_cirq = self.build_heisenberg_hamiltonian(self.qubits)
# 计算哈密顿量矩阵
self.h_matrix = self.hamiltonian_cirq.matrix(qubit_order=self.qubits)
相应地,修改`build_heisenberg_hamiltonian`方法,接收`qubits`参数。
另外,在`compute_energy`方法中,纯态计算时,注意`state_vector`的顺序与`qubit_order`一致即可(因为我们使用相同的量子比特顺序构建哈密顿量矩阵)。
最后,由于我们修改了能量计算方式,在混合态计算中,我们使用的是整个密度矩阵与整个哈密顿量矩阵的乘积的迹。这种方法在量子比特数较多时内存消耗大,但12量子比特可以接受。
修改后的完整代码:
</think>### 错误分析
根据错误信息,问题在于尝试将单量子比特噪声通道 (`cirq.depolarize`) 应用到两个量子比特上。在 Cirq 中:
- `cirq.depolarize` 是单量子比特噪声通道
- `cirq.two_qubit_depolarize` 是双量子比特噪声通道
### 修改后的完整代码
```python
import os
import numpy as np
import tensorflow as tf
from multiprocessing import Pool
import matplotlib.pyplot as plt
import time
import utils
import argparse
import cirq
import qsimcirq
class VqeTrainerNew:
def __init__(self, n_cir_parallel, n_runs, max_iteration, n_qubit, hamiltonian, noise_param=None):
self.n_qubit = n_qubit
self.max_iteration = max_iteration
self.n_cir_parallel = n_cir_parallel
self.n_runs = n_runs
self.hamiltonian_ = hamiltonian
# 创建量子比特列表
self.qubits = [cirq.LineQubit(i) for i in range(n_qubit)]
# 构建海森堡模型哈密顿量
self.hamiltonian_cirq = self.build_heisenberg_hamiltonian()
# 计算哈密顿量矩阵(用于高效计算期望值)
self.h_matrix = self.hamiltonian_cirq.matrix(qubit_order=self.qubits)
# 噪声参数设置
if noise_param is None:
self.noise = False
else:
self.noise = True
self.two_qubit_depolarizing_p = noise_param['two_qubit_channel_depolarizing_p']
self.single_qubit_depolarizing_p = noise_param['single_qubit_channel_depolarizing_p']
self.bit_flip_p = noise_param['bit_flip_p']
# qsim模拟器配置(使用GPU)
options = {'t': 16, 'f': 2, 'g': 0} # t: 线程数, f: GPU模式
self.simulator = qsimcirq.QSimSimulator(options)
def build_heisenberg_hamiltonian(self):
"""构建海森堡模型哈密顿量(Cirq PauliSum格式)"""
ham = 0.0
# ZZ相互作用
for i in range(self.n_qubit - 1):
ham += self.hamiltonian_['hzz'] * cirq.Z(self.qubits[i]) * cirq.Z(self.qubits[i+1])
if self.hamiltonian_['pbc'] and self.n_qubit > 2:
ham += self.hamiltonian_['hzz'] * cirq.Z(self.qubits[-1]) * cirq.Z(self.qubits[0])
# XX相互作用
for i in range(self.n_qubit - 1):
ham += self.hamiltonian_['hxx'] * cirq.X(self.qubits[i]) * cirq.X(self.qubits[i+1])
if self.hamiltonian_['pbc'] and self.n_qubit > 2:
ham += self.hamiltonian_['hxx'] * cirq.X(self.qubits[-1]) * cirq.X(self.qubits[0])
# YY相互作用
for i in range(self.n_qubit - 1):
ham += self.hamiltonian_['hyy'] * cirq.Y(self.qubits[i]) * cirq.Y(self.qubits[i+1])
if self.hamiltonian_['pbc'] and self.n_qubit > 2:
ham += self.hamiltonian_['hyy'] * cirq.Y(self.qubits[-1]) * cirq.Y(self.qubits[0])
# 外场项
for i in range(self.n_qubit):
if self.hamiltonian_['hz'] != 0:
ham += self.hamiltonian_['hz'] * cirq.Z(self.qubits[i])
if self.hamiltonian_['hx'] != 0:
ham += self.hamiltonian_['hx'] * cirq.X(self.qubits[i])
if self.hamiltonian_['hy'] != 0:
ham += self.hamiltonian_['hy'] * cirq.Y(self.qubits[i])
return ham
def create_circuit(self, param, structure):
"""构建Cirq量子电路(含噪声)"""
circuit = cirq.Circuit()
param_idx = 0
applied_gates = [] # 记录已应用的量子门用于调试
# 主电路
for gate in structure:
gate_type = gate[0]
q1 = self.qubits[gate[1]]
if gate_type == "cx":
q2 = self.qubits[gate[2]]
circuit.append(cirq.CNOT(q1, q2))
applied_gates.append(f"CNOT({q1}, {q2})")
if self.noise:
# 使用正确的双量子比特去极化噪声
circuit.append(cirq.two_qubit_depolarize(self.two_qubit_depolarizing_p).on(q1, q2))
applied_gates.append(f"Depolarize2({q1}, {q2})")
elif gate_type == "cz":
q2 = self.qubits[gate[2]]
circuit.append(cirq.CZ(q1, q2))
applied_gates.append(f"CZ({q1}, {q2})")
if self.noise:
circuit.append(cirq.two_qubit_depolarize(self.two_qubit_depolarizing_p).on(q1, q2))
applied_gates.append(f"Depolarize2({q1}, {q2})")
elif gate_type == "ry":
circuit.append(cirq.ry(param[param_idx]).on(q1))
applied_gates.append(f"Ry({q1})")
if self.noise:
circuit.append(cirq.depolarize(self.single_qubit_depolarizing_p).on(q1))
applied_gates.append(f"Depolarize1({q1})")
param_idx += 1
elif gate_type == "rz":
circuit.append(cirq.rz(param[param_idx]).on(q1))
applied_gates.append(f"Rz({q1})")
if self.noise:
circuit.append(cirq.depolarize(self.single_qubit_depolarizing_p).on(q1))
applied_gates.append(f"Depolarize1({q1})")
param_idx += 1
elif gate_type == "rx":
circuit.append(cirq.rx(param[param_idx]).on(q1))
applied_gates.append(f"Rx({q1})")
if self.noise:
circuit.append(cirq.depolarize(self.single_qubit_depolarizing_p).on(q1))
applied_gates.append(f"Depolarize1({q1})")
param_idx += 1
else:
raise ValueError(f"Invalid gate type: {gate_type}")
# 添加比特翻转噪声
if self.noise:
for q in self.qubits:
circuit.append(cirq.bit_flip(self.bit_flip_p).on(q))
applied_gates.append(f"BitFlip({q})")
# 打印应用的量子门(调试用)
# print("Circuit gates:", "\n".join(applied_gates))
return circuit
def compute_energy(self, param, structure):
"""计算能量期望值(纯态/混合态)"""
circuit = self.create_circuit(param, structure)
if self.noise:
# 混合态模拟(密度矩阵)
result = self.simulator.simulate(circuit)
density_matrix = result.final_density_matrix
# 计算期望值: Tr[ρH]
energy = np.trace(density_matrix @ self.h_matrix).real
noise_status = "With Noise"
else:
# 纯态模拟
result = self.simulator.simulate(circuit)
state_vector = result.state_vector()
# 计算期望值: <ψ|H|ψ>
energy = np.vdot(state_vector, self.h_matrix @ state_vector).real
noise_status = "No Noise"
# 打印每次能量计算结果
# print(f"Computed energy ({noise_status}): {energy:.6f}")
return energy
def parameter_shift_grad(self, param, structure):
"""使用参数移位法计算梯度"""
grad = np.zeros_like(param)
base_energy = self.compute_energy(param, structure)
for i in range(len(param)):
# 创建参数移位拷贝
param_plus = param.copy()
param_minus = param.copy()
param_plus[i] += np.pi / 2
param_minus[i] -= np.pi / 2
# 计算能量差值
e_plus = self.compute_energy(param_plus, structure)
e_minus = self.compute_energy(param_minus, structure)
# 计算梯度分量
grad[i] = (e_plus - e_minus) / 2
return grad
def get_param_num(self, cir):
param_num = 0
for i in range(len(cir)):
if cir[i][0] in ['rx', 'ry', 'rz']:
param_num += 1
return param_num
def train_circuit(self, circuit_and_seed):
single_circuit = circuit_and_seed[0]
seed = circuit_and_seed[1]
np.random.seed(seed)
tf.random.set_seed(seed)
param_num = self.get_param_num(single_circuit[0])
L = single_circuit[1]
par = np.random.normal(loc=0, scale=1 / (8 * (L + 2)), size=param_num)
param = tf.Variable(initial_value=tf.convert_to_tensor(par, dtype=tf.float32))
energy_epoch = []
opt = tf.keras.optimizers.Adam(0.05)
if param_num > 0:
for i in range(self.max_iteration):
# 计算能量和梯度
with tf.GradientTape() as tape:
energy = self.compute_energy(param.numpy(), single_circuit[0])
# 使用参数移位法计算梯度
grad = self.parameter_shift_grad(param.numpy(), single_circuit[0])
# 应用梯度
opt.apply_gradients([(tf.convert_to_tensor(grad), param)])
energy_epoch.append(energy)
print(f"Seed {seed}, Iter {i}: Energy = {energy:.6f}") # 打印每次优化能量
else:
energy = self.compute_energy(param.numpy(), single_circuit[0])
energy_epoch = [energy for _ in range(self.max_iteration)]
return energy, param.numpy(), energy_epoch
def draw(self, loss_list, exp_dir, circuit_id, best_index):
'''绘图函数保持不变'''
epochs = range(1, len(loss_list) + 1)
plt.figure(figsize=(10, 6))
plt.plot(epochs, loss_list, label='Training Loss', marker='o', markersize=1, color='blue')
plt.title('Training Loss Over Epochs')
plt.xlabel('Epochs')
plt.ylabel('Loss')
plt.legend()
plt.grid()
plt.tight_layout()
save_path_img = f'{exp_dir}/circuit_train_curve/'
if not os.path.exists(save_path_img):
os.makedirs(save_path_img)
plt.savefig(save_path_img + f'circuit_{circuit_id}_{best_index}.png')
plt.close()
def batch_train_parallel(self, device_name, task_name, run_id):
start_time = time.time()
exp_dir = f'result/cir_sample/{device_name}_{task_name}/training/run_{run_id}/'
samples = utils.load_pkl(f'{exp_dir}samples.pkl')[:10]
work_queue = []
for i in range(len(samples)):
print(f'circuit id: {i}')
work_queue.extend([[samples[i], j] for j in range(self.n_runs)])
result = []
for work in work_queue:
print(f"Starting training with seed {work[1]} for circuit {work[0][1]}")
temp_train = self.train_circuit(work)
result.append(temp_train)
energy, param, energy_epoch = [], [], []
for part in result:
energy.append(part[0])
param.append(part[1])
energy_epoch.append(part[2])
energy_f, param_f, energy_epoch_f = [], [], []
for i in range(len(samples)):
index0 = i * self.n_runs
index1 = index0 + self.n_runs
best_index = np.argmin(energy[index0:index1]) + index0
energy_f.append(energy[best_index])
param_f.append(param[best_index])
energy_epoch_f.append(energy_epoch[best_index])
self.draw(energy_epoch[best_index], exp_dir, i, best_index)
if not os.path.exists(exp_dir):
os.makedirs(exp_dir)
utils.save_pkl(energy_f, f'{exp_dir}energy.pkl')
utils.save_pkl(param_f, f'{exp_dir}param.pkl')
utils.save_pkl(energy_epoch_f, f'{exp_dir}energy_epoch.pkl')
end_time = time.time()
duration = end_time - start_time
print(f"Run time: {int(duration // 3600)} hours {int((duration % 3600) // 60)} minutes")
utils.save_pkl(duration, exp_dir + 'duration.pkl')
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument("--noise", type=bool, default=True, help="是否考虑噪声")
parser.add_argument("--two_qubit_depolarizing_p", type=float, default=0.01, help="双qubit去极化噪声大小")
parser.add_argument("--single_qubit_depolarizing_p", type=float, default=0.001, help="单qubit去极化噪声大小")
parser.add_argument("--bit_flip_p", type=float, default=0.01, help="比特翻转噪声大小")
args = parser.parse_args()
device_name = 'grid_16q'
task_name = 'Heisenberg_12'
qubit = 12
noise_param = None
if args.noise:
noise_param = {
'two_qubit_channel_depolarizing_p': args.two_qubit_depolarizing_p,
'single_qubit_channel_depolarizing_p': args.single_qubit_depolarizing_p,
'bit_flip_p': args.bit_flip_p
}
hamiltonian = {
'pbc': True,
'hzz': 1, 'hxx': 1, 'hyy': 1,
'hx': 0, 'hy': 0, 'hz': 1
}
trainer = VqeTrainerNew(
n_cir_parallel=1,
n_runs=1,
max_iteration=1000,
n_qubit=qubit,
hamiltonian=hamiltonian,
noise_param=noise_param
)
all_start_time = time.time()
run_id = 0
trainer.batch_train_parallel(device_name, task_name, run_id)
all_end_time = time.time()
print(f"Total time: {all_end_time - all_start_time:.2f} seconds")
```
### 关键修改说明:
1. **噪声通道修复**:
- 使用 `cirq.two_qubit_depolarize` 替代 `cirq.depolarize` 用于双量子比特门
- 单量子比特门继续使用 `cirq.depolarize`
- 添加了应用门调试信息,便于验证
2. **量子比特管理优化**:
- 在初始化时创建 `self.qubits` 列表
- 所有量子门操作引用此列表中的量子比特对象
- 避免在每次电路构建时重复创建量子比特
3. **哈密顿量计算优化**:
- 提前计算哈密顿量矩阵 `self.h_matrix`
- 纯态能量计算:`<ψ|H|ψ> = np.vdot(state_vector, H @ state_vector)`
- 混合态能量计算:`Tr(ρH) = np.trace(density_matrix @ H)`
- 无需遍历Pauli项,计算更高效
4. **调试信息增强**:
- 添加每个种子的训练开始提示
- 打印每次优化的能量值
- 可选打印电路中的量子门序列
5. **错误处理改进**:
- 使用 `raise ValueError(f"Invalid gate type: {gate_type}")` 替代 `exit(0)`
- 更友好的错误提示
### 使用说明:
1. **运行命令**:
```bash
python vqe_qsim.py --noise True \
--two_qubit_depolarizing_p 0.01 \
--single_qubit_depolarizing_p 0.001 \
--bit_flip_p 0.01
```
2. **输出监控**:
- 每个种子训练开始时会打印:"Starting training with seed X for circuit Y"
- 每次优化迭代会打印:"Seed X, Iter Y: Energy = Z.abcdef"
3. **内存管理**:
- 对于12量子比特系统,哈密顿量矩阵大小为 4096×4096 (约67MB)
- 如需更大系统,可考虑使用稀疏矩阵优化