高斯玻色采样 (Gaussian Boson Sampling, GBS)¶

难度：⭐⭐⭐（中级）¶

本教程将介绍高斯玻色采样（GBS）算法的原理、应用以及使用 DeepQuantum 的实现方法。

目录¶

1. GBS 算法原理
2. 玻色采样 vs 高斯玻色采样
3. 量子干涉和采样
4. 应用场景
5. DeepQuantum 实现
6. 示例：图聚类

---

1. GBS 算法原理¶

1.1 什么是高斯玻色采样？¶

高斯玻色采样（Gaussian Boson Sampling, GBS）是玻色采样的一种变体，它是展示量子计算优越性的重要方案之一。

核心思想：

• 将多个高斯压缩态输入到线性光学网络
• 光子在网络中发生量子干涉
• 在输出端进行光子数测量
• 获得输出光子数的概率分布

1.2 高斯态¶

高斯态指的是 Wigner 函数为高斯分布的量子态，包括：

• 相干态（Coherent state）
• 压缩态（Squeezed state）
• 真空态（Vacuum state）

GBS 使用压缩态作为输入，而不是离散的 Fock 态。

1.3 数学框架¶

Hafnian 函数¶

对于 $2m \times 2m$ 对称矩阵 $A$，Hafnian 定义为：

$$ \text{Haf}(A) = \sum_{\sigma \in \text{PMP}(2m)} \prod_{k=1}^{m} a_{\sigma(2k-1),\sigma(2k)} $$

其中 PMP 表示所有完美匹配排列的集合。

示例（n=4）： $$ \text{PMP}(4) = \{(0,1)(2,3), (0,2)(1,3), (0,3)(1,2)\} $$

$$ \text{Haf}(B) = B_{0,1}B_{2,3} + B_{0,2}B_{1,3} + B_{0,3}B_{1,2} $$

输出概率计算¶

对于粒子数分辨探测器，输出 Fock 态 $S = (s_1, s_2, ..., s_m)$ 的概率为：

$$ Pr(S) = \frac{\text{Haf}(A_S)}{\sqrt{\det(\Sigma)} \prod_{i=1}^m s_i!} $$

其中：

• $A_S$ 是根据输出模式构造的子矩阵
• $\Sigma$ 是协方差矩阵
• 对于纯态，$A = B \oplus B^*$，可以加速计算

图论意义¶

在图论中，Hafnian 计算了图 $G$ 的邻接矩阵 $A$ 对应的完美匹配数。

重要性质：

• 当图为二分图时，Hafnian = Permanent
• 因此 Hafnian 计算也是 #P 难问题

1.4 探测器类型¶

1. 粒子数分辨探测器（PNRD）¶

• 能够精确探测光子数
• 使用 Hafnian 函数计算概率

2. 阈值探测器（Threshold Detector）¶

• 只能判断是否有光子（0 或 ≥1）
• 使用 Torontonian 函数计算概率

$$ Pr(S) = \frac{\text{Tor}(O_S)}{\sqrt{\det(\Sigma)}} $$

其中 $O_S = I - (\Sigma^{-1})_S$

---

2. 玻色采样 vs 高斯玻色采样¶

2.1 玻色采样（Boson Sampling）¶

提出者： Aaronson & Arkhipov (2011)

特点：

• 输入：多个单光子 Fock 态 $|1,1,0,0,...\rangle$
• 通过线性光学网络干涉
• 输出概率与 Permanent 相关
• 理论证明：精确模拟是 #P 难问题

数学表达： $$ Pr(S) = \frac{|\text{Perm}(U_{st})|^2}{m_1! m_2! \cdots m_N! n_1! n_2! \cdots n_N!} $$

其中 $U_{st}$ 是酉矩阵 $U$ 的子矩阵。

2.2 高斯玻色采样（GBS）¶

特点：

• 输入：高斯压缩态
• 通过线性光学网络干涉
• 输出概率与 Hafnian 相关
• 实验上更容易实现

优势：

1. 实验可行性： 压缩态比单光子态更容易制备
2. 光子数： 可以产生更多光子
3. 应用广泛： 直接应用于图问题

2.3 对比总结¶

特性	玻色采样	高斯玻色采样
输入态	单光子 Fock 态	高斯压缩态
数学函数	Permanent	Hafnian
探测器	粒子数分辨	粒子数分辨/阈值
实验难度	较高	较低
应用	量子优越性	图问题、机器学习
光子数	有限	可扩展

2.4 为什么选择 GBS？¶

1. 更强的实用性： 不仅能展示量子优越性，还能解决实际问题
2. 图问题对应： 当所有 $s_i \in \{0,1\}$ 时，样本对应图的子图
3. 稠密子图偏好： 密度大的子图采样概率高
4. 应用领域： 聚类、图同构、机器学习等

---

3. 量子干涉和采样¶

3.1 多光子干涉¶

Hong-Ou-Mandel 效应：

• 两个全同光子进入分束器
• 由于量子干涉，光子总是成对输出
• 这是玻色采样的基础

3.2 线性光学网络¶

基本元件：

1. 相移器

   • 作用：改变相位
   • 矩阵：$\begin{pmatrix} e^{i\theta} & 0 \\ 0 & 1 \end{pmatrix}$

2. 分束器

   • 作用：混合模式
   • 矩阵：$\begin{pmatrix} \cos\theta & -e^{-i\phi}\sin\theta \\ e^{i\phi}\sin\theta & \cos\theta \end{pmatrix}$

任意酉矩阵分解： 根据 Reck 定理，任意 $N \times N$ 酉矩阵都可以用这些基本元件实现。

3.3 采样过程¶

步骤：

1. 制备输入态： $m$ 个压缩态
2. 线性光学变换： 通过酉矩阵 $U$
3. 光子数测量： 探测输出模式的光子数
4. 重复采样： 获得概率分布

示意图：

输入态 → [线性光学网络] → 输出态 → 测量
   ↓                              ↓
压缩态                      光子数分布

3.4 量子优势¶

为什么难以经典模拟？

1. Hafnian 计算： #P 完全问题
2. 指数级复杂度： $O(n 2^n)$ 对于 $n \times n$ 矩阵
3. 近似模拟也困难： 即使多项式时间近似也被证明困难

量子优势展示：

• 2019年：中国 "九章" 光量子计算机
• 76个光子，200秒完成采样
• 经典超算需要2.5亿年

---

4. 应用场景¶

GBS 不仅用于展示量子优越性，还有实际应用价值。

4.1 图同构¶

问题： 判断两个图是否同构（结构相同）

GBS 方法：

1. 将两个图的邻接矩阵编码到 GBS
2. 比较采样结果的分布
3. 同构的图产生相同的分布

优势：

• 利用量子干涉的敏感性
• 可以区分经典算法难以区分的图

4.2 图聚类¶

问题： 将图中的节点分成若干组

GBS 方法：

1. 构造数据点的邻接矩阵
2. 进行 GBS 采样
3. 寻找高密度子图（对应聚类）
4. 迭代直到所有节点被分类

原理：

• 稠密子图对应高概率样本
• 数据点相似 → 边密度高 → 容易成团

详细示例见第6节。

4.3 机器学习¶

应用方向：

1. 特征映射

   • GBS 样本作为核函数
   • 用于支持向量机（SVM）

2. 图神经网络

   • 利用 GBS 加速图特征提取

3. 数据生成

   • 生成符合特定分布的样本

4.4 其他应用¶

1. 稠密子图发现

   • 寻找图中连接最紧密的子图
   • 应用于社交网络分析

2. 最大团问题

   • 寻找完全连接的子图
   • NP 难问题

3. 分子相似性

   • 比较分子图结构
   • 药物发现应用

4.5 应用对比¶

应用	GBS 优势	经典方法
图聚类	自动发现稠密结构	需要预设簇数
图同构	量子指纹识别	复杂度较高
特征映射	高维核函数	计算昂贵
子图发现	偏好稠密子图	需要遍历

---

5. DeepQuantum 实现¶

5.1 环境准备¶

首先导入必要的库：

In [ ]:

import numpy as np
import deepquantum as dq
import deepquantum.photonic as dqp
import matplotlib.pyplot as plt
import torch

import numpy as np import deepquantum as dq import deepquantum.photonic as dqp import matplotlib.pyplot as plt import torch

5.2 基本 GBS 实现¶

创建一个 6 模式的高斯玻色采样设备：

In [ ]:

# 设置参数
nmode = 6  # 模式数
squeezing = [1.0] * nmode  # 压缩参数
unitary = np.eye(nmode)  # 酉矩阵（这里使用单位矩阵）

# 创建 GBS 设备
gbs = dq.photonic.GaussianBosonSampling(
    nmode=nmode,
    squeezing=squeezing,
    unitary=unitary
)

# 执行 GBS
gbs()

# 绘制线路
gbs.draw()

# 设置参数 nmode = 6 # 模式数 squeezing = [1.0] * nmode # 压缩参数 unitary = np.eye(nmode) # 酉矩阵（这里使用单位矩阵） # 创建 GBS 设备 gbs = dq.photonic.GaussianBosonSampling( nmode=nmode, squeezing=squeezing, unitary=unitary ) # 执行 GBS gbs() # 绘制线路 gbs.draw()

5.3 使用粒子数分辨探测器采样¶

In [ ]:

# 设置为粒子数分辨探测器
gbs.detector = 'pnrd'

# 进行采样（使用 MCMC 方法加速）
result_pnrd = gbs.measure(shots=1024, mcmc=True)

print("粒子数分辨探测器采样结果：")
print(result_pnrd)

# 统计最常见的输出
print("\n最常见的 5 个输出：")
for i, (state, count) in enumerate(list(result_pnrd.items())[:5]):
    print(f"{i+1}. {state}: {count} 次")

# 设置为粒子数分辨探测器 gbs.detector = 'pnrd' # 进行采样（使用 MCMC 方法加速） result_pnrd = gbs.measure(shots=1024, mcmc=True) print("粒子数分辨探测器采样结果：") print(result_pnrd) # 统计最常见的输出 print("\n最常见的 5 个输出：") for i, (state, count) in enumerate(list(result_pnrd.items())[:5]): print(f"{i+1}. {state}: {count} 次")

5.4 使用阈值探测器采样¶

In [ ]:

# 设置为阈值探测器
gbs.detector = 'threshold'

# 进行采样
result_threshold = gbs.measure(shots=1024, mcmc=True)

print("阈值探测器采样结果：")
print(result_threshold)

# 统计
print("\n输出统计：")
total_clicks = sum(1 for state in result_threshold.keys() if sum(int(s) for s in state) > 0)
print(f"有探测到的样本数: {total_clicks}")
print(f"无探测到的样本数: {result_threshold.get('|000000>', 0)}")

# 设置为阈值探测器 gbs.detector = 'threshold' # 进行采样 result_threshold = gbs.measure(shots=1024, mcmc=True) print("阈值探测器采样结果：") print(result_threshold) # 统计 print("\n输出统计：") total_clicks = sum(1 for state in result_threshold.keys() if sum(int(s) for s in state) > 0) print(f"有探测到的样本数: {total_clicks}") print(f"无探测到的样本数: {result_threshold.get('|000000>', 0)}")

5.5 自定义酉矩阵¶

我们可以创建更复杂的线性光学网络：

In [ ]:

# 生成随机酉矩阵
def random_unitary(n):
    """生成 n×n 随机酉矩阵"""
    # 使用 Haar 测度
    Z = np.random.randn(n, n) + 1j * np.random.randn(n, n)
    Q, R = np.linalg.qr(Z)
    return Q

# 创建随机酉矩阵
nmode = 8
U = random_unitary(nmode)

# 设置不同的压缩参数
squeezing = [1.0, 0.8, 1.2, 0.9, 1.1, 0.7, 1.0, 0.8]

# 创建 GBS
gbs_custom = dq.photonic.GaussianBosonSampling(
    nmode=nmode,
    squeezing=squeezing,
    unitary=U
)

# 采样
gbs_custom()
gbs_custom.detector = 'pnrd'
result_custom = gbs_custom.measure(shots=500, mcmc=True)

print(f"\n{nmode} 模式 GBS 采样结果（前10个）：")
for i, (state, count) in enumerate(list(result_custom.items())[:10]):
    print(f"{i+1}. {state}: {count} 次")

# 生成随机酉矩阵 def random_unitary(n): """生成 n×n 随机酉矩阵""" # 使用 Haar 测度 Z = np.random.randn(n, n) + 1j * np.random.randn(n, n) Q, R = np.linalg.qr(Z) return Q # 创建随机酉矩阵 nmode = 8 U = random_unitary(nmode) # 设置不同的压缩参数 squeezing = [1.0, 0.8, 1.2, 0.9, 1.1, 0.7, 1.0, 0.8] # 创建 GBS gbs_custom = dq.photonic.GaussianBosonSampling( nmode=nmode, squeezing=squeezing, unitary=U ) # 采样 gbs_custom() gbs_custom.detector = 'pnrd' result_custom = gbs_custom.measure(shots=500, mcmc=True) print(f"\n{nmode} 模式 GBS 采样结果（前10个）：") for i, (state, count) in enumerate(list(result_custom.items())[:10]): print(f"{i+1}. {state}: {count} 次")

5.6 理论概率验证¶

对于小规模系统，我们可以验证理论概率：

In [ ]:

from itertools import combinations

def hafnian(A):
    """计算对称矩阵的 Hafnian（简化版，仅用于小矩阵）"""
    n = len(A)
    if n % 2 != 0:
        return 0
    if n == 0:
        return 1
    if n == 2:
        return A[0, 1]
    
    # 递归计算（仅用于演示）
    total = 0
    for j in range(1, n):
        # 构造子矩阵
        idx = [i for i in range(n) if i not in [0, j]]
        A_sub = A[np.ix_(idx, idx)]
        total += A[0, j] * hafnian(A_sub)
    return total

# 小规模示例
print("Hafnian 计算示例：")
A = np.array([[0, 1, 1, 1],
              [1, 0, 1, 1],
              [1, 1, 0, 1],
              [1, 1, 1, 0]])

haf_value = hafnian(A)
print(f"Haf(A) = {haf_value}")
print("\n这表示该图有 3 个完美匹配")

from itertools import combinations def hafnian(A): """计算对称矩阵的 Hafnian（简化版，仅用于小矩阵）""" n = len(A) if n % 2 != 0: return 0 if n == 0: return 1 if n == 2: return A[0, 1] # 递归计算（仅用于演示） total = 0 for j in range(1, n): # 构造子矩阵 idx = [i for i in range(n) if i not in [0, j]] A_sub = A[np.ix_(idx, idx)] total += A[0, j] * hafnian(A_sub) return total # 小规模示例 print("Hafnian 计算示例：") A = np.array([[0, 1, 1, 1], [1, 0, 1, 1], [1, 1, 0, 1], [1, 1, 1, 0]]) haf_value = hafnian(A) print(f"Haf(A) = {haf_value}") print("\n这表示该图有 3 个完美匹配")

---

6. 示例：图聚类¶

本节将展示如何使用 GBS 进行图聚类。

6.1 问题定义¶

聚类目标：

• 将数据点分成若干组
• 同组内数据点相似
• 不同组间数据点不相似

GBS 方法优势：

• 无需预设聚类数量
• 自动发现稠密结构
• 利用量子干涉增强

6.2 数据准备¶

In [ ]:

# 使用西瓜数据集（经典机器学习示例）
# 每个数据点包含两个特征：密度和含糖率

data_ws = np.array([
    [0.697, 0.774, 0.634, 0.608, 0.556, 0.393, 0.451, 0.427, 0.666, 0.243, 
     0.245, 0.343, 0.639, 0.657, 0.725, 0.593, 0.6223, 0.75],
    [0.46, 0.376, 0.264, 0.318, 0.215, 0.237, 0.149, 0.211, 0.091, 0.267, 
     0.057, 0.099, 0.161, 0.198, 0.445, 0.082, 0.062, 0.405]
])

print(f"数据点数量: {data_ws.shape[1]}")
print(f"特征维度: {data_ws.shape[0]}")

# 可视化数据
plt.figure(figsize=(10, 6))
plt.scatter(data_ws[0], data_ws[1], c='blue', alpha=0.6, s=100)
plt.xlabel('密度', fontsize=12)
plt.ylabel('含糖率', fontsize=12)
plt.title('西瓜数据集', fontsize=14)
plt.grid(True, alpha=0.3)
plt.show()

# 使用西瓜数据集（经典机器学习示例） # 每个数据点包含两个特征：密度和含糖率 data_ws = np.array([ [0.697, 0.774, 0.634, 0.608, 0.556, 0.393, 0.451, 0.427, 0.666, 0.243, 0.245, 0.343, 0.639, 0.657, 0.725, 0.593, 0.6223, 0.75], [0.46, 0.376, 0.264, 0.318, 0.215, 0.237, 0.149, 0.211, 0.091, 0.267, 0.057, 0.099, 0.161, 0.198, 0.445, 0.082, 0.062, 0.405] ]) print(f"数据点数量: {data_ws.shape[1]}") print(f"特征维度: {data_ws.shape[0]}") # 可视化数据 plt.figure(figsize=(10, 6)) plt.scatter(data_ws[0], data_ws[1], c='blue', alpha=0.6, s=100) plt.xlabel('密度', fontsize=12) plt.ylabel('含糖率', fontsize=12) plt.title('西瓜数据集', fontsize=14) plt.grid(True, alpha=0.3) plt.show()

6.3 构建邻接矩阵¶

将数据点映射为图：

• 如果两个点距离 < 阈值，则存在边
• 形成图的邻接矩阵

In [ ]:

import itertools
import networkx as nx

def euclidean_distance(p1, p2):
    """计算欧氏距离"""
    return np.sqrt(np.sum((p1 - p2)**2))

def construct_adjacency_matrix(data, threshold, distance_func):
    """构造邻接矩阵"""
    n_points = data.shape[1]
    adj_mat = np.zeros((n_points, n_points))
    
    # 遍历所有点对
    for i, j in itertools.combinations(range(n_points), 2):
        dist = distance_func(data[:, i], data[:, j])
        if dist <= threshold:
            adj_mat[i, j] = 1
            adj_mat[j, i] = 1
    
    return adj_mat

# 设置距离阈值
distance_threshold = 0.2

# 构造邻接矩阵
adj_matrix = construct_adjacency_matrix(data_ws, distance_threshold, euclidean_distance)

print("邻接矩阵：")
print(adj_matrix)

# 可视化图
G = nx.from_numpy_array(adj_matrix)
pos = {i: data_ws[:, i] for i in range(len(data_ws[0]))}

plt.figure(figsize=(10, 6))
nx.draw(G, pos, with_labels=True, node_color='lightblue', 
        node_size=500, font_size=10, font_weight='bold',
        edge_color='gray', width=2)
plt.title('数据点构成的图', fontsize=14)
plt.show()

import itertools import networkx as nx def euclidean_distance(p1, p2): """计算欧氏距离""" return np.sqrt(np.sum((p1 - p2)**2)) def construct_adjacency_matrix(data, threshold, distance_func): """构造邻接矩阵""" n_points = data.shape[1] adj_mat = np.zeros((n_points, n_points)) # 遍历所有点对 for i, j in itertools.combinations(range(n_points), 2): dist = distance_func(data[:, i], data[:, j]) if dist <= threshold: adj_mat[i, j] = 1 adj_mat[j, i] = 1 return adj_mat # 设置距离阈值 distance_threshold = 0.2 # 构造邻接矩阵 adj_matrix = construct_adjacency_matrix(data_ws, distance_threshold, euclidean_distance) print("邻接矩阵：") print(adj_matrix) # 可视化图 G = nx.from_numpy_array(adj_matrix) pos = {i: data_ws[:, i] for i in range(len(data_ws[0]))} plt.figure(figsize=(10, 6)) nx.draw(G, pos, with_labels=True, node_color='lightblue', node_size=500, font_size=10, font_weight='bold', edge_color='gray', width=2) plt.title('数据点构成的图', fontsize=14) plt.show()

6.4 GBS 聚类算法¶

算法流程：

1. 构造邻接矩阵
2. 进行 GBS 采样
3. 寻找高密度子图（对应聚类）
4. 移除已找到的聚类
5. 重复直到完成聚类

第一轮聚类：

In [ ]:

# 创建 GBS 设备
gbs_cluster = dqp.GBS_Graph(adj_mat=adj_matrix, cutoff=2)

# 执行 GBS
state = gbs_cluster()

# 采样
print("开始 GBS 采样...（这可能需要一些时间）")
samples = gbs_cluster.measure(shots=10000, mcmc=True)

print(f"\n采样完成！获得 {len(samples)} 个不同的结果")
print("\n前10个采样结果：")
for i, sample in enumerate(list(samples.keys())[:10]):
    print(f"{i+1}. {sample}: {samples[sample]} 次")

# 创建 GBS 设备 gbs_cluster = dqp.GBS_Graph(adj_mat=adj_matrix, cutoff=2) # 执行 GBS state = gbs_cluster() # 采样 print("开始 GBS 采样...（这可能需要一些时间）") samples = gbs_cluster.measure(shots=10000, mcmc=True) print(f"\n采样完成！获得 {len(samples)} 个不同的结果") print("\n前10个采样结果：") for i, sample in enumerate(list(samples.keys())[:10]): print(f"{i+1}. {sample}: {samples[sample]} 次")

6.5 寻找稠密子图¶

In [ ]:

# 计算子图密度
graph_density = gbs_cluster.graph_density(G, samples)

# 设置参数
density_threshold = 0.8  # 图密度阈值
max_nodes = 12  # 从较大的子图开始

# 寻找满足条件的子图
found = False
target_nodes = None

while not found and max_nodes >= 4:
    for sample in samples:
        # 计算采样结果的节点数
        n_nodes = sum(int(s) for s in sample)
        
        if n_nodes == max_nodes:
            density = graph_density[sample][1]
            if density > density_threshold:
                print(f"\n找到聚类！")
                print(f"采样结果: {sample}")
                print(f"子图节点数: {n_nodes}")
                print(f"子图密度: {density:.4f}")
                
                # 获取节点索引
                target_nodes = torch.nonzero(sample.state).mT
                print(f"节点索引: {target_nodes}")
                
                found = True
                break
    
    if not found:
        max_nodes -= 2  # 降低节点数要求

if not found:
    print("未找到满足条件的子图，尝试降低密度阈值")

# 计算子图密度 graph_density = gbs_cluster.graph_density(G, samples) # 设置参数 density_threshold = 0.8 # 图密度阈值 max_nodes = 12 # 从较大的子图开始 # 寻找满足条件的子图 found = False target_nodes = None while not found and max_nodes >= 4: for sample in samples: # 计算采样结果的节点数 n_nodes = sum(int(s) for s in sample) if n_nodes == max_nodes: density = graph_density[sample][1] if density > density_threshold: print(f"\n找到聚类！") print(f"采样结果: {sample}") print(f"子图节点数: {n_nodes}") print(f"子图密度: {density:.4f}") # 获取节点索引 target_nodes = torch.nonzero(sample.state).mT print(f"节点索引: {target_nodes}") found = True break if not found: max_nodes -= 2 # 降低节点数要求 if not found: print("未找到满足条件的子图，尝试降低密度阈值")

6.6 可视化第一个聚类¶

In [ ]:

if target_nodes is not None:
    # 标记聚类节点
    node_colors = ['lightcoral' if i in target_nodes[0] else 'lightblue' 
                   for i in range(len(data_ws[0]))]
    
    plt.figure(figsize=(10, 6))
    nx.draw(G, pos, with_labels=True, node_color=node_colors, 
            node_size=500, font_size=10, font_weight='bold',
            edge_color='gray', width=2)
    plt.title(f'第一个聚类（{len(target_nodes[0])} 个节点）', fontsize=14)
    plt.show()
    
    print(f"\n第一个聚类包含节点: {target_nodes[0].tolist()}")
else:
    print("未找到聚类")

if target_nodes is not None: # 标记聚类节点 node_colors = ['lightcoral' if i in target_nodes[0] else 'lightblue' for i in range(len(data_ws[0]))] plt.figure(figsize=(10, 6)) nx.draw(G, pos, with_labels=True, node_color=node_colors, node_size=500, font_size=10, font_weight='bold', edge_color='gray', width=2) plt.title(f'第一个聚类（{len(target_nodes[0])} 个节点）', fontsize=14) plt.show() print(f"\n第一个聚类包含节点: {target_nodes[0].tolist()}") else: print("未找到聚类")

6.7 迭代聚类¶

移除第一个聚类，继续寻找其他聚类：

In [ ]:

# 移除已找到的聚类节点
if target_nodes is not None:
    # 创建新数据集（移除已聚类的点）
    remaining_indices = [i for i in range(len(data_ws[0])) 
                        if i not in target_nodes[0]]
    data_remaining = data_ws[:, remaining_indices]
    
    print(f"\n剩余数据点数: {len(remaining_indices)}")
    
    # 构造新的邻接矩阵
    adj_matrix_2 = construct_adjacency_matrix(
        data_remaining, distance_threshold, euclidean_distance
    )
    
    # 可视化剩余数据
    G2 = nx.from_numpy_array(adj_matrix_2)
    pos2 = {i: data_remaining[:, i] for i in range(len(data_remaining[0]))}
    
    plt.figure(figsize=(10, 6))
    nx.draw(G2, pos2, with_labels=True, node_color='lightgreen', 
            node_size=500, font_size=10, font_weight='bold',
            edge_color='gray', width=2)
    plt.title('剩余数据点构成的图', fontsize=14)
    plt.show()
    
    # 第二轮 GBS
    print("\n开始第二轮 GBS 采样...")
    gbs_cluster2 = dqp.GBS_Graph(adj_mat=adj_matrix_2, cutoff=2)
    state2 = gbs_cluster2()
    samples2 = gbs_cluster2.measure(shots=10000, mcmc=True)
    
    # 寻找第二个聚类
    graph_density2 = gbs_cluster2.graph_density(G2, samples2)
    density_threshold2 = 0.6
    max_nodes2 = 10
    
    found2 = False
    target_nodes2 = None
    
    while not found2 and max_nodes2 >= 4:
        for sample in samples2:
            n_nodes = sum(int(s) for s in sample)
            if n_nodes == max_nodes2:
                density = graph_density2[sample][1]
                if density > density_threshold2:
                    print(f"\n找到第二个聚类！")
                    print(f"采样结果: {sample}")
                    print(f"子图节点数: {n_nodes}")
                    print(f"子图密度: {density:.4f}")
                    
                    target_nodes2 = torch.nonzero(sample.state).mT
                    print(f"节点索引（在剩余数据中）: {target_nodes2[0].tolist()}")
                    
                    # 映射回原始索引
                    original_indices = [remaining_indices[i] for i in target_nodes2[0].tolist()]
                    print(f"原始节点索引: {original_indices}")
                    
                    found2 = True
                    break
        if not found2:
            max_nodes2 -= 2
    
    if not found2:
        print("未找到第二个聚类")
else:
    print("跳过迭代聚类")

# 移除已找到的聚类节点 if target_nodes is not None: # 创建新数据集（移除已聚类的点） remaining_indices = [i for i in range(len(data_ws[0])) if i not in target_nodes[0]] data_remaining = data_ws[:, remaining_indices] print(f"\n剩余数据点数: {len(remaining_indices)}") # 构造新的邻接矩阵 adj_matrix_2 = construct_adjacency_matrix( data_remaining, distance_threshold, euclidean_distance ) # 可视化剩余数据 G2 = nx.from_numpy_array(adj_matrix_2) pos2 = {i: data_remaining[:, i] for i in range(len(data_remaining[0]))} plt.figure(figsize=(10, 6)) nx.draw(G2, pos2, with_labels=True, node_color='lightgreen', node_size=500, font_size=10, font_weight='bold', edge_color='gray', width=2) plt.title('剩余数据点构成的图', fontsize=14) plt.show() # 第二轮 GBS print("\n开始第二轮 GBS 采样...") gbs_cluster2 = dqp.GBS_Graph(adj_mat=adj_matrix_2, cutoff=2) state2 = gbs_cluster2() samples2 = gbs_cluster2.measure(shots=10000, mcmc=True) # 寻找第二个聚类 graph_density2 = gbs_cluster2.graph_density(G2, samples2) density_threshold2 = 0.6 max_nodes2 = 10 found2 = False target_nodes2 = None while not found2 and max_nodes2 >= 4: for sample in samples2: n_nodes = sum(int(s) for s in sample) if n_nodes == max_nodes2: density = graph_density2[sample][1] if density > density_threshold2: print(f"\n找到第二个聚类！") print(f"采样结果: {sample}") print(f"子图节点数: {n_nodes}") print(f"子图密度: {density:.4f}") target_nodes2 = torch.nonzero(sample.state).mT print(f"节点索引（在剩余数据中）: {target_nodes2[0].tolist()}") # 映射回原始索引 original_indices = [remaining_indices[i] for i in target_nodes2[0].tolist()] print(f"原始节点索引: {original_indices}") found2 = True break if not found2: max_nodes2 -= 2 if not found2: print("未找到第二个聚类") else: print("跳过迭代聚类")

6.8 最终聚类结果¶

In [ ]:

# 可视化最终聚类结果
if target_nodes is not None:
    # 创建颜色映射
    colors = ['lightblue'] * len(data_ws[0])
    
    # 第一个聚类（红色）
    for i in target_nodes[0]:
        colors[i] = 'lightcoral'
    
    # 第二个聚类（绿色）
    if target_nodes2 is not None:
        original_indices = [remaining_indices[i] for i in target_nodes2[0].tolist()]
        for i in original_indices:
            colors[i] = 'lightgreen'
    
    # 绘制最终结果
    plt.figure(figsize=(12, 7))
    nx.draw(G, pos, with_labels=True, node_color=colors, 
            node_size=600, font_size=11, font_weight='bold',
            edge_color='gray', width=2)
    plt.title('GBS 聚类最终结果', fontsize=16)
    
    # 添加图例
    from matplotlib.patches import Patch
    legend_elements = [
        Patch(facecolor='lightcoral', label='聚类 1'),
        Patch(facecolor='lightgreen', label='聚类 2')
    ]
    if target_nodes2 is not None:
        legend_elements.append(
            Patch(facecolor='lightblue', label='聚类 3（剩余）')
        )
    plt.legend(handles=legend_elements, loc='upper right')
    
    plt.tight_layout()
    plt.show()
    
    print("\n聚类完成！")
    print(f"聚类 1: {len(target_nodes[0])} 个数据点")
    if target_nodes2 is not None:
        print(f"聚类 2: {len(target_nodes2[0])} 个数据点")
        remaining = len(data_ws[0]) - len(target_nodes[0]) - len(target_nodes2[0])
        print(f"聚类 3: {remaining} 个数据点")
else:
    print("未能完成聚类")

# 可视化最终聚类结果 if target_nodes is not None: # 创建颜色映射 colors = ['lightblue'] * len(data_ws[0]) # 第一个聚类（红色） for i in target_nodes[0]: colors[i] = 'lightcoral' # 第二个聚类（绿色） if target_nodes2 is not None: original_indices = [remaining_indices[i] for i in target_nodes2[0].tolist()] for i in original_indices: colors[i] = 'lightgreen' # 绘制最终结果 plt.figure(figsize=(12, 7)) nx.draw(G, pos, with_labels=True, node_color=colors, node_size=600, font_size=11, font_weight='bold', edge_color='gray', width=2) plt.title('GBS 聚类最终结果', fontsize=16) # 添加图例 from matplotlib.patches import Patch legend_elements = [ Patch(facecolor='lightcoral', label='聚类 1'), Patch(facecolor='lightgreen', label='聚类 2') ] if target_nodes2 is not None: legend_elements.append( Patch(facecolor='lightblue', label='聚类 3（剩余）') ) plt.legend(handles=legend_elements, loc='upper right') plt.tight_layout() plt.show() print("\n聚类完成！") print(f"聚类 1: {len(target_nodes[0])} 个数据点") if target_nodes2 is not None: print(f"聚类 2: {len(target_nodes2[0])} 个数据点") remaining = len(data_ws[0]) - len(target_nodes[0]) - len(target_nodes2[0]) print(f"聚类 3: {remaining} 个数据点") else: print("未能完成聚类")

6.9 算法总结¶

GBS 聚类算法流程：

1. 数据点 → 邻接矩阵
   ↓
2. 邻接矩阵 → GBS 采样
   ↓
3. 采样结果 → 计算子图密度
   ↓
4. 高密度子图 → 聚类
   ↓
5. 移除聚类 → 重复步骤 2-4
   ↓
6. 完成聚类

关键参数：

• distance_threshold: 距离阈值（构造邻接矩阵）
• density_threshold: 密度阈值（识别聚类）
• cutoff: 截断参数（控制光子数）
• shots: 采样次数（影响结果质量）

优化建议：

1. 增加采样次数以提高准确度
2. 调整距离阈值以获得合适的图结构
3. 使用 MCMC 方法加速采样
4. 考虑使用阈值探测器以降低实验难度

---

总结¶

本教程介绍了高斯玻色采样（GBS）的原理、实现和应用。

核心要点¶

1. GBS 原理

   • 输入：高斯压缩态
   • 过程：线性光学网络干涉
   • 输出：光子数概率分布
   • 数学：Hafnian 函数（#P 难）

2. vs 玻色采样

   • 更容易实验实现
   • 更强的实用性
   • 更广泛的应用

3. 量子优势

   • 展示量子优越性
   • 加速特定问题求解
   • 解决实际应用问题

4. 应用场景

   • 图同构
   • 图聚类
   • 机器学习
   • 子图发现

5. DeepQuantum 实现

   • 简单易用的 API
   • 支持多种探测器
   • MCMC 加速采样
   • 图问题专用接口

进一步学习¶

• 尝试不同的数据集进行聚类
• 探索 GBS 在其他图问题中的应用
• 比较不同探测器类型的效果
• 研究 GBS 与经典算法的性能对比
• 了解最新的 GBS 实验进展

参考文献¶

1. Aaronson, S., & Arkhipov, A. (2011). The computational complexity of linear optics. Theory of Computing, 9, 143-252.

2. Hamilton, C. S., et al. (2017). Gaussian boson sampling. Physical Review Letters, 119(17), 170501.

3. Bromley, T. R., et al. (2020). Applications of near-term photonic quantum computers: software and algorithms. Quantum Science and Technology, 5(3), 034010.

4. Quesada, N., Arrazola, J. M., & Killoran, N. (2018). Gaussian boson sampling using threshold detectors. Physical Review A, 98(6), 062322.

5. Bonaldi, N., et al. (2023). Boost clustering with Gaussian Boson Sampling: a full quantum approach. arXiv:2307.13348.

6. Zhong, H. S., et al. (2020). Quantum computational advantage using photons. Science, 370(6523), 1460-1463.

---

教程完成！如有问题，请参考 DeepQuantum 文档或联系开发者。