计算量子核矩阵

我们来计算这些量子核。如前所述，主要思想是运用量子计算机来估计数据点 $x_i$ 和 $x_j$ 之间的相似性，这种相似性体现在它们对应量子特征态 $|\phi(x_i)\rangle$ 和 $|\phi(x_j)\rangle$ 的内积中。我们通过特征映射电路 $U_{\phi(x)}$ 来表示这些态，这些电路作用于初始的 $|0\rangle^{\otimes n}$ 态：$|\phi(x)\rangle = U_{\phi(x)}|0\rangle^{\otimes n}$。核矩阵的单个元素通常定义为该内积的平方模：

$$k(x_i, x_j) = |\langle \phi(x_j) | \phi(x_i) \rangle|^2$$

通过模拟态向量 $|\phi(x_i)\rangle$ 和 $|\phi(x_j)\rangle$ 直接计算此内积仅在量子比特数量非常少时可行，因为希尔伯特空间维度会指数级增长（$2^n$）。量子计算提供了高效估计此值的方法，即使维度很高时也是如此。

使用量子电路估计核矩阵元素

有几种常用的量子电路设计用于估计 $k(x_i, x_j)$。

1. 重叠法（逆电路法）

这可以说是最直接且常用的方法，尤其是在模拟中以及与变分算法结合时。主要思想是准备态 $|\phi(x_i)\rangle$，然后施加 $x_j$ 对应的特征映射的逆操作，表示为 $U_{\phi(x_j)}^\dagger$。

步骤如下：

1. 初始化：将 $n$ 量子比特系统初始化为全零态：$|0\rangle^{\otimes n}$。
2. 施加第一个特征映射：施加对应于第一个数据点 $x_i$ 的幺正电路 $U_{\phi(x_i)}$。态变为 $|\psi_1\rangle = U_{\phi(x_i)}|0\rangle^{\otimes n} = |\phi(x_i)\rangle$。
3. 施加第二个逆特征映射：施加对应于第二个数据点 $x_j$ 的逆（共轭转置）电路 $U_{\phi(x_j)}^\dagger$。最终态为 $|\psi_{final}\rangle = U_{\phi(x_j)}^\dagger U_{\phi(x_i)}|0\rangle^{\otimes n}$。
4. 测量：对所有 $n$ 个量子比特进行计算基测量。
5. 估计概率：重复步骤1-4多次（例如 $N_{\text{shots}}$ 次），并计数结果为全零态 $|0\rangle^{\otimes n}$ 的次数。估计概率为 $P_{\text{est}}(|0\rangle^{\otimes n}) = \frac{\text{计数}(|0\rangle^{\otimes n})}{N_{\text{shots}}}$。

为什么这样做可行？从最终态 $|\psi_{final}\rangle$ 测量到全零态 $|0\rangle^{\otimes n}$ 的概率由玻恩规则给出：

$$P(|0\rangle^{\otimes n}) = |\langle 0|^{\otimes n} | \psi_{final} \rangle|^2 = |\langle 0|^{\otimes n} U_{\phi(x_j)}^\dagger U_{\phi(x_i)} |0\rangle^{\otimes n}|^2$$

识别出 $\langle 0|^{\otimes n} U_{\phi(x_j)}^\dagger$ 是 $\langle \phi(x_j)|$ 的左矢，而 $U_{\phi(x_i)} |0\rangle^{\otimes n}$ 是 $|\phi(x_i)\rangle$ 的右矢，我们得到：

$$P(|0\rangle^{\otimes n}) = |\langle \phi(x_j) | \phi(x_i) \rangle|^2 = k(x_i, x_j)$$

因此，所需的核值就是运行这个组合电路后测量到全零态的概率。这种方法需要一个量子电路，其深度大致等于 $U_{\phi(x_i)}$ 和 $U_{\phi(x_j)}$ 深度之和。

以下是说明电路结构的图表：

使用重叠法估计 $k(x_i, x_j) = |\langle \phi(x_j) | \phi(x_i) \rangle|^2$ 的电路图。系统被初始化，施加特征映射 $U_{\phi(x_i)}$，接着施加逆特征映射 $U_{\phi(x_j)}^\dagger$，最后测量所有量子比特。 $|0...0\rangle$ 结果的概率估计核矩阵元素。

2. SWAP测试变体

另一种常用方法涉及一个辅助量子比特（辅助比特）。虽然有时被称为“SWAP测试”，但通常用于核估计的电路比用于态比较的完整SWAP测试略为简单。

1. 初始化：首先为 $|\phi(x_i)\rangle$ 准备 $n$ 个量子比特，为 $|\phi(x_j)\rangle$ 准备 $n$ 个量子比特，以及一个辅助量子比特。初始化为 $|0\rangle_{\text{anc}} |0\rangle^{\otimes n} |0\rangle^{\otimes n}$。
2. 准备态：对第二个寄存器施加 $U_{\phi(x_i)}$，对第三个寄存器施加 $U_{\phi(x_j)}$。态为 $|0\rangle_{\text{anc}} |\phi(x_i)\rangle |\phi(x_j)\rangle$。
3. 哈达玛门：对辅助比特施加一个哈达玛门：$\frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)_{\text{anc}} |\phi(x_i)\rangle |\phi(x_j)\rangle$。
4. 受控SWAP门：在第二个和第三个寄存器之间施加一个受辅助量子比特控制的SWAP操作。态变为 $\frac{1}{\sqrt{2}}(|0\rangle_{\text{anc}} |\phi(x_i)\rangle |\phi(x_j)\rangle + |1\rangle_{\text{anc}} |\phi(x_j)\rangle |\phi(x_i)\rangle)$。
5. 哈达玛门：对辅助比特再施加一个哈达玛门。
6. 测量辅助比特：在计算基下测量辅助比特。
7. 估计概率：重复步骤1-6共 $N_{\text{shots}}$ 次，并估计辅助比特测量为态 $|0\rangle$ 的概率 $P_{\text{est}}(|0\rangle_{\text{anc}})$。

在辅助比特上施加第二个哈达玛门后，$|0\rangle_{\text{anc}}$ 分量的系数是 $\frac{1}{2}(|\phi(x_i)\rangle |\phi(x_j)\rangle + |\phi(x_j)\rangle |\phi(x_i)\rangle)$。测量到 $|0\rangle$ 的概率是该分量在辅助比特 $|0\rangle$ 态上的投影的平方范数，其简化为：

$$P(|0\rangle_{\text{anc}}) = \frac{1}{2} (1 + |\langle \phi(x_j) | \phi(x_i) \rangle|^2)$$

因此，核值可以提取为：

$$k(x_i, x_j) = |\langle \phi(x_j) | \phi(x_i) \rangle|^2 = 2 \cdot P(|0\rangle_{\text{anc}}) - 1$$

这种方法需要 $2n+1$ 个量子比特，并涉及一个可能复杂的受控SWAP操作，这会显著增加电路深度和门数量，相比于重叠法，尤其是当SWAP操作需要分解为原生门时。然而，它提供了与内积平方直接相关的估计值。

实际应用中，重叠法（逆电路）通常更受青睐，因为它所需的量子比特数较少且电路可能更浅，使其更适合近期的量子设备和模拟器。

统计估计与采样噪声

特别要注意的是，两种方法都通过概率估计间接得到核值。因为量子测量本质上是概率性的，我们需要多次执行相关电路（$N_{\text{shots}}$ 次），并使用目标结果的频率（例如 $|0\rangle^{\otimes n}$ 或 $|0\rangle_{\text{anc}}$）来估计真实概率。

这会引入统计噪声或“采样”噪声。核矩阵元素估计 $k_{\text{est}}(x_i, x_j)$ 的准确性取决于 $N_{\text{shots}}$。估计的标准差通常按 $O(1/\sqrt{N_{\text{shots}}})$ 缩放。实现高精度需要大量的测量次数，这会增加计算成本。

构建完整核矩阵

经典基于核的算法，如支持向量机（SVM），需要整个格拉姆矩阵 $K$，其中元素 $K_{ij} = k(x_i, x_j)$ 表示训练集 $\{x_1, x_2, ..., x_m\}$ 中第 $i$ 个和第 $j$ 个数据点之间的核评估。

要使用量子方法构建这个 $m \times m$ 矩阵：

1. 遍历数据对：遍历所有唯一的数据点对 $(x_i, x_j)$，其中 $1 \le i \le j \le m$。由于 $k(x_i, x_j) = k(x_j, x_i)$（因为 $|\langle \phi(x_j) | \phi(x_i) \rangle|^2 = |\langle \phi(x_i) | \phi(x_j) \rangle|^2$），我们只需计算包含对角线的上（或下）三角部分。这相当于 $m(m+1)/2$ 个唯一的对。
2. 估计每个元素：对于每对 $(x_i, x_j)$，构建并执行所选的量子电路（例如，重叠法电路 $U_{\phi(x_j)}^\dagger U_{\phi(x_i)}$）共 $N_{\text{shots}}$ 次。
3. 计算 $K_{ij}$：根据测量统计数据计算核值 $K_{ij}$（例如，对于重叠法， $K_{ij} = P_{\text{est}}(|0\rangle^{\otimes n})$）。
4. 填充矩阵：将计算值存储到 $K_{ij}$ 和 $K_{ji}$ 中。

总成本涉及运行 $\approx m^2/2$ 个不同的量子电路，每个运行 $N_{\text{shots}}$ 次。这种数据集大小 $m$ 的二次缩放，结合每次估计的成本（$N_{\text{shots}}$），是量子核算法整体运行时间中的一个重要因素。

实现说明：模拟器对比硬件

核矩阵的计算方式在很大程度上取决于执行后端：

• 态向量模拟器：如果量子比特数 $n$ 足够小（通常小于30），经典模拟器可以存储完整的量子态向量。在这种情况下，你可以直接模拟电路 $U_{\phi(x_i)}$ 和 $U_{\phi(x_j)}$ 来得到态向量 $|\phi(x_i)\rangle$ 和 $|\phi(x_j)\rangle$，然后以经典方式计算它们的内积 $\langle \phi(x_j) | \phi(x_i) \rangle$。这完全避免了采样噪声，但受限于经典内存。这对于调试和小型测试来说非常有效。
• Qasm模拟器/真实硬件：对于更大的量子比特数或在实际量子处理器上运行时，你必须使用上述基于电路的估计方法（重叠法或SWAP测试）。
  • Qasm模拟器：它们通过采样结果来模拟量子计算的概率性质，模仿真实硬件但没有物理噪声。你需要指定 $N_{\text{shots}}$。
  • 量子硬件：在硬件上执行会引入物理噪声（退相干、门错误、读出错误），除了固有的采样噪声之外。估计概率 $P_{\text{est}}$ 将受到这种噪声的干扰。第7章（硬件考量和误差缓解）中讨论的技术对于从硬件获取有意义的结果变得必不可少。

以下伪代码概述了在需要测量次数的后端上使用重叠法计算单个核矩阵元素的过程：

# 类似于Python的伪代码
import numpy as np
# 假设存在 QuantumCircuit 构建器和 Backend 运行器
# （例如来自 Qiskit, PennyLane, Cirq 等）

def build_feature_map_circuit(data_point, num_qubits):
  """为给定数据点构建量子电路 U_phi(x)。"""
  qc = QuantumCircuit(num_qubits)
  # ... 根据数据点和所选的特征映射策略添加门 ...
  # 示例：ZZ特征映射、泡利特征映射或自定义电路
  # qc.rx(data_point[0], 0)
  # qc.ry(data_point[1], 1)
  # qc.cz(0, 1)
  # ...
  return qc

def estimate_kernel_entry_overlap(x_i, x_j, feature_map_builder, num_qubits, backend, num_shots):
  """使用重叠法估计 k(xi, xj)。"""

  # 1. 构建 U_phi(xi)
  circuit_i = feature_map_builder(x_i, num_qubits)

  # 2. 构建 U_phi(xj) 并获取其逆
  circuit_j = feature_map_builder(x_j, num_qubits)
  try:
    circuit_j_dagger = circuit_j.inverse()
  except Exception as e:
    print(f"Warning: 无法自动反转 xj 的电路。请确保构建器支持反转。错误：{e}")
    # 或者如果特定特征映射需要，手动处理反转
    raise e

  # 3. 组合电路：U_phi(xj)^dagger * U_phi(xi)
  # 注意：电路组合顺序可能因库而异（qc1.compose(qc2) 对比 qc2(qc1)）
  # 假设 qc1.compose(qc2) 先应用 qc2，然后应用 qc1
  measurement_circuit = circuit_j_dagger.compose(circuit_i)

  # 4. 添加计算基测量
  # 确保测量发生在所有幺正操作*之后*
  measurement_circuit.measure_all() # 或者如果需要，测量特定量子比特

  # 5. 在后端执行
  # 转译可能在此处隐式或显式发生
  job = backend.run(measurement_circuit, shots=num_shots)
  result = job.result()
  counts = result.get_counts(measurement_circuit)

  # 6. 计算概率 P(|0...0>)
  zero_string = '0' * num_qubits
  prob_zero = counts.get(zero_string, 0) / num_shots

  kernel_value = prob_zero
  return kernel_value

# --- 示例用法 ---
# 定义你的特征映射逻辑
# def my_feature_map(data_point, num_qubits): ...

# num_qubits = 4
# data_point_1 = np.array([0.1, 0.2, 0.3, 0.4])
# data_point_2 = np.array([0.5, 0.6, 0.7, 0.8])

# backend = Aer.get_backend('qasm_simulator') # Qiskit模拟器示例
# num_shots = 8192

# k_12 = estimate_kernel_entry_overlap(
#     data_point_1,
#     data_point_2,
#     my_feature_map,
#     num_qubits,
#     backend,
#     num_shots
# )

# 打印(f"估计的核值 k(x1, x2): {k_12}")

这个过程构成了应用量子核方法的核心计算过程。在研究QSVM等特定算法之前，了解这些元素是如何估计的、相关成本以及模拟与硬件执行之间的差异是十分基础的。