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

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

**课程目标：**

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 详解](qsp-tutorial.md)的记号：$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 详解](qsp-tutorial.md)"理论推导"一节的**剥层引理**与**补多项式**：

**第一步：求补多项式 $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$ 的关系。