论文背景

问题: 给定一个势函数 $V:\mathbb{R}^d→\mathbb{R}$ 和一个 inverse temperature $\beta>0$, 目标是从 Gibbs 分布 $$\begin{equation} σ \propto e^{-βV} \end{equation}$$ 采样 (有时称为 classical Gibbs sampling).

看上去这个概率分布是直接给出的, 但是高维的时候直接采样会遇到高维积分等本质难题.

基础知识补充

随机过程基础知识(发现这些全都用不上)

基本定义

可测空间: $(S, \mathcal{S})$, 其中 $S$ 是集合, $\mathcal{S}$ 是 $S$ 上的 σ-代数.

测度: 测度是一个特殊的可测映射 $μ:(S, \mathcal{S})→([0,+∞], \mathcal{B}([0,+∞]))$, 满足:

1. $μ(∅) = 0$.
2. 可数可加性: 对任意可数个两两不交的集合 $A_1, A_2, ... ∈ \mathcal{S}$, 有 $$\begin{equation} μ\left(\bigcup_{i=1}^{∞} A_i\right) = \sum_{i=1}^{∞} μ(A_i). \end{equation}$$

测度空间: $(S, \mathcal{S}, μ)$, 其中 $μ$ 是定义在 $(S, \mathcal{S})$ 上的测度.

概率空间: $(\Omega, \mathcal{F}, \mathbb{P})$, 其中 $\mathbb{P}$ 是定义在 $(\Omega, \mathcal{F})$ 上的概率测度, 满足 $\mathbb{P}(\Omega) = 1$.

随机变量: 定义在概率空间 $(\Omega, \mathcal{F}, \mathbb{P})$ 上, 值域为可测空间 $(S, \mathcal{S})$ 的可测映射 $X:(\Omega, \mathcal{F})→(S, \mathcal{S})$.

随机变量(狭义): 定义在概率空间 $(\Omega, \mathcal{F}, \mathbb{P})$ 上, 值域为实数集 $\mathbb{R}$ 的可测映射 $X:(\Omega, \mathcal{F})→(\mathbb{R}, \mathcal{B}(\mathbb{R}))$.

概率分布: 随机变量 $X$ 的概率分布是测度 $μ$ 定义为概率测度的 push-forward: $$\begin{equation} μ(A) = \mathbb{P}(X^{-1}(A)), \quad ∀ A ∈ \mathcal{S}. \end{equation}$$

Markov process

$$\begin{equation} \mathbf{P}\left\{X_{t+1}=y\mid H_{t-1}\cap\{X_t=x\}\right\}=\mathbf{P}\left\{X_{t+1}=y\mid X_t=x\right\}=P(x,y). \end{equation}$$

因为“马尔可夫性”是一种关于“信息”的性质，而不是关于“路径值”的性质。 所以必须用 filtration（信息随时间增长的 σ-域流）才能精确定义。

布朗运动

布朗运动 = 从 0 开始、增量是独立的高斯随机变量，并且路径连续的随机过程.

路径虽然连续但是几乎处处不可微.

重要的性质是: $(\mathrm{d}W_t)^2 ≈ \mathrm{d}t$, 其中 $\mathrm{d}W_t$ 的意思是 $W_{t+\mathrm{d}t} - W_t$.

随机微积分 (Itô calculus)

最基础的构造:

$$\begin{equation} \int_0^t f(s)\mathrm{d}W_s := \lim_{n→∞} \sum_{i=0}^{n-1} f(t_i)(W_{t_{i+1}} - W_{t_i}), \end{equation}$$

这里必须在左端点取值, 这是 Itô 积分的定义. (与前向欧拉相关联)

SDE

SDE 给的是样本轨迹 / path dynamics.

SDE 描述了一个随机过程的演化. 其中 $b$ 称为 drift coefficient, $σ$ 称为 diffusion coefficient.

概率分布

一个随机变量的概率分布是样本集 $\Omega$ 上的概率测度 $\mathbb{P}$ 由随机变量 push-forward 到值域上的测度.

具体形式化:

设 $(\Omega, \mathcal{F}, \mathbb{P})$ 是概率空间, $X:\Omega→ S$ 是值域为可测空间 $(S, \mathcal{S})$ 的随机变量. 则 $X$ 的概率分布是测度 $\mu$ 定义为 $$\begin{equation} \mu(A) = \mathbb{P}(X^{-1}(A)), \quad ∀ A ∈ \mathcal{S}. \end{equation}$$

平稳分布 (stationary distribution)

对于一个随机过程 $X_t$ (通常是马尔可夫过程), 如果存在一个概率分布 $π$ 满足: 当 $X_0$ 服从 $π$ 时, 对任意的 $t≥0$, $X_t$ 仍然服从 $π$. 则称 $π$ 是该随机过程的平稳分布.

Fokker-Planck equation

Fokker–Planck 给的是分布 dynamics / law 的演化.

给定 SDE $$\begin{equation} \mathrm{d}X_t = b(X_t)\mathrm{d}t + σ(X_t)\mathrm{d}W_t, \end{equation}$$ Fokker-Planck 方程描述了 $X_t$ 的概率密度函数 $ρ(x,t)$ 随时间的演化: $$\begin{equation} \partial_t ρ(t,x) = L^* ρ(t,x), \end{equation}$$ 其中 $L^*$ 是生成元 $L$ 的伴随算子.

平稳分布是 Fokker-Planck 方程的稳态解, 即满足 $L^* ρ = 0$ 的概率密度函数 $ρ$.

Detailed balance condition

Langevin dynamics

Langevin dynamics 是一种马尔可夫过程, 它的平稳分布是目标分布 $σ$. 它的 SDE 形式为 $$\begin{equation} dX_t = -∇V(X_t)dt + \sqrt{2/β}dW_t, \end{equation}$$ 其中 $W_t$ 是标准布朗运动. 它的 Fokker-Planck 方程为 $$\begin{equation} \frac{∂ρ}{∂t} = ∇⋅(ρ∇V) + \frac{1}{β}Δρ. \end{equation}$$

解决总体思路

Gibbs 采样的分布其实是 overdamped Langevin equation $$\begin{equation} dX_t = -∇V(X_t)dt + \sqrt{2/β}dW_t \end{equation}$$ 的平稳分布.

下面的 Fokker-Planck equation 描述了概率密度函数 $ρ(x,t)$ 的演化: $$\begin{equation} \frac{∂ρ}{∂t} = \mathcal{L}(\rho) := ∇⋅(ρ∇V) + \frac{1}{β}Δρ. \end{equation}$$ 里面的 $\mathcal{L}$ 是半群生成元 (semigroup generator). (P.S. 这篇文章里面用的记号和大部分常用教材不一致, 通常这里是 $\mathcal{L}^*$) $\mathcal{L}$ 的标准 $L^2$ 伴随算子是 $$\begin{equation} \mathcal{L}^† = -∇V⋅∇ + \frac{1}{β}Δ. \end{equation}$$

一个重要的性质是 $\mathcal{L}^†$ 关于 $σ$-weighted inner product 自伴.

Gibbs 测度就是这个 dynamics 的平稳分布, 即满足 $\mathcal{L}(\sigma) = 0$ 的概率密度函数 $\sigma$.

论文里面定义了一个 Witten Laplacian $$\begin{equation} \mathcal{H} = -σ^{-1/2}∘ \mathcal{L} ∘ σ^{1/2} = -\frac{1}{β}Δ + \left(\frac{\beta \Vert∇V\Vert^2}{4} - \frac{1}{2} ΔV\right), \end{equation}$$ 它关于标准 $L^2$ 内积自伴.

(它是 metastability 研究中的标准工具)

在 detailed balance 情况下(等价于 $\mathcal{L}^†$ 关于 $σ$-inner product 是自伴的), $\mathcal{H}$ 在 $L^2(\mathrm{d}x)$ 上的 spectra 与 $-\mathcal{L}$ (或者 $-\mathcal{L}^†$) 在 $L^2(σ\mathrm{d}x)$ 上的 spectra 一致.

$\mathcal{H}$ 的 kernel 就是 encoded Gibbs state $|\sqrt{σ}⟩$, 满足 $⟨x|\sqrt{σ}⟩=\sqrt{\sigma(x)}$ (这是一个标准 $L^2$-inner product 下的归一化态).

这里说 encoded 是因为并不是直接得到 Gibbs distribution, 而是把它的平方根编码到一个量子态里面, 使得“测量=抽样”.

在理论模型里面, 我们有 $|\sqrt{σ}⟩∈L^2(\mathbb{R}^d, \mathrm{d}x)$, $|x⟩$ 其实不是 $L^2(\mathbb{R}^d, \mathrm{d}x)$ 里面的一个态, 它其实是 Dirac 泛函 (分布) $δ_x$, 并且这里的 $x$ 是一个参量, 表示位置, $⟨x|\sqrt{σ}⟩$ 也不是通常的内积, 结果不是一个值, 而是表示输入不同 $x$ 时的函数值, 从而 $⟨x|\sqrt{σ}⟩$ 表示的是位置表象下的波函数 $\sqrt{σ(x)}$.

按照论文的 assumption, 当 Langevin dynamics 是 ergodic 的时候, $|\sqrt{\sigma}⟩$ 是 $\mathcal{H}$ 的唯一基态, 它有 spectral gap (谱隙, 表示离它的本征值, 也就是 $0$, 最近的本征值和它的差距) $\text{Gap}(\mathcal{H}) = \text{Gap}(\mathcal{L}^†)$.

所以现在, 我们就把采样 Gibbs distribution 的问题, 转化成了制备 $\mathcal{H}$ 的基态 $|\sqrt{σ}⟩$ 的问题.

在 $\mathcal{H}$ 的形式化定义里面, 需要对 $V$ 求二阶导, 而 Langevin dynamics 只需要对 $V$ 求一阶导. 下面的分解顺便解决了这个问题: $$\begin{equation} \mathcal{H} = \sum_{j=1}^d L_j^† L_j, \quad L_j = -i\frac{1}{\sqrt{β}} ∂_{x_j} -i \frac{\sqrt{β}}{2} ∂_{x_j}V, ∀ j=1,2,...,d. \end{equation}$$

注意到每个 $L_j$ 都零化 $|\sqrt{σ}⟩$, 也就是说 $|\sqrt{\sigma}⟩$ 是每一项 $L_j^† L_j$ 的基态 (但此时未必是唯一的基态), 只要将每一项 $L_j^† L_j$ 的基态空间作交, 就得到 $|\sqrt{σ}⟩$. 把所有的 $L_i$ 放到一起, 定义 $$\begin{equation} \mathbb{L}:= [L_1^T, L_2^T, ..., L_d^T]^T, \end{equation}$$ 这样就有 $\mathcal{H} = \mathbb{L}^†\mathbb{L}$. (构造 $\mathbb{L}$ 只需要一阶导数)

进一步, 寻找 $\mathcal{H}$ 的基态 (encoded Gibbs state) $|\sqrt{\sigma}⟩$ 相当于准备 $\mathbb{L}$ 的唯一的对应于 singular value $0$ 的 right singular vector. $\mathbb{L}$ 的 singular value gap 是 $\text{Gap}(\mathbb{L}) = \sqrt{\text{Gap}(\mathcal{H})} = 1 / \sqrt{C_{\text{PI}}}$. 这是 operator level 的 quantum speedup 的来源(我觉得这个是一个需要理解的关键点, 似乎这个技术可以扩展). 这里 speedup 是源于 QSVT 算法的特性, 它的时间复杂度与 singular value gap 成反比, 所以如果 singular value gap 变大, 就可以加速制备基态的过程(详情可见 https://doi.org/10.1145/3313276.3316366. 的 3.4).

并且我觉得还有一个需要理解的问题是, 如果进行这样的分解(类似于取平方根, 关键在于让 spectral gap 取根号)就可以加速, 为什么不直接在传统方法(也就是直接解 Fokker-Planck 方程)上进行类似的分解呢? (感觉主要需要去分析两个复杂度和 $C_{\text{PI}}$ 有关的 Classical algorithm: ULA 和 Proximal (https://chewisinho.github.io))

ULA (Unadjusted Langevin Algorithm) 就是直接对 Langevin dynamics 进行 Euler-Maruyama 离散化, 我感觉就是模拟动力学过程, 直观上是梯度下降 + 高斯噪声. 这时 spectral gap 是和动力学过程本身有关的, 没有办法取根号.

Proximal 算法是一个 Gibbs 采样的增广变量技巧, $C_{\text{PI}}$ 出现在迭代的分析里面. (还没有仔细看)

回到算法, 算法的主要步骤是:

1. 准备一个 warm start 态 $|\phi⟩$ (与目标态有较大重叠的态), 这一点容易做到, 因为我们不可能对于目标分布没有一点点知识.

2. 用量子电路对离散化的矩阵 $\mathbb{L}$ 进行 block-encoding, 也就是 把 $\mathbb{L}$ 嵌入到一个更大的酉矩阵里面.

3. 使用 QSVT 去 filter out $|\phi⟩$ 中的非零 singular value 部分, 得到目标态 $|g⟩≈|\sqrt{σ}⟩$.

4. 在计算基上测量 $|g⟩$, 得到 Gibbs 采样结果.

下面具体要说明一下在量子计算框架(我们的框架是门模型(qubit)框架)下, 第四步是什么意思: 上述步骤得到的量子态 $|g⟩$ 假设是 $n$ 个 qubit 的, 它在计算基上的表示为 $$\begin{equation} |g⟩ ≈ \sum_{k\in\{0,1\}^n} g_k |k⟩, \end{equation}$$ 其中 $|k⟩$ 是 computational basis. 一次测量就是对每个 qubit 进行一次标准读出, 每次会得到一个 $n$ 位的二进制字符串 $k$. Born 规则告诉我们, 有分布 $$\begin{equation} P(\text{输出} k) = |g_k|^2, \end{equation}$$ 所以第四步采样的本质是在概率分布 $\{p_k = |g_k|^2\}$ 中抽样一个索引 $k$. 下面要做的就是解码, 在做这一点之前, 首先要进行一些概念上的澄清.

在理论上, $|\sqrt{σ}⟩$ 是一个连续的对象, 在空间 $L^2(\mathbb{R}^d, \mathrm{d}x)$ 里面. 但是在门模型量子电路上, 我们处理的是有限维的 $(\mathbb{C}^{2})^{⊗n}≅ \mathbb{C}^{2^n}$. 为了让连续的函数进入我们的量子电路, 需要以下两个步骤:

1. 截断 (truncation): 把 $\mathbb{R}^d$ 截断到一个有限的区域 (通常是一个大立方体) $Ω = [-a,a]^d$ 里面.

2. 离散化 (discretization): 把 $Ω$ 的每一维都划分成 $N=2^q$ 个等距小格, 于是 $Ω$ 被划分成 $N^d$ 个小立方体单元格, 或者网格点, 这时就得到一个有限维空间 $\mathbb{C}^{N^d}$, 从而可以嵌入到 $n=qd$ 个 qubit 的计算空间里面.

这时, 每个索引 $k$ (一个 $n$ 位的二进制字符串) 可以被看作是 $d$ 维空间里面的一个网格点的坐标, 记为 $x_k∈Ω$. 于是我们就可以把 $|g_k|^2$ 看作是位置 $x_k$ 处的概率密度值. 而在这样的计算基下, $|g⟩$ 近似于截断和离散化之后的 $$\begin{equation} |\sqrt{σ}⟩ ≈ \sum_{k} \sqrt{σ(x_k)} |x_k⟩, \end{equation}$$ 从而 $|g_k|^2≈\sigma(x_k)$. 也就是说对于得到的 $|g⟩$ 进行计算基测量, 得到的索引 $k$ 就是从截断和离散化之后的 Gibbs 分布中采样得到的结果, 从而这时测量就自动实现了 Gibbs 采样.