Polar Express：用 Chebyshev 逼近把 Muon 的矩阵正交化提速一倍

在 Muon 优化器：矩阵正交化驱动的梯度更新中，我们建立了 msign 算子的数学骨架：把梯度矩阵 $G$ 投影到最近的正交矩阵 $\text{msign}(G) = UV^T$ ，并用 Newton-Schulz 迭代避免完整的 SVD。这套方案已经在 Kimi K2 上实现了 2× 训练加速，但它有一个隐疾：Newton-Schulz 在矩阵奇异值跨度大时，前十几步几乎不动。

ICLR 2026 Honorable Mention 论文 The Polar Express: Optimal Matrix Sign Methods and Their Application to the Muon Algorithm (Amsel et al.) 把这个问题钉到了棺材里。他们证明：对任意奇异值区间 $[\ell, u]$ ，存在唯一的区间最优 odd quintic 多项式——而 Newton-Schulz 的 $\frac{3}{2}x - \frac{1}{2}x^3$ 不是这个最优解。把每一步的多项式都换成对应区间的最优解，复合下来在 $T$ 步 deg-5 多项式空间中全局最优。

实测上，Polar Express 比 Newton-Schulz 快约 2 倍达到机器精度；在 GPT-2 Large + FineWeb 1B tokens 的训练中，val loss 在所有 learning rate 上一致优于 Newton-Schulz 和此前流行的 Jordan / You 启发式。本文从 Newton-Schulz 的痛点出发，把 Jordan 启发式、You 六步法、Polar Express 三家解法放在同一个数学框架下对比，详解其背后的 Chebyshev 一致逼近 + Remez 算法，并逐行解读 80 行开源实现。

痛点：Newton-Schulz 的"前 17 步几乎不动"

Newton-Schulz 迭代

Z_{k+1} = \frac{3}{2}Z_k - \frac{1}{2}Z_k Z_k^T Z_k

在标量上等价于 $f(x) = \frac{3}{2}x - \frac{1}{2}x^3$ 。它是把所有奇异值"推向 1"的一个标量映射，作用于每个奇异值独立。

固定点分析： $f(x) = x \Rightarrow x = 0, \pm 1$ 。在 $(0, 1)$ 内 $f(x) > x$ ，所以奇异值单调上升趋于 1。问题出在收敛速度上。 $f$ 在 $x = 1$ 附近的 Taylor 展开

f(1 - \epsilon) = 1 - \frac{3}{2}\epsilon^2 + O(\epsilon^3)

二次收敛——靠近 1 时极快。但在 $x \to 0$ 时

f(\epsilon) \approx \frac{3}{2}\epsilon

只有 1.5 倍的线性增长。如果某个奇异值是 $\sigma = 10^{-6}$ （接近实际深度学习梯度矩阵的最小奇异值），要把它推到 $\sigma \approx 1$ 需要 $\log_{1.5}(10^6) \approx 34$ 步线性增长——而每步矩阵乘法都很贵。

论文 Figure 4 用合成矩阵（ $\sigma \in [10^{-6}, 1]$ 均匀对数分布）展示：Newton-Schulz 前 17 步几乎不动，第 18 步之后才进入二次收敛区域，加起来需要 25-30 步。

这就是问题：经典数值分析时代追求"双精度收敛到机器精度"，迭代次数本来就多；深度学习时代每步都要走 GPU 矩阵乘，前 17 步浪费的 FLOPs 就是真金白银。

修补尝试：Jordan 启发式与 You 六步法

社区早就意识到这个痛点。Jordan（Muon 原作者）在博客里给出一个度-5 多项式

p_{\text{Jordan}}(x) = 3.4445\,x - 4.7750\,x^3 + 2.0315\,x^5

代替 Newton-Schulz 的 $\frac{3}{2}x - \frac{1}{2}x^3$ 。这个系数是用近似梯度下降在某个奇异值分布上优化出来的——本质是启发式调参。它在前几步比 Newton-Schulz 快得多，但代价致命：收敛性丢了。论文 Figure 4 显示 Jordan 迭代到第 11 步左右就卡死在误差 $\approx 0.3$ 不动，无论再迭代多少次都过不去。

You 进一步把每一步用不同的 deg-5 系数，称作"六步法"。这是 Jordan 思路的扩展，但仍然不收敛——只是把 plateau 的高度从 0.3 压低一点，最终上限和 Jordan 同级。

两家的共同问题：用启发式调参代替严格的最优逼近。他们直觉对了（区间最优多项式应该自适应），但缺一个数学工具——能严格回答"给定区间 $[\ell, u]$ ，对常数 1 的最优 odd quintic 逼近是什么"。

这个工具是 Chebyshev 一致逼近理论里的 Remez 算法。

区间最优多项式：Chebyshev 与 Remez

一致逼近问题

设 $f: [\ell, u] \to \mathbb{R}$ 是目标函数（这里 $f \equiv 1$ ）， $\mathcal{P}_d^{\text{odd}}$ 是次数 $\leq d$ 的奇多项式空间（即 $\text{span}\{x, x^3, x^5, \ldots, x^d\}$ ， $d$ 为奇数）。 $L^\infty$ 意义下的最佳逼近是

p^\star = \arg\min_{p \in \mathcal{P}_d^{\text{odd}}} \max_{x \in [\ell, u]} |f(x) - p(x)|

为什么是 $L^\infty$ 而非 $L^2$ ？因为我们要的是任意奇异值都收敛——最差情况上界比平均误差更重要。

等振荡定理（Equioscillation Theorem）

Chebyshev 1854 年的关键定理： $p^\star$ 是 $f$ 的 $L^\infty$ 最佳次数 $\leq d$ 多项式逼近，当且仅当误差函数 $f - p^\star$ 在 $[\ell, u]$ 上至少 $d + 2$ 个点处取得 $\pm \max\|f - p^\star\|$ 且符号交替。

直观图像：最优逼近误差曲线像一条"波浪"，在区间内的极值点处误差正负交替触顶。如果某段误差比其他段小，那一定不是最优——我们可以重新分配多项式系数让所有极值更均衡。

对 odd quintic $p(x) = ax + bx^3 + cx^5$ ， $\dim \mathcal{P}_5^{\text{odd}} = 3$ ，所以等振荡要求 $d + 2 = 4$ 个交替点。设这 4 个点为 $\ell, q, r, u$ （边界 + 两个内部极值），则有：

\begin{aligned} 1 - p^\star(\ell) &= +E \\ 1 - p^\star(q) &= -E \\ 1 - p^\star(r) &= +E \\ 1 - p^\star(u) &= -E \end{aligned}

写成矩阵形式（把 $(a, b, c, E)$ 视为未知数）：

\begin{bmatrix} \ell & \ell^3 & \ell^5 & 1 \\ q & q^3 & q^5 & -1 \\ r & r^3 & r^5 & 1 \\ u & u^3 & u^5 & -1 \end{bmatrix} \begin{bmatrix} a \\ b \\ c \\ E \end{bmatrix} = \begin{bmatrix} 1 \\ 1 \\ 1 \\ 1 \end{bmatrix} \tag{1}

只要知道极值点位置 $(q, r)$ ，解一次 4×4 线性方程就拿到 $(a, b, c)$ 。但极值点位置本身依赖于 $(a, b, c)$ ——是个鸡生蛋问题。

Remez 算法：交换迭代

Remez (1934) 给出了解法：

初始化 一组候选极值点（论文用四分位 $q_0 = (3\ell + u)/4, r_0 = (\ell + 3u)/4$ ）循环：
用当前 $(q_k, r_k)$ 解线性系统 (1) 得 $(a, b, c, E_k)$
计算误差函数 $1 - p(x)$ 的真实极值点位置： $p'(x) = a + 3bx^2 + 5cx^4 = 0$ ，令 $y = x^2$ 解二次方程得新的 $(q_{k+1}, r_{k+1})$
收敛判据： $|E_{k+1} - E_k| < \epsilon$ 则停

实际中通常 3-5 轮就收敛。optimal_quintic 的 Python 实现只有 20 行：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
def optimal_quintic(l, u):
    if 1 - 5e-6 <= l / u:
        return (15/8)/u, (-10/8)/(u**3), (3/8)/(u**5)  # near-converged shortcut
    q = (3*l + u) / 4
    r = (l + 3*u) / 4
    E, old_E = inf, None
    while not old_E or abs(old_E - E) > 1e-15:
        old_E = E
        LHS = np.array([
            [l, l**3, l**5, 1],
            [q, q**3, q**5, -1],
            [r, r**3, r**5, 1],
            [u, u**3, u**5, -1],
        ])
        a, b, c, E = np.linalg.solve(LHS, np.ones(4))
        q, r = np.sqrt((-3*b + np.array([-1, 1]) *
                        sqrt(9*b**2 - 20*a*c)) / (10*c))
    return float(a), float(b), float(c)

完美对应上面三个步骤：LHS 是等振荡线性系统，np.linalg.solve 是步骤 1，sqrt(...) 一行是步骤 2 的二次方程显式解（ $p'(x) = 0$ 化为 $y = x^2$ 的二次方程，闭式根公式）。

近似收敛的 shortcut（if 1 - 5e-6 <= l / u）值得一说：当区间已经收缩到 $[u(1-5\times 10^{-6}), u]$ ，多项式逼近的数值条件数会爆炸，Remez 迭代失效。论文证明此时极限多项式恰好是 $\frac{15}{8}x - \frac{10}{8}x^3 + \frac{3}{8}x^5$ （即标准 deg-5 Newton-Schulz）经 $u$ 缩放后的形式，直接返回。

区间复合：从单步最优到全局最优

单步最优只解决"给定区间 $[\ell, u]$ 怎么找最好的多项式"，但 Muon 迭代是要复合 $T$ 个多项式。一个自然的问题：贪心地每步取当前区间的最优多项式，复合下来还是不是全局最优？

直觉上不一定——贪心算法在很多问题上得不到全局最优。但 Polar Express 论文的 Theorem 4.1 给出了肯定回答：

Theorem 4.1（区间最优的递归结构）。设 $\ell_1 = \ell, u_1 = u$ 。定义
$p_t = \arg\min_{p \in \mathcal{P}_d^{\text{odd}}} \max_{x \in [\ell_t, u_t]} |1 - p(x)|$
并令 $\ell_{t+1} = p_t(\ell_t)$ ， $u_{t+1} = 2 - \ell_{t+1}$ 。则
$p^\star = p_T \circ p_{T-1} \circ \cdots \circ p_1$
在所有 $T$ 个 deg- $d$ 奇多项式复合中全局最优。

证明的核心引理是：在区间 $[\ell_t, u_t]$ 上最优的 $p_t$ ，其最小值恰在 $\ell_t$ 处，最大值在某个内部点 $u_{t+1}$ 处。由等振荡对称性 $p_t(\ell_t) + p_t(u_{t+1}) = 2$ 给出 $u_{t+1} = 2 - p_t(\ell_t)$ 。这意味着：

计算下一轮区间不需要矩阵运算，只需 1 次标量求值 $\ell_{t+1} = p_t(\ell_t)$ 。

收敛速率（Theorem 4.3）：对 $d = 2q + 1$ ，

\|\text{polar}(M) - X_T\|_2 \leq |1 - \ell^2|^{(q+1)^T}

deg-3（ $q=1$ ）二次收敛，deg-5（ $q=2$ ）三次收敛——这是 Newton-Schulz 也具有的渐近阶。Polar Express 的优势不在渐近，而在早期收敛：从 $\ell \approx 10^{-6}$ 出发不再"前 17 步几乎不动"。

把 Remez 求解和区间递推串起来就是 Algorithm 1 的离线阶段（optimal_composition）：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
def optimal_composition(l, num_iters, safety_factor_eps=0, cushion=0):
    u = 1
    coefficients = []
    for iter in range(num_iters):
        a, b, c = optimal_quintic(max(l, cushion*u), u)
        # ... safety factor + cushion rescaling (see below)
        coefficients.append((a, b, c))
        l = a*l + b*l**3 + c*l**5   # Theorem 4.1 区间递推
        u = 2 - l                    # 等振荡对称
    return coefficients

l = a*l + b*l**3 + c*l**5 一行直接对应 $\ell_{t+1} = p_t(\ell_t)$ ，纯标量运算，毫秒级。整个离线阶段用 float64 跑一次得到 $T$ 个系数三元组 $(a, b, c)$ ，烤进 coeffs_list 全局变量，线上推理直接查表用。

三家解法横向对照

把 Jordan 启发式、You 六步法、Polar Express 放在同一框架下比较：

方法	每步多项式	收敛性	早期表现	渐近阶	bfloat16 友好性
Newton-Schulz (deg-5)	固定 $p(x) = \frac{15}{8}x - \frac{10}{8}x^3 + \frac{3}{8}x^5$	√ 全局收敛	极慢（前 17 步几乎不动）	三次	需注意
Jordan 启发式	固定 $(3.4445, -4.7750, 2.0315)$	× 卡在误差 0.3	快	—	未验证
You 六步法	6 个固定 deg-5 多项式	× 仍不收敛	快（前 6 步）	—	未验证
Polar Express (deg-5, $T=8$ )	8 个区间最优多项式（Remez 求解）	√ Theorem 4.1 保证最优	快（5-6 步到机器精度）	三次	√ 显式工程设计

三家的本质分歧：

Newton-Schulz：用"无脑通用"多项式——对任何区间都用同一个。代价是初期慢，因为没有针对当前区间优化。
Jordan / You：意识到要"自适应"，但用启发式调参——快但失去收敛保证。卡在 0.3 不动是因为某些奇异值落在多项式的"死区"。
Polar Express：用严格的最优逼近理论。每一步多项式都是当前区间 $[\ell_t, u_t]$ 上的 $L^\infty$ 最佳逼近，且组合后全局最优。

这里有个有趣的辩证：Newton-Schulz 的"固定多项式"是数值分析时代的产物——那时假设输入分布未知，必须用对任意输入都好的多项式。深度学习时代输入分布是已知的（梯度矩阵的奇异值有典型分布），可以预先针对这个分布优化。Polar Express 的本质是把输入分布的先验知识编进多项式系数。

苏剑林在 msign 算子的 Newton-Schulz 迭代（下）中也指出了同一方向：

Newton-Schulz 迭代的优化已经形成了一个小的研究方向，主要思路是利用输入分布的先验信息来调整多项式系数。

Polar Express 给出了这个方向上的"理论最优"答案。Jordan 和 You 是这个方向上的早期尝试，Polar Express 是定论。

bfloat16 工程取舍：safety factor 与 cushion

理论最优在 float64 下完美，搬到 bfloat16（深度学习训练的实际精度）会出意外。论文 Section 4.4 给出三个工程干预：

Safety factor

理论上每步要求奇异值始终留在 $[\ell_t, u_t]$ 内。但 bfloat16 的 rounding 可能让某个 $\sigma > u_t$ 一点点，落到 $[u_t, u_t + \epsilon]$ 区间。在那里，多项式 $p_t$ 超出了它的"设计区域"，会把误差放大而非缩小。Nakatsukasa & Higham (2012) 给出的常数：每次 rounding 误差 $\epsilon$ 经过 $p_t$ 后被放大到约 $25\epsilon$ ，多步累积就发散。

干预方法：把每个 $p_t$ 替换为 $p_t(x / 1.01)$ ，等价于把"安全区间"略微往下搬一点。代码实现极简：

1
2
3
4
if iter < num_iters - 1:
    a /= safety_factor          # safety_factor = 1.01
    b /= safety_factor**3       # 注意三次方
    c /= safety_factor**5       # 五次方

为什么不同次幂的除数不同？因为 $p(x/sf) = \frac{a}{sf}x + \frac{b}{sf^3}x^3 + \frac{c}{sf^5}x^5$ 。最后一步（iter == num_iters - 1）不应用 safety factor——它不需要为下一步预留误差余量。

Cushion（自由下界保护）

等振荡多项式有一个数值病态行为：当 $\ell_t / u_t$ 极小（比如 $10^{-6}$ ）时，多项式在 $u_t$ 端的"压缩比" $p(u_t) / u_t$ 会非常接近 0。这意味着大奇异值会被压得极小，下一轮的 $u_{t+1}$ 就变得非常接近 $\ell_{t+1}$ ——区间在 1 附近过度收缩，加上 bfloat16 的精度损失，后续迭代基本失效。

cushion 是个软限制：当 $\ell_t < c \cdot u_t$ （论文用 $c = 0.02$ ）时，Remez 求解用区间 $[c \cdot u_t, u_t]$ 而非 $[\ell_t, u_t]$ ，然后重缩放系数让 $p(\ell_t) + p(u_t) = 2$ （恢复等振荡对称）：

1
2
3
4
5
if cushion*u > l:
    pl = a*l + b*l**3 + c*l**5
    pu = a*u + b*u**3 + c*u**5
    rescaler = 2 / (pl + pu)
    a *= rescaler; b *= rescaler; c *= rescaler

代价是稍微减慢收敛——但论文说不显著。

预生成系数表

最后一个工程取舍：所有 $T = 10$ 个多项式系数在程序加载时一次性算好，固化进 coeffs_list 全局变量。线上推理只跑矩阵乘法循环：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
coeffs_list = optimal_composition(
    l=1e-3, num_iters=10,
    safety_factor_eps=1e-2, cushion=0.02
)

@torch.compile
def PolarExpress(G: torch.Tensor, steps: int) -> torch.Tensor:
    X = G.bfloat16()
    if G.size(-2) > G.size(-1): X = X.mT      # FLOPs 优化：让 X 总是宽矩阵
    X = X / (X.norm(dim=(-2, -1), keepdim=True) * 1.01 + 1e-7)
    hs = coeffs_list[:steps] + list(repeat(coeffs_list[-1], steps - len(coeffs_list)))
    for a, b, c in hs:
        A = X @ X.mT
        B = b * A + c * A @ A
        X = a * X + B @ X    # X <- aX + b X X^T X + c (XX^T)^2 X = aX + bX^3 + cX^5
    if G.size(-2) > G.size(-1): X = X.mT
    return X

这意味着 Polar Express 在线推理的每步开销与 Newton-Schulz 完全相同——3 次小矩阵乘 + 1 次大矩阵乘。唯一新增的是查表读 $(a, b, c)$ ，可以忽略。

Horner 形式的等价性值得一验。设 $X = U\Sigma V^T$ 为 SVD，则

XX^T = U\Sigma^2 U^T, \quad XX^TX = U\Sigma^3 V^T = X(X^TX) = U\Sigma V^T \cdot V\Sigma V^T

所以 $X X^T X$ 在 SVD 意义下就是把 $X$ 的奇异值立方（同时保留奇异向量）。同理 $(XX^T)^2 X$ 是奇异值五次方。因此

aX + bXX^TX + c(XX^T)^2X = U(a\Sigma + b\Sigma^3 + c\Sigma^5)V^T

恰是把多项式 $p(x) = ax + bx^3 + cx^5$ 作用到每个奇异值上。

FLOPs 计数：设 $X \in \mathbb{R}^{m \times n}$ 且 $m \leq n$ （转置后保证）， $X X^T$ 是 $m \times m$ 小矩阵。每步 3 次 $m \times m$ 矩阵乘 + 1 次 $m \times m \times n$ 矩阵乘，约 $(2m + 2n)m^2$ FLOPs，与 Newton-Schulz 持平。

实验：GPT-2 上的实战

论文 Section 5 给出三组实验：合成矩阵收敛、真实梯度收敛、GPT-2 端到端训练。

合成矩阵（论文 Figure 4）： $\sigma$ 在 $[10^{-6}, 1]$ 对数均匀分布，左右奇异向量随机。结果：

迭代步数	Newton-Schulz	Jordan	You (6 步)	Polar Express ( $\ell=10^{-3}$ )
5 步	误差 $\approx 0.95$	$\approx 0.4$	$\approx 0.35$	$\approx 0.3$
10 步	$\approx 0.7$	卡在 0.3	$\approx 0.3$	$\approx 0.05$
15 步	$\approx 0.05$	卡在 0.3	卡在 0.3	机器精度
25 步	机器精度	卡在 0.3	卡在 0.3	机器精度

最反直觉的发现：Jordan 启发式即使迭代 100 步也卡在误差 0.3——它从一开始就不收敛。Polar Express 第 15 步已经到机器精度。

GPT-2 Large + FineWeb 1B tokens（论文 Figure 1）：

优化器	最佳 val loss	学习率
AdamW	4.172	1e-4
Muon (Newton-Schulz)	—	—
Muon (Jordan)	3.398	2e-2
Muon (You)	3.400	2e-2
Muon (Polar Express)	3.340	2e-2

差距看似不大（3.340 vs 3.398 = 0.058），但这是在所有 learning rate 上稳定的优势。论文 Figure 6 把 lr 横扫一遍，Polar Express 在每个 lr 点都领先。

实际工程意义：

减少 Muon 内 msign 子例程的迭代步数（论文实验用 5 步，业界常规是 5-7 步）
在 lr 容忍度上有更宽的甜区——重要，因为 lr 是大模型训练最不希望调的超参
训练曲线更稳定，方差更小

后续工作 / 与其他方向的接口

Polar Express 把"Muon 子例程优化"这条线推到了理论最优。但它不是 Muon 工程化的终点——还有几条平行的方向：

流式幂迭代（苏剑林在基于流式幂迭代的 Muon 实现中讨论的方案）：不是在每步训练中跑完整的 Newton-Schulz / Polar Express，而是把 SVD 估计平摊到训练循环里——本质是把多项式迭代用幂迭代替代。两条路线互补：Polar Express 优化"单次正交化"，流式幂迭代优化"训练循环里的均摊成本"。
MuP（Maximal Update Parametrization）框架：Muon 在 MuP 下天然满足缩放稳定性，使得 hyperparameter 可以从小模型迁移到大模型。Polar Express 不改变 Muon 的 MuP 性质，因为它只优化 msign 子例程，不动外层。
流形上的最速下降：苏剑林在流形上的最速下降系列中把 Muon 视为 Stiefel 流形上的优化器——Polar Express 是其投影算子的更高效实现，没有改变流形结构本身。
Halley 迭代 / QDWH 等经典数值分析方法：它们渐近收敛阶更高（Halley 是三次，QDWH 是六次），但依赖矩阵求逆或 QR 分解——GPU 不友好。Polar Express 在 GPU 时代的胜利是保留了"只用矩阵乘"的特性，同时获得了区间最优收敛速度。

小结

Polar Express 的核心贡献用一句话概括：用 Chebyshev 一致逼近 + Remez 算法，把 Newton-Schulz 的"固定多项式"换成"区间最优多项式"，并证明贪心区间复合在 $T$ 步 deg-5 多项式空间中全局最优。这个改动只增加了离线一次性的多项式求解成本（毫秒级），线上推理 FLOPs 与 Newton-Schulz 完全相同，但前期收敛速度快约 2 倍。

从 Newton-Schulz → Jordan / You → Polar Express 的演进，展示了一条非常清晰的研究路径：

观察痛点（NS 前 17 步几乎不动）
启发式尝试（Jordan 调参快但不收敛）
发现根本工具（Chebyshev / Remez 是这个问题的数学语言）
严格最优化（Polar Express）

每一步都是站在前人肩膀上的小步前进，但累积下来是从"启发式 + 失败"到"理论最优 + GPT-2 端到端胜利"的质变。这也是数学工具进入工程的典型方式：以经典理论（Chebyshev 1854）的翻译为载体，进入新的工程语境（bfloat16 GPU）。

对深度学习优化器设计的启示：矩阵正交化子例程的最优化已经接近天花板（理论最优 + 工程实现完备）。下一个 2× 加速大概率来自其他方向——流式幂迭代均摊、动量与正交化的联合设计、或者跳出 Muon 框架的全新优化器范式。

参考文献

Amsel, N., Persson, D., Musco, C., & Gower, R. M. (2025). The Polar Express: Optimal Matrix Sign Methods and Their Application to the Muon Algorithm. ICLR 2026 (Honorable Mention). arXiv:2505.16932
代码: github.com/NoahAmsel/PolarExpress
苏剑林. msign 算子的 Newton-Schulz 迭代（上）. https://kexue.fm/archives/10922
苏剑林. msign 算子的 Newton-Schulz 迭代（下）. https://kexue.fm/archives/10996
苏剑林. 流形上的最速下降：3. Muon + Stiefel. https://kexue.fm/archives/11221
苏剑林. 基于流式幂迭代的 Muon 实现：1. 初识. https://kexue.fm/archives/11654
Jordan, K., et al. (2024). Muon: An optimizer for hidden layers in neural networks. arXiv preprint.
Chen, J., & Chow, E. (2014). A stable scaling of Newton-Schulz for improving the sign function computation of a Hermitian matrix. Argonne National Laboratory Tech Report ANL/MCS-P5059-0114.
Nakatsukasa, Y., & Freund, R. W. (2016). Computing fundamental matrix decompositions accurately via the matrix sign function in two iterations. SIAM Review, 58(3).
Nakatsukasa, Y., & Higham, N. J. (2012). Backward stability of iterations for computing the polar decomposition. SIAM J. Matrix Anal. Appl., 33(2).
Higham, N. J. (2008). Functions of Matrices: Theory and Computation. SIAM.
Achieser, N. I. (1956). Theory of Approximation. Frederick Ungar Publishing.
Trefethen, L. N. (2013). Approximation Theory and Approximation Practice. SIAM.

痛点：Newton-Schulz 的"前 17 步几乎不动"#

修补尝试：Jordan 启发式与 You 六步法#

区间最优多项式：Chebyshev 与 Remez#

一致逼近问题#

等振荡定理（Equioscillation Theorem）#

Remez 算法：交换迭代#

区间复合：从单步最优到全局最优#

三家解法横向对照#

bfloat16 工程取舍：safety factor 与 cushion#

Safety factor#

Cushion（自由下界保护）#

预生成系数表#

实验：GPT-2 上的实战#

后续工作 / 与其他方向的接口#

小结#

相关概念#

参考文献#