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}[⟨0|U_{\text{QSP}}(\vec{\Phi };x)|0⟩]=P(x),\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,\dots ,\tilde{d}$，$\tilde{d}\gtrsim 2d$）上定义最小二乘损失

$$\mathcal{L}(\vec{\Phi })=\frac{1}{\tilde{d}}\sum _{j=1}^{\tilde{d}}|\mathrm{Re}⟨0|U_{\text{QSP}}(\vec{\Phi };x_j)|0⟩-P(x_j){|}^2.$$

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

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

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

能力边界。双精度浮点下，这套流程稳定到 $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}=\frac{p_d}{q_{d-1}}$（最高次切比雪夫系数之比）确定最后一个相位 ${\Phi }_d$，且剥除后的 $(P^{\prime },Q^{\prime })$ 是 $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 }$ 应给出完全一致的 $⟨0|U_{\text{QSP}}|0⟩$ 曲线。

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

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

这带来两条实务推论：

1. 误差不是均匀爆炸，而是从端点开始。验证相位时，务必单独检查 $x=\pm 1$ 及其邻域（例如 $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,\epsilon =10^{-3}$。由 minimax 理论，逼近 $\mathrm{sign}(x)$（过渡带 $[\frac{1}{\kappa },1]$）的奇多项式次数 $d\approx \kappa \log (12/\epsilon )\approx 500$ 量级（快速收缩纹波多项式）；
2. 生成目标 $P$（切比雪夫系数形式）；
3. 调用相位求解器（优化法，度数从 $d/8$ 递增）得 $\vec{\Phi }$；
4. 验证：在含 $x=\pm (1-10^{-4})$ 的密网格上计算 $\mathrm{Re}⟨0|U_{\text{QSP}}(\vec{\Phi };x)|0⟩$，比对 $\Vert \cdot -P{\Vert }_{\infty }\le 10^{-12}$（这验证的是”相位正确”，不是”多项式好”——后者由第 1 步的理论保证）；
5. 把 $\vec{\Phi }$ 交给 QSVT 组装层（块编码 + 交替投影序列）。

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

本词条总结

• 求相位 = 解一个”方程数=未知数、解离散有限”的非线性系统；可行性精确为零损失，算法应当收敛到零而不是停在梯度小处。
• 优化法（切比雪夫节点 + 牛顿/LM + 度数递增）简单稳健到 $d\sim 10^3$；确定性构造（补多项式求根/谱分解 + 剥层，或精确 arccos 变体）免初值、可到 $d\sim 10^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)=ZW(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$ 的关系。

---

参考资料

• 本词条整理自《量子计算算法教程》原文：ch05-hamiltonian-qsp/qsp-phase-factors
• 内容遵循 CC BY-NC-SA 4.0 许可协议