QSP 相位因子的数值计算详解：从优化问题到确定性算法

前置阅读：量子信号处理（QSP）详解——本篇接续其"第二步：多项式到相位参数的转换"，把那一步从"存在性"推进到"怎么算"。

课程目标：

  1. 把"给定目标多项式 \(P\)，求相位向量 \(\vec\Phi\)"明确为一个非线性方程组/优化问题，并理解它为什么不是"随便交给 scipy 就能解好"的。

  2. 掌握优化法（Haah 2019；Dong–Meng–Whaley–Lin 2021）：切比雪夫节点采样、度数递增、牛顿迭代的病态与对策。

  3. 掌握确定性构造的思路：补多项式（complementary polynomial）求根 + 逐层剥除（Fejér–Riesz 谱分解路线），以及精确 arccos 变体。

  4. 理解数值病态的根源（端点附近的相位"挤压"）与高精度算术的出场时机。

  5. 会用现成工具（QSPPACK、pyqsp、Qualtran 子程序）跑通一个小例子并自行验证。

1. 问题的形状：\(d+1\) 个未知数、\(d+1\) 个方程、以及它的坏脾气

回顾 QSP 详解的记号：\(d\) 阶序列 \(U_{\text{QSP}}(\vec\Phi;x) = S(\Phi_0)\prod_{k=1}^d W(x)S(\Phi_k)\)，其 \((0,0)\) 元素是 \(P(x)\)。给定目标多项式（例如逼近 \(e^{-i\tau x}\)、\(\mathrm{sign}(x)\) 或 \(1/(\kappa x)\) 截断式的 \(P\)），求相位问题是：

\[ \text{找到 } \vec\Phi \in \mathbb{R}^{d+1} \text{ 使得 } \mathrm{Re}\big[\langle 0|U_{\text{QSP}}(\vec\Phi;x)|0\rangle\big] = P(x),\quad \forall x \in [-1,1]. \]

对比两边切比雪夫展开 \(P = \sum_k p_k T_k\)：系数方程恰好 \(d+1\) 个（受奇偶性约束后减半），未知相位 \(d+1\) 个（对称化后同样减半）——方程数与未知数相当，且解集是离散有限的（QSP 详解中的"相位解离散性"命题）。这既说明问题"刚好可解"，也说明它不是凸优化：损失地形有大量对称的局部结构，通用求解器从随机初值出发很容易陷入"系数对了模长、错了相位"的陷阱。

数值上还有三重坏脾气：目标 \(P\) 常常是minimax 最优多项式（通带内有等幅纹波），可行解处损失函数的 Hessian 病态；高次（\(d \gtrsim 500\)）时相位对扰动极其敏感；约定不统一（\(W(x)\) 与 \(e^{i\theta Z}\) 两套信号约定、相位排在左还是右）让不同工具的输出不可直接混用——这是复现文献数值时的头号踩坑点。

2. 方法一：优化法——把存在性当成零点去找

优化法的关键洞察是：可行性是精确的。若 \((P, Q)\) 是可行对（满足 \(|P|^2 + (1-x^2)|Q|^2 = 1\)），则存在 \(\vec\Phi\) 使损失严格为零——我们找的不是"近似最优"，而是一个已知的零点。这决定了算法形态：用迭代法收敛到零，而不是停在梯度小的地方。

损失函数。在切比雪夫节点 \(x_j = \cos\frac{(2j-1)\pi}{4\tilde d}\)（\(j = 1, \ldots, \tilde d\)，\(\tilde d \gtrsim 2d\)）上定义最小二乘损失

\[ \mathcal L(\vec\Phi) \;=\; \frac{1}{\tilde d} \sum_{j=1}^{\tilde d} \Big| \mathrm{Re}\,\langle 0|U_{\text{QSP}}(\vec\Phi;x_j)|0\rangle - P(x_j) \Big|^2 . \]

选切比雪夫节点不是习惯问题：它们是切比雪夫插值的最优采样点，能使离散损失以机器精度逼近连续 \(\ell_2\) 损失，避免在 \(x = \pm1\) 附近的信息欠采样。

迭代与雅可比。\(\langle 0|U_{\text{QSP}}|0\rangle\) 对每个 \(\Phi_k\) 的偏导有精确公式（对 \(S(\Phi_k)\) 求导即乘一个 \(iZ\)，把 \(iZ\) 插进乘积链的对应位置），因此可以用牛顿法/Levenberg–Marquardt 而无需数值微分。雅可比在解附近病态——但方向"好"：奇异值谱呈几何衰减，衰减方向恰与"整体相位平移"这类无害冗余对齐。实践对策是 LM 阻尼 + 在残差进入 \(10^{-8}\) 后切换到更高精度的线性求解。

度数递增（progressive phase factor）。直接对大 \(d\) 冷启动必败。正确姿势是阶梯式：先解 \(d_{\text{small}}\)（如 \(d=2,4,8,\ldots\) 或按问题尺度），把低次解补零相位后作为下一次的初值。这一策略惊人地有效——低次的相位结构携带了高频纹波位置的"骨架"，逐级加倍通常一路收敛到底。

能力边界。双精度浮点下，这套流程稳定到 \(d \sim 10^2\) 量级毫无压力、到 \(d \sim 10^3\) 需要小心实现（阻尼、递增、对称化），再往上就会撞墙——原因见第 4 节。此时要么换方法二，要么用 mpmath 之类的多精度算术。

3. 方法二：确定性构造——剥层与求根

优化法是"猜–验证"，确定性构造是"按图纸拆"。它源自 Gilyén 论文与 Haah (2019) 的递归算法，核心是 QSP 详解"理论推导"一节的剥层引理与补多项式：

第一步：求补多项式 \(Q\)。 可行性方程 \(|P|^2 + (1-x^2)|Q|^2 = 1\) 把 \(Q\) 的系数钉死（在 \(Q(0)\) 相位选择的分歧意义下唯一）。把 \(x = \cos\theta\) 代入并用 \(T_k(\cos\theta) = \cos k\theta\)，方程变成三角多项式的恒等式；在其"半边谱"上做 Fejér–Riesz 谱分解（把非负三角多项式分解为 \(|Q(e^{i\theta})|^2\)）即得 \(Q\)。

第二步：逐层剥除。 剥层引理说：若 \((P, Q)\) 是 \(d\) 层可行对，则 \(e^{2i\Phi_d} = \dfrac{p_d}{q_{d-1}}\)（最高次切比雪夫系数之比）确定最后一个相位 \(\Phi_d\)，且剥除后的 \((P', Q')\) 是 \(d-2\) 层可行对。于是从最高层往下剥，每层解一个辐角方程——每一步只做一次除法与一次 arctan，总共 \(O(d)\) 次剥除。加上求 \(Q\) 的谱分解，整体复杂度 \(O(d^2)\)（朴素实现），快速多项式乘法可压到近线性。

精确 arccos 变体。 上述除法 \(p_d/q_{d-1}\) 在数值上偶尔需要取模长（理论值为 1，浮点下有抖动）。后续工作（Chao 等 2020 及相关改进）把"每层一次剥除"重写为"每层一次精确 \(\arccos\)"的稳定递推，把可稳定计算的次数推到 \(d \sim 10^4\) 量级，成为高次数场景的标准做法。

方法对比：

优化法

确定性构造

初值

需要度数递增策略

免初值，直接从系数出发

复杂度

每次牛顿 \(O(d^2)\)，共 \(O(\log d)\) 级

\(O(d^2)\)（可加速）

稳定上限（双精度）

\(d \sim 10^3\)

\(d \sim 10^4\)

灵活性

可加正则化、软约束（如固定若干相位）

只解决标准问题

实现门槛

中（雅可比易错）

高（谱分解、剥层细节多）

实践建议：小中规模（\(d \lesssim 500\)）用优化法（实现简单、可交叉验证），大规模或需要批量计算时用确定性构造，两者互为对方的验证器——同一 \(P\) 两法得到的 \(\vec\Phi\) 应给出完全一致的 \(\langle 0|U_{\text{QSP}}|0\rangle\) 曲线。

4. 病态从哪来：端点挤压与高精度算术

为什么高次相位难算？几何图像是端点挤压：\(W(x)\) 的两条谱支在 \(x = \pm1\) 处简并（\(W(\pm1) = \pm I\)，间隙关闭）。要把 \(\langle 0|U_{\text{QSP}}|0\rangle\) 在端点附近也钉在目标多项式上，相邻相位被迫彼此靠近——对 \(e^{-i\tau x}\) 这类目标，最优相位序列中相邻相位差可以小到 \(O(\rho^{-d})\)（\(\rho\) 为 Bernstein 椭圆参数，见 QSP 详解"逼近精度"节）。\(d\) 越大、椭圆越贴近 \([-1,1]\)，挤压越狠：相位的有效数字被指数级消耗。

这带来两条实务推论：

  1. 误差不是均匀爆炸，而是从端点开始。验证相位时，务必单独检查 \(x = \pm1\) 及其邻域（例如 \(x = 1 - 10^{-4}\)），不要只看网格中段。

  2. 超过双精度能力时换多精度算术，而不是加迭代。\(d \gtrsim 10^3\!\sim\!10^4\) 后牛顿迭代"不收敛"通常不是初值问题，而是浮点分辨率不够——用 mpmath（Python）把工作精度提到 \(O(d)\) 位十进制，同一套算法立即恢复收敛。

另一个工程细节：对称化。对偶/对称目标（绝大多数应用），相位满足 \(\Phi_k = \Phi_{d-k}\)，把自由参数减半——既降低优化维度，也让确定性构造的剥层成对进行，数值更稳。多数工具默认输出对称化序列，读文献表格时注意对方是否做了这一步。

5. 工具与一个小例子

  • QSPPACK（Dong & Lin，MATLAB）：确定性构造 + 优化法的参考实现，文档里有一系列应用（模拟、求逆、Gibbs）的目标多项式与相位表。

  • pyqsp（Python 礏区实现）：pyqsp.angle_sequence 一行求相位，适合快速原型；约定与 QSPPACK 略有差异，混用时先做端点对齐检查。

  • Qualtran / pyLIQTR：容错资源估计框架里的 QSP 子程序，把相位求解、块编码组装、代价统计串成流水线。

五步走一个小例子（求逆滤波）：目标 \(\mathrm{sign}\) 函数近似（求逆的核心部件）。

  1. 设 \(\kappa = 10,\ \varepsilon = 10^{-3}\)。由 minimax 理论，逼近 \(\mathrm{sign}(x)\)（过渡带 \([\tfrac1\kappa, 1]\)）的奇多项式次数 \(d \approx \kappa \log(12/\varepsilon) \approx 500\) 量级（快速收缩纹波多项式）；

  2. 生成目标 \(P\)（切比雪夫系数形式）；

  3. 调用相位求解器（优化法，度数从 \(d/8\) 递增）得 \(\vec\Phi\)；

  4. 验证：在含 \(x = \pm(1 - 10^{-4})\) 的密网格上计算 \(\mathrm{Re}\langle 0|U_{\text{QSP}}(\vec\Phi;x)|0\rangle\)，比对 \(\| \cdot - P\|_\infty \le 10^{-12}\)（这验证的是"相位正确"，不是"多项式好"——后者由第 1 步的理论保证）；

  5. 把 \(\vec\Phi\) 交给 QSVT 组装层（块编码 + 交替投影序列）。

第 4 步的自验证（chebfun 风格地在网格上算函数值）成本忽略不计，却能在几秒内抓住约定错、奇偶错、挤压爆炸三类最常见事故——请把它当作固定动作。

本课总结

  • 求相位 = 解一个"方程数=未知数、解离散有限"的非线性系统；可行性精确为零损失，算法应当收敛到零而不是停在梯度小处。

  • 优化法（切比雪夫节点 + 牛顿/LM + 度数递增）简单稳健到 \(d\sim10^3\)；确定性构造（补多项式求根/谱分解 + 剥层，或精确 arccos 变体）免初值、可到 \(d\sim10^4\)，两者互为验证器。

  • 病态根源是端点谱简并导致的相位挤压：有效数字随 \(d\) 指数消耗；对策是端点邻域专项检查与多精度算术，而非盲目迭代。

  • 生态：QSPPACK / pyqsp / Qualtran；注意约定差异与对称化状态。

  • 本篇补全了 QSP→QSVT 流水线的"最后一公里"：从"存在相位"到"把相位算出来、验正确"。

习题

  1. 手推 \(d = 2\) 的系数方程：设 \(P(x) = p_0 + p_2(2x^2 - 1)\)（偶），写出 \(U_{\text{QSP}}\) 的 \((0,0)\) 元素与 \((\Phi_0+\Phi_1+\Phi_2,\ \Phi_0-\Phi_1+\Phi_2)\) 的关系，对给定 \((p_0, p_2)\)（如 \(\cos(2x)\) 的截断）解出三个相位并数值验证。

  2. 实现题：用自动微分（如 JAX）实现 \(\mathcal L(\vec\Phi)\) 与其梯度，对 \(P = T_8\) 从零初值跑 LM；再实现"度数从 2 递增到 64"策略，比较收敛率。

  3. 证明对称化断言：若 \(P\) 是偶（奇）多项式且可行，则存在 \(\Phi_k = \Phi_{d-k}\) 的解（提示：对 \(U_{\text{QSP}}\) 取 \(x \to -x\)，利用 \(W(-x) = Z\, W(x)\, Z\) 并把 \(Z\) 吸收进相位重参数化）。

  4. 定量感受挤压：取 \(P_d \to e^{-i\tau x}\)（\(\tau = 50\)），用任一工具求 \(d = 64, 128, 256\) 的相位，画相邻相位差 \(|\Phi_{k+1} - \Phi_k|\) 的最小值随 \(d\) 的变化；解释其与 Bernstein 椭圆参数 \(\rho \approx d/\tau\) 的关系。