矩阵函数近似中的暴力美学

科学空间|Scientific Spaces 登录 打赏公式天象链接时光博览归档 渴望成为一个小飞侠

欢迎订阅

个性邮箱

天象信息

观测ISS

LaTeX

关于博主

欢迎访问“科学空间”，这里将与您共同探讨自然科学，回味人生百态；也期待大家的分享～

千奇百怪 Everything 天文探索 Astronomy 数学研究 Mathematics 物理化学 Phy-chem 信息时代 Big-Data 生物自然 Biology 图片摄影 Photograph 问题百科 Questions 生活/情感 Life-Feeling 资源共享 Resources 首页 数学研究 矩阵函数近似中的暴力美学 25 Jun 矩阵函数近似中的暴力美学 By 苏剑林 | 2026-06-25 | 2664位读者 |

此前，我们在《矩阵平方根和逆平方根的高效计算》和《矩阵r次方根和逆r次方根的高效计算》探讨过矩阵幂的高效计算，它们将原本用于计算 msign msign 的Newton-Schulz迭代推广到了矩阵幂中。而近日发布到arXiv上的《Muonp: Muon with Fractional Spectral Powers》则提供了另外一个构造迭代来计算矩阵幂的思路，也给了笔者一些启发。

但不管是笔者还是Muonp的思路，总的来说都是“能用，但不够通用”——比如它们只能用来计算矩阵的有理次方，并且复杂度随着最简分数的分子分母增大而增大，这显然不够科学。为了克服这些缺陷，本文提出一套理论上可以拟合任意矩阵函数的通用近似框架。

两种函数 #

说到矩阵函数，实际上它有两种略微不同的含义，这里简单介绍一下。

第一种是对特征值操作的矩阵函数（EIG型）。设矩阵 M 𝑀 的特征值分解为 QΛ Q −1 𝑄 𝛬 𝑄 − 1 ，那么定义 f[M]≜Qf(Λ) Q −1 𝑓 [ 𝑀 ] ≜ 𝑄 𝑓 ( 𝛬 ) 𝑄 − 1 ，其中 f(Λ) 𝑓 ( 𝛬 ) 是指将 f 𝑓 运算逐元素应用到对角线上。这实际上就是我们通常说的矩阵函数，比如矩阵指数、矩阵对数都属于此类，它们通常用幂级数来定义。

第二种是对奇异值操作的矩阵函数（SVD型）。设矩阵 M 𝑀 的奇异值分解为 UΣ V ⊤ 𝑈 𝛴 𝑉 ⊤ ，那么定义 f{M}≜Uf(Σ) V ⊤ 𝑓 { 𝑀 } ≜ 𝑈 𝑓 ( 𝛴 ) 𝑉 ⊤ 。Muon的 msign msign 运算，我们在《通过msign来计算奇异值裁剪mclip（上）》和《通过msign来计算奇异值裁剪mclip（下）》讨论的 mclip mclip 运算，都属于此类。

EIG型矩阵函数只对方阵成立，并且只有在复数域才能保证特征值分解的存在性，所以 f 𝑓 的定义通常还需要延拓到复数（除非将输入矩阵限定在实对称矩阵内）；SVD型矩阵函数对任意形状的矩阵都可以定义，任意实矩阵都有实数范围内的奇异值分解，且奇异值总是非负，这些特性使得SVD型矩阵函数的定义和分析相对简化一些。

这两类矩阵函数各有各的应用场景，它们的计算难度也没有本质不同，但在计算细节上会有所区别。对于幂来说，特征值的任意正整数次幂都可以直接计算，但奇异值我们只能高效计算奇数次幂：

[M ] n = {M } 2n+1 = Q Λ n Q −1 =(QΛ Q −1 ) n = M n U Σ 2n+1 V ⊤ =UΣ V ⊤ (V Σ 2 V ⊤ ) n =M( M ⊤ M ) n (1) (1) [ 𝑀 ] 𝑛 =

𝑄 𝛬 𝑛 𝑄 − 1 = ( 𝑄 𝛬 𝑄 − 1 ) 𝑛 = 𝑀 𝑛

{ 𝑀 } 2 𝑛 + 1 =

𝑈 𝛴 2 𝑛 + 1 𝑉 ⊤ = 𝑈 𝛴 𝑉 ⊤ ( 𝑉 𝛴 2 𝑉 ⊤ ) 𝑛 = 𝑀 ( 𝑀 ⊤ 𝑀 ) 𝑛

即 [M ] n [ 𝑀 ] 𝑛 跟我们用矩阵乘法定义的 M n 𝑀 𝑛 是一致的，而 {M } 2n+1 { 𝑀 } 2 𝑛 + 1 需要借助 ( M ⊤ M ) n ( 𝑀 ⊤ 𝑀 ) 𝑛 来间接计算。所以当我们想要通过多项式来逼近矩阵函数时，EIG型函数可以选择任意次的多项式，SVD型函数则只能选择奇次多项式。

现有方法 #

这一节我们以SVD型分数幂 {M } 1/3 { 𝑀 } 1 / 3 为例，介绍现有的两种近似计算方案。简单起见，假设 M 𝑀 的奇异值都在 [0,1] [ 0 , 1 ] 内。

第一种思路是利用恒等式 {M } 1/3 =M( M ⊤ M ) −1/3 { 𝑀 } 1 / 3 = 𝑀 ( 𝑀 ⊤ 𝑀 ) − 1 / 3 ，然后代入《矩阵r次方根和逆r次方根的高效计算》的结果，得到

G 0 =M, P 0 = M ⊤ M G t+1 = G t ( a t+1 I+ b t+1 P t + c t+1 P 2 t ) P t+1 =( a t+1 I+ b t+1 P t + c t+1 P 2 t ) 3 P t (2) (3) 𝐺 0 = 𝑀 , 𝑃 0 = 𝑀 ⊤ 𝑀

(2) 𝐺 𝑡 + 1 = 𝐺 𝑡 ( 𝑎 𝑡 + 1 𝐼 + 𝑏 𝑡 + 1 𝑃 𝑡 + 𝑐 𝑡 + 1 𝑃 𝑡 2 ) (3) 𝑃 𝑡 + 1 = ( 𝑎 𝑡 + 1 𝐼 + 𝑏 𝑡 + 1 𝑃 𝑡 + 𝑐 𝑡 + 1 𝑃 𝑡 2 ) 3 𝑃 𝑡

最终将有 lim t→∞ G t ={M } 1/3 lim 𝑡 → ∞ 𝐺 𝑡 = { 𝑀 } 1 / 3 。注意 P t 𝑃 𝑡 的迭代是独立的，这个方案的原理是选择适当的 a t , b t , c t 𝑎 𝑡 , 𝑏 𝑡 , 𝑐 𝑡 让 P t →I 𝑃 𝑡 → 𝐼 ，然后单独分离出来的 a t+1 I+ b t+1 P t + c t+1 P 2 t 𝑎 𝑡 + 1 𝐼 + 𝑏 𝑡 + 1 𝑃 𝑡 + 𝑐 𝑡 + 1 𝑃 𝑡 2 的连乘就能趋于 P −1/3 0 𝑃 0 − 1 / 3 。它的主要问题是要先显式计算 M ⊤ M 𝑀 ⊤ 𝑀 ，导致条件数平方，计算精度也随之下降。

第二种思路来自《Muonp: Muon with Fractional Spectral Powers》，它将 m 1/3 𝑚 1 / 3 的计算视为求方程 x 3 −m=0 𝑥 3 − 𝑚 = 0 的根，那么我们只需要设计一个奇次多项式迭代来求该方程的根，就可以写出对应的矩阵迭代。论文考虑的是不动点迭代

x t+1 = x t +c(m− x 3 t )⇔ X t+1 = X t +c(M− X t X ⊤ t X t ) (4) (4) 𝑥 𝑡 + 1 = 𝑥 𝑡 + 𝑐 ( 𝑚 − 𝑥 𝑡 3 ) ⇔ 𝑋 𝑡 + 1 = 𝑋 𝑡 + 𝑐 ( 𝑀 − 𝑋 𝑡 𝑋 𝑡 ⊤ 𝑋 𝑡 )

通过简单分析可以证明，保证任意 m∈[0,1] 𝑚 ∈ [ 0 , 1 ] 都收敛的条件是 c≤2/3 𝑐 ≤ 2 / 3 ，如果取固定值可以考虑 c=2/3 𝑐 = 2 / 3 。这个思路给笔者最大的启发是将输入 M 𝑀 加入到每一步迭代过程中，而不单单依赖于上一步结果 X t 𝑋 𝑡 ，但它仍留有很多的未解决的问题。比如当迭代格式有多参数时，如何逐一确定这些参数，是否可以每一步选择不同的参数，提高收敛速度，这些问题的答案都不得而知。

除此之外，上述两种思路有一个共同的问题，就是计算量强依赖于所求幂次的最简分数形式，比如，如果我们想要求 0.33 0.33 次幂，最简分数是 33/100 33 / 100 ，那么这两种思路都需要我们设法求某个矩阵的 100 100 次方，这是很可观的计算量，但实际上 0.33 0.33 跟 1/3 1 / 3 相差无几，不应该有这么大区别才对。

暴力美学 #

所以，我们需要一个更通用的框架，能够为指定矩阵函数推导一个可用的近似，并且对于相近的矩阵函数能够自动给出相近的结果，最好理论上还能支持幂函数外的任意函数。这便是本文的目标。

不失一般性，考虑SVD型矩阵函数 f{M} 𝑓 { 𝑀 } ，它将奇异值 m 𝑚 变成 f(m) 𝑓 ( 𝑚 ) ，我们需要构造一个奇次多项式迭代来逼近 f(m) 𝑓 ( 𝑚 ) 。假设每一步同时依赖于当前值 x t 𝑥 𝑡 和输入 m 𝑚 （一般地，我们还可以考虑将 x 1 ,⋯, x t−1 𝑥 1 , ⋯ , 𝑥 𝑡 − 1 都加入到迭代中），每一步迭代的阶次不超过3，那么可以构造出一般的迭代格式

x t+1 = c t+1,1 x t + c t+1,2 m+ c t+1,3 x 3 t + c t+1,4 m 3 + c t+1,5 x 2 t m+ c t+1,6 x t m 2 (5) (5) 𝑥 𝑡 + 1 = 𝑐 𝑡 + 1 , 1 𝑥 𝑡 + 𝑐 𝑡 + 1 , 2 𝑚 + 𝑐 𝑡 + 1 , 3 𝑥 𝑡 3 + 𝑐 𝑡 + 1 , 4 𝑚 3 + 𝑐 𝑡 + 1 , 5 𝑥 𝑡 2 𝑚 + 𝑐 𝑡 + 1 , 6 𝑥 𝑡 𝑚 2

其中 c t+1 =( c t+1,1 ,⋯, c t+1,6 ) 𝑐 𝑡 + 1 = ( 𝑐 𝑡 + 1 , 1 , ⋯ , 𝑐 𝑡 + 1 , 6 ) 是待求参数，每步有6个参数，允许每一步都取不同值，以提高近似程度。这里的思路其实很直接，就是把我们能想到的项都加进来，不用强求解释性，由接下来的拟合结果来决定它们是否有用。这样做虽然没有从理论上求得解析解的优雅感，但却有着一种通用计算的暴力美感。

接下来该怎么求 c t+1 𝑐 𝑡 + 1 呢？朴素的做法是像之前的文章《Muon优化器赏析：从向量到矩阵的本质跨越》、《msign算子的Newton-Schulz迭代（上）》一样，先固定一个迭代步数 T 𝑇 ，这样就能将整个迭代视为一个模型，接着再选择一个回归目标，就可以用梯度优化器去端到端训练。

然后，尽管这样做理论上有机会得到更接近全局最优的解，但实践中往往会面临很多的困难——多步迭代后模型的非线性、非凸性，以及初始化的随机性，这都会明显影响效果。

贪心策略 #

为了稳定效果，笔者提出一种逐层求解的贪心策略——每步都假设 c 1 , c 2 ,⋯, c t 𝑐 1 , 𝑐 2 , ⋯ , 𝑐 𝑡 都已知的前提下，去求 c t+1 𝑐 𝑡 + 1 ，优化目标为

c ∗ t+1 = argmin c t+1 d( x t+1 ,f)s.t.∥ c t+1 ∥ ∞ ≤B (6) (6) 𝑐 𝑡 + 1 ∗ = argmin 𝑐 𝑡 + 1 ⁡ 𝑑 ( 𝑥 𝑡 + 1 , 𝑓 ) s.t. ‖ 𝑐 𝑡 + 1 ‖ ∞ ≤ 𝐵

这里 ∥ c t+1 ∥ ∞ ≤B ‖ 𝑐 𝑡 + 1 ‖ ∞ ≤ 𝐵 约束了 c t+1 𝑐 𝑡 + 1 的最大绝对值都不能超过 B 𝐵 ，避免极端放大后缩小导致的精度丢失，我们也可以考虑给每一个 c t+1,i 𝑐 𝑡 + 1 , 𝑖 都选择不同的 B i 𝐵 𝑖 。 d( x t+1 ,f) 𝑑 ( 𝑥 𝑡 + 1 , 𝑓 ) 是成本函数，我们有两种选择

L 2 : L ∞ : ∫ b a ( x t+1 −f(m) ) 2 dm max m∈[a,b] | x t+1 −f(m)| (7) (8) (7) 𝐿 2 : ∫ 𝑎 𝑏 ( 𝑥 𝑡 + 1 − 𝑓 ( 𝑚 ) ) 2 𝑑 𝑚 (8) 𝐿 ∞ : max 𝑚 ∈ [ 𝑎 , 𝑏 ] | 𝑥 𝑡 + 1 − 𝑓 ( 𝑚 ) |

其中 [a,b] [ 𝑎 , 𝑏 ] 是我们关心的奇异值区间。很明显， L 2 𝐿 2 关心的是平均误差， L ∞ 𝐿 ∞ 关心的是最大误差，后面我们将会看到，在这两个选择下我们都可以精确求解。

贪心策略虽然没法保证最优，但也有它的独特优势。首先，Polar Express证明了，如果 f(m)=1 𝑓 ( 𝑚 ) = 1 （即计算 msign msign ），那么贪心解即最优解，这表明在某些情况下贪心解也能很接近最优解；其次，贪心解每一步是精确的、递进的，确保每一步都在更接近目标，我们可以预计算足够多的 c ∗ 𝑐 ∗ ，然后按需截断前面若干步，都是一个很好的近似。

最后，由于贪心解的每一步都在近似目标，所以它的波动还是可控的，不会出现某一步剧烈放大然后下一步再缩小的现象，这对于数值计算尤其是低精度计算尤为重要。我们还可以继续添加约束来实现我们期望的性质，比如 ∥ c t+1 ∥ ∞ ≤B ‖ 𝑐 𝑡 + 1 ‖ ∞ ≤ 𝐵 就是很基本的一个约束，如果有必要，还可以添加更多。

逐一求解 #

之所以选择这两个成本函数，主要原因还是它们对应的优化问题都可精确求解，接下来我们演示它们的求解过程。

首先是 L 2 𝐿 2 ，它需要算一个积分，我们用离散近似后，它就变成了一个有限样本的回归问题，注意到 x t+1 𝑥 𝑡 + 1 关于 c t+1 𝑐 𝑡 + 1 是线性的，所以这只不过是一个线性回归问题，精确可解！即便在此基础上再增加一些线性约束，那也只不过是一个凸二次规划，依然精确可解并且求解工具很成熟。

然后是 L ∞ 𝐿 ∞ ，它需要对整个区间的误差取 max max ，这一点我们也是离散化近似。接下来，我们可以通过一个巧妙的变换，将它转化成一个线性规划问题：引入新变量 z 𝑧 ，我们有

min c t+1 max m∈[a,b] | x t+1 −f(m)|= min c t+1 ,z {z ∣ ∣ −z≤ x t+1 −f(m)≤z,∀m∈[a,b]} (9) (9) min 𝑐 𝑡 + 1 max 𝑚 ∈ [ 𝑎 , 𝑏 ] | 𝑥 𝑡 + 1 − 𝑓 ( 𝑚 ) | = min 𝑐 𝑡 + 1 , 𝑧 { 𝑧 | − 𝑧 ≤ 𝑥 𝑡 + 1 − 𝑓 ( 𝑚 ) ≤ 𝑧 , ∀ 𝑚 ∈ [ 𝑎 , 𝑏 ] }

如果将 [a,b] [ 𝑎 , 𝑏 ] 离散化成 N 𝑁 个点，那么这些样本将转化成关于 7 7 个未知数的 2N 2 𝑁 个线性不等式约束，在这些约束下求 z 𝑧 的最小值。线性规划的求解比二次规划更简单也更成熟，自然也没有什么困难。

能将 L ∞ 𝐿 ∞ 情形转化为线性规划，也是贪心策略的优势之一。如果我们试图用梯度下降类的优化器去联合优化各个 c 𝑐 ，那 L ∞ 𝐿 ∞ 的 max max 将会是优化的核心障碍，因为梯度难以通过 max max 有效传播，加上多步迭代后的函数已经极度非线性，使得求全局解的努力难以奏效。而通过“贪心解+线性规划”的转化，我们可以稳定地获得一个有效的解。

参考实现 #

如果只局限在Numpy和Scipy内，对于带边界约束的线性回归，我们可以用scipy.optimize.lsq_linear来求解，至于一般的线性规划，我们可以用scipy.optimize.linprog来求解。但考虑到代码的通用性和简洁性，以及随时加自定义约束的可能性，我们建议基于CVXPY等专门的凸优化工具来求解。

下面以 f(m)= m 1/3 𝑓 ( 𝑚 ) = 𝑚 1 / 3 为例，给出基于CVXPY的一个参考实现：

import numpy as np import cvxpy as cp

N = 10000 # 离散化点数 B = 10 # 参数边界 m = np.linspace(0, 1.01, N) # 我们关心的奇异值范围是0～1，但参数估计时考虑大一点，以保证稳定性

x, f = m, m(1 / 3) coefs = [] for t in range(10): A = np.array([x, m, x3, m3, x2 * m, x * m**2]).T c = cp.Variable(6) objective = cp.Minimize(cp.sum_squares(A @ c - f)) # L2

objective = cp.Minimize(cp.max(cp.abs(A @ c - f))) # L∞

constraints = [cp.abs(c) <= B] problem = cp.Problem(objective, constraints) result = problem.solve() coefs.append(c.value.round(3)) x = A @ c.value.round(3) print(f'iter {t + 1}, max error:', np.abs(x - f).max()) print(f'iter {t + 1}, mse error:', np.square(x - f).mean())

coefs = np.array(coefs) 效果比较 #

下面我们比较 L 2 𝐿 2 贪心解、 L ∞ 𝐿 ∞ 贪心解与Muonp的固定步长迭代 x t+1 = x t + 2 3 (m− x 3 t ) 𝑥 𝑡 + 1 = 𝑥 𝑡 + 2 3 ( 𝑚 − 𝑥 𝑡 3 ) 在立方根 f(m)= m 1/3 𝑓 ( 𝑚 ) = 𝑚 1 / 3 上的表现。所有方法都从 x 0 =m 𝑥 0 = 𝑚 出发，离散区间为 [0,1.01] [ 0 , 1.01 ] ，参数边界 B=10 𝐵 = 10 ， T 𝑇 为总迭代步数。

L 2 𝐿 2 贪心解迭代10步：最大误差约 8.5× 10 −2 8.5 × 10 − 2 ，均方误差约 5.4× 10 −5 5.4 × 10 − 5 ；

L ∞ 𝐿 ∞ 贪心解迭代10步：最大误差约 4.6× 10 −2 4.6 × 10 − 2 ，均方误差约 8.2× 10 −4 8.2 × 10 − 4 ；

Muonp迭代10步：最大误差约 1.4× 10 −1 1.4 × 10 − 1 ，均方误差约 6.3× 10 −4 6.3 × 10 − 4 ；

Muonp迭代20步：最大误差约 1.0× 10 −1 1.0 × 10 − 1 ，均方误差约 1.3× 10 −4 1.3 × 10 − 4 。

对比图如下：

立方根近似方法对比

从数值上看， L ∞ 𝐿 ∞ 贪心解在最大误差上优势明显， L 2 𝐿 2 次之，固定步长的Muonp方案在 T=20 𝑇 = 20 时，最大误差还不如 T=10 𝑇 = 10 的 L 2 𝐿 2 贪心解；从图上可以看到，Muonp方案的劣势区间主要是0附近，这是因为 c=2/3 𝑐 = 2 / 3 虽然兼顾了 [0,1] [ 0 , 1 ] 整个区间，但对于0附近来说太小了，这是静态步长的主要缺陷。

通过将代码中的 1/3 1 / 3 换成 1/5 1 / 5 ，我们就可以得到5次方根的近似，且迭代阶次不变，但如果用Muonp方案或者笔者之前的 r 𝑟 次方根算法，则至少需要将迭代阶次提升到5，这便是本文框架通用性的体现。

一些结果 #

下面给出求 {M } 0 { 𝑀 } 0 （ msign msign ）、 {M } 1/2 { 𝑀 } 1 / 2 、 {M } 1/3 { 𝑀 } 1 / 3 、 {M } 1/4 { 𝑀 } 1 / 4 的 L ∞ 𝐿 ∞ 贪心解的迭代系数。所有系数均保留三位小数，参数边界 B=10 𝐵 = 10 。其中 m 1/2 , m 1/3 , m 1/4 𝑚 1 / 2 , 𝑚 1 / 3 , 𝑚 1 / 4 的拟合区间为 [0,1.01] [ 0 , 1.01 ] ； m 0 𝑚 0 （即常数1）的拟合区间取 [0.001,1.01] [ 0.001 , 1.01 ] 。

每步迭代格式为

X t+1 = c t+1,1 X t + c t+1,2 M+ c t+1,3 X t X ⊤ t X t + c t+1,4 M M ⊤ M+ c t+1,5 X t X ⊤ t M+ c t+1,6 X t M ⊤ M (10) (10) 𝑋 𝑡 + 1 =

𝑐 𝑡 + 1 , 1 𝑋 𝑡 + 𝑐 𝑡 + 1 , 2 𝑀 + 𝑐 𝑡 + 1 , 3 𝑋 𝑡 𝑋 𝑡 ⊤ 𝑋 𝑡

+ 𝑐 𝑡 + 1 , 4 𝑀 𝑀 ⊤ 𝑀 + 𝑐 𝑡 + 1 , 5 𝑋 𝑡 𝑋 𝑡 ⊤ 𝑀 + 𝑐 𝑡 + 1 , 6 𝑋 𝑡 𝑀 ⊤ 𝑀

其中 X 0 =M 𝑋 0 = 𝑀 ，假设 M 𝑀 的奇异值都已经归一化到 [0,1] [ 0 , 1 ] 内。注意 M ⊤ M 𝑀 ⊤ 𝑀 和 M M ⊤ M 𝑀 𝑀 ⊤ 𝑀 是第一步必须计算的，可以存下来， c t+1,3 X t X ⊤ t X t 𝑐 𝑡 + 1 , 3 𝑋 𝑡 𝑋 𝑡 ⊤ 𝑋 𝑡 与 c t+1,5 X t X ⊤ t M 𝑐 𝑡 + 1 , 5 𝑋 𝑡 𝑋 𝑡 ⊤ 𝑀 可以合并成 X t X ⊤ t ( c t+1,3 X t + c t+1,5 M) 𝑋 𝑡 𝑋 𝑡 ⊤ ( 𝑐 𝑡 + 1 , 3 𝑋 𝑡 + 𝑐 𝑡 + 1 , 5 𝑀 ) ，所以这个迭代虽然有六项，但实际相比Muonp的迭代只增加一步矩阵乘法 X t M ⊤ M 𝑋 𝑡 𝑀 ⊤ 𝑀 。

各系数如下：

{M } 0 {M } 1/2 {M } 1/3 {M } 1/4 t 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 c t,1 2.564 2.668 2.670 2.487 2.350 2.095 1.761 1.547 1.503 1.499 0.931 1.366 1.495 1.261 1.145 1.112 1.071 1.060 1.035 1.027 1.199 1.168 1.897 1.702 1.525 1.392 1.292 1.245 1.173 1.180 1.385 1.398 2.155 1.957 1.789 1.640 1.517 1.405 1.345 1.258 c t,2 2.564 1.243 −0.793 0.018 0.017 0.001 −0.000 −0.001 −0.000 −0.000 0.931 0.368 −0.070 0.402 0.370 0.429 0.535 0.427 0.651 0.506 1.199 1.717 −0.485 −0.140 0.087 0.406 0.384 0.809 0.517 1.091 1.385 1.910 −0.986 −0.663 −0.099 0.020 0.324 0.343 0.651 0.748 c t,3 −1.256 −1.987 −0.791 −0.635 −0.616 −0.582 −0.537 −0.508 −0.501 −0.499 −0.245 −5.209 −5.241 −5.486 −4.675 −4.937 −5.013 −4.697 −5.296 −4.613 −0.405 −4.360 −4.022 −4.187 −4.142 −4.142 −3.870 −4.187 −3.433 −4.355 −0.517 −3.967 −3.272 −3.329 −3.540 −3.559 −3.590 −3.395 −3.514 −3.030 c t,4 −1.256 0.715 3.138 0.033 −0.013 0.001 0.001 0.001 −0.000 0.000 −0.245 0.216 0.379 −0.066 1.648 1.000 1.129 1.472 0.695 1.888 −0.405 0.767 2.809 2.457 2.415 2.750 3.122 3.013 3.409 2.894 −0.517 1.031 3.804 3.625 3.340 3.267 3.190 3.280 3.224 3.640 c t,5 −1.256 7.158 1.727 0.019 −0.003 0.000 0.001 0.001 0.000 0.000 −0.245 10.000 10.000 10.000 10.000 10.000 10.000 10.000 10.000 10.000 −0.405 10.000 10.000 10.000 10.000 10.000 10.000 10.000 9.255 10.000 −0.517 10.000 9.311 9.393 9.381 9.537 9.429 9.260 9.145 8.201 c t,6 −1.256 −10.000 −4.170 −0.060 0.002 −0.002 −0.002 −0.001 0.000 −0.000 −0.245 −5.723 −5.557 −5.143 −7.530 −6.635 −6.762 −7.288 −6.123 −7.839 −0.405 −8.393 −9.167 −8.854 −8.921 −9.500 −10.000 −10.000 −10.000 −9.933 −0.517 −9.554 −10.000 −10.000 −10.000 −10.000 −10.000 −10.000 −10.000 −10.000

𝑡

𝑐 𝑡 , 1

𝑐 𝑡 , 2

𝑐 𝑡 , 3

𝑐 𝑡 , 4

𝑐 𝑡 , 5

𝑐 𝑡 , 6

1

2.564

2.564

− 1.256

− 1.256

− 1.256

− 1.256

2

2.668

1.243

− 1.987

0.715

7.158

− 10.000

3

2.670

− 0.793

− 0.791

3.138

1.727

− 4.170

4

2.487

0.018

− 0.635

0.033

0.019

− 0.060

{ 𝑀 } 0

5

2.350

0.017

− 0.616

− 0.013

− 0.003

0.002

6

2.095

0.001

− 0.582

0.001

0.000

− 0.002

7

1.761

− 0.000

− 0.537

0.001

0.001

− 0.002

8

1.547

− 0.001

− 0.508

0.001

0.001

− 0.001

9

1.503

− 0.000

− 0.501

− 0.000

0.000

0.000

10

1.499

− 0.000

− 0.499

0.000

0.000

− 0.000

1

0.931

0.931

− 0.245

− 0.245

− 0.245

− 0.245

2

1.366

0.368

− 5.209

0.216

10.000

− 5.723

3

1.495

− 0.070

− 5.241

0.379

10.000

− 5.557

4

1.261

0.402

− 5.486

− 0.066

10.000

− 5.143

{ 𝑀 } 1 / 2

5

1.145

0.370

− 4.675

1.648

10.000

− 7.530

6

1.112

0.429

− 4.937

1.000

10.000

− 6.635

7

1.071

0.535

− 5.013

1.129

10.000

− 6.762

8

1.060

0.427

− 4.697

1.472

10.000

− 7.288

9

1.035

0.651

− 5.296

0.695

10.000

− 6.123

10

1.027

0.506

− 4.613

1.888

10.000

− 7.839

1

1.199

1.199

− 0.405

− 0.405

− 0.405

− 0.405

2

1.168

1.717

− 4.360

0.767

10.000

− 8.393

3

1.897

− 0.485

− 4.022

2.809

10.000

− 9.167

4

1.702

− 0.140

− 4.187

2.457

10.000

− 8.854

{ 𝑀 } 1 / 3

5

1.525

0.087

− 4.142

2.415

10.000

− 8.921

6

1.392

0.406

− 4.142

2.750

10.000

− 9.500

7

1.292

0.384

− 3.870

3.122

10.000

− 10.000

8

1.245

0.809

− 4.187

3.013

10.000

− 10.000

9

1.173

0.517

− 3.433

3.409

9.255

− 10.000

10

1.180

1.091

− 4.355

2.894

10.000

− 9.933

1

1.385

1.385

− 0.517

− 0.517

− 0.517

− 0.517

2

1.398

1.910

− 3.967

1.031

10.000

− 9.554

3

2.155

− 0.986

− 3.272

3.804

9.311

− 10.000

4

1.957

− 0.663

− 3.329

3.625

9.393

− 10.000

{ 𝑀 } 1 / 4

5

1.789

− 0.099

− 3.540

3.340

9.381

− 10.000

6

1.640

0.020

− 3.559

3.267

9.537

− 10.000

7

1.517

0.324

− 3.590

3.190

9.429

− 10.000

8

1.405

0.343

− 3.395

3.280

9.260

− 10.000

9

1.345

0.651

− 3.514

3.224

9.145

− 10.000

10

1.258

0.748

− 3.030

3.640

8.201

− 10.000

文章小结 #

本文提出了一种基于贪心策略的通用矩阵函数近似框架——不追求严格的可解释性，而是直接构造简单的多项式迭代，基于贪心策略每一步都去回归目标，然后选择适当的成本函数，转化成二次规划或线性规划问题，精确、稳定地求解迭代参数，最终获得一个有效的计算方案。

转载到请包括本文地址：https://spaces.ac.cn/archives/11787

更详细的转载事宜请参考：《科学空间FAQ》

如果您还有什么疑惑或建议，欢迎在下方评论区继续讨论。

如果您觉得本文还不错，欢迎分享/打赏本文。打赏并非要从中获得收益，而是希望知道科学空间获得了多少读者的真心关注。当然，如果你无视它，也不会影响你的阅读。再次表示欢迎和感谢！

打赏

如果您需要引用本文，请参考：

苏剑林. (Jun. 25, 2026). 《矩阵函数近似中的暴力美学 》[Blog post]. Retrieved from https://spaces.ac.cn/archives/11787

@online{kexuefm-11787, title={矩阵函数近似中的暴力美学}, author={苏剑林}, year={2026}, month={Jun}, url={\url{https://spaces.ac.cn/archives/11787}}, }

分类：数学研究    标签：函数, 迭代, 近似, 矩阵, SVD 5 评论 < 强制间隔投影（Margin-Enforcing Projection） | 让炼丹更科学一些（七）：步长调度与权重平均 > 你也许还对下面的内容感兴趣 矩阵参数的奇异值熵越高越好吗？ DeepSeek V4的tid2eid是怎么来的？ 直接以FID为Loss：从梯度计算到流式训练 如何更科学地估计矩阵的谱范数？ MuP之上：4. 坚守参数的稳定性 基于流式幂迭代的Muon实现：5. 延伸 基于流式幂迭代的Muon实现：4. 原理 基于流式幂迭代的Muon实现：3. 雕琢 基于流式幂迭代的Muon实现：2. 加速 基于流式幂迭代的Muon实现：1. 初识 发表你的看法 追风随影

June 26th, 2026

苏神，你好！如果考虑不动点法求根，参考论文中采用的方法在10步的近似效果不好，是否跟选取的迭代格式有关系？在数值计算中，不动点迭代可以选取不同的格式，不同格式的收敛速度也是不一致的(有线性的，超线性的)，理论上可能存在收敛速度更快的迭代格式？

回复评论 苏剑林 发表于 June 26th, 2026

是存在这种可能性，但你得告诉我构造方式我才能比较呀。

比如对于立方根来说， x t+1 = x t +c(m− x 3 t ) 𝑥 𝑡 + 1 = 𝑥 𝑡 + 𝑐 ( 𝑚 − 𝑥 𝑡 3 ) 算是最简单的迭代格式，确实可以考虑其他更高阶的多项式迭代，但具体怎么构造、系数怎么选择，并没有人给出过一般的结果，我也无从验证。

回复评论 苏剑林 发表于 June 28th, 2026

@追风随影|comment-29522

不好意思，我不小心删了你上一条评论。这里补充回复一下。

你上一条评论，大致的意思就是牛顿法，让 c∝1/m 𝑐 ∝ 1 / 𝑚 ，自适应步长来加速。然而我们最终目的，不是真的求标量 m 𝑚 的 1/3 1 / 3 次方，而是求矩阵 M 𝑀 的 {M } 1/3 { 𝑀 } 1 / 3 ，也就是每个奇异值都取 1/3 1 / 3 次方。这种情况下，我们顶多只能根据最大的奇异值去调 c 𝑐 ，但这样的 c 𝑐 对于剩余奇异值来说依然是次优的。如果要给每个奇异值都设置最优步长，那就要引入逆矩阵运算，这就超出了多项式迭代的范围。

回复评论 追风随影

June 29th, 2026

苏神，你好！非常感谢你提供的问题背景解答，以下是我思路的一个整理：现在的思路消除了矩阵求逆，但可能会引入更高的复杂度，优势是是可以保证不动点迭代的根不变。 设矩阵 M∈ R m×n 𝑀 ∈ 𝑅 𝑚 × 𝑛 的奇异值分解为

M=UΣ V T ,Σ=diag( σ 1 ,…, σ r ),0≤ σ i ≤1. (11) (11) 𝑀 = 𝑈 Σ 𝑉 𝑇 , Σ = diag ⁡ ( 𝜎 1 , … , 𝜎 𝑟 ) , 0 ≤ 𝜎 𝑖 ≤ 1.

建议的迭代格式为

X t+1 = X t +(M− X t X T t X t ) p t ( X T t X t ) . (12) (12) 𝑋 𝑡 + 1 = 𝑋 𝑡 + ( 𝑀 − 𝑋 𝑡 𝑋 𝑡 𝑇 𝑋 𝑡 ) 𝑝 𝑡 ( 𝑋 𝑡 𝑇 𝑋 𝑡 ) .

其中 p t 𝑝 𝑡 是普通多项式，例如

p t (s)= a t,0 + a t,1 s+⋯+ a t,d s d . (13) (13) 𝑝 𝑡 ( 𝑠 ) = 𝑎 𝑡 , 0 + 𝑎 𝑡 , 1 𝑠 + ⋯ + 𝑎 𝑡 , 𝑑 𝑠 𝑑 .

在标量奇异值层面，对应

x i,t+1 = x i,t +( σ i − x 3 i,t ) p t ( x 2 i,t ) . (14) (14) 𝑥 𝑖 , 𝑡 + 1 = 𝑥 𝑖 , 𝑡 + ( 𝜎 𝑖 − 𝑥 𝑖 , 𝑡 3 ) 𝑝 𝑡 ( 𝑥 𝑖 , 𝑡 2 ) .

因此，第 i 𝑖 个奇异值在第 t 𝑡 步的有效步长为

c eff i,t = p t ( x 2 i,t ). (15) (15) 𝑐 𝑖 , 𝑡 e f f = 𝑝 𝑡 ( 𝑥 𝑖 , 𝑡 2 ) .

当 x i,t 𝑥 𝑖 , 𝑡 接近目标值 σ 1/3 i 𝜎 𝑖 1 / 3 时，有

x 2 i,t ≈ σ 2/3 i , (16) (16) 𝑥 𝑖 , 𝑡 2 ≈ 𝜎 𝑖 2 / 3 ,

所以

c eff i,t ≈ p t ( σ 2/3 i ). (17) (17) 𝑐 𝑖 , 𝑡 e f f ≈ 𝑝 𝑡 ( 𝜎 𝑖 2 / 3 ) .

理想情况下，希望

p t (s)≈ 1 3s ,s= σ 2/3 i . (18) (18) 𝑝 𝑡 ( 𝑠 ) ≈ 1 3 𝑠 , 𝑠 = 𝜎 𝑖 2 / 3 .

这样就可以在不使用逆矩阵的情况下，用多项式 p t 𝑝 𝑡 近似每个奇异值的局部最优步长。

取初值

X 0 =M. (19) (19) 𝑋 0 = 𝑀 .

假设某一步有

X t =Udiag( x 1,t ,…, x r,t ) V T . (20) (20) 𝑋 𝑡 = 𝑈 diag ⁡ ( 𝑥 1 , 𝑡 , … , 𝑥 𝑟 , 𝑡 ) 𝑉 𝑇 .

则

X T t X t =Vdiag( x 2 1,t ,…, x 2 r,t ) V T , (21) (21) 𝑋 𝑡 𝑇 𝑋 𝑡 = 𝑉 diag ⁡ ( 𝑥 1 , 𝑡 2 , … , 𝑥 𝑟 , 𝑡 2 ) 𝑉 𝑇 ,

并且

X t X T t X t =Udiag( x 3 1,t ,…, x 3 r,t ) V T . (22) (22) 𝑋 𝑡 𝑋 𝑡 𝑇 𝑋 𝑡 = 𝑈 diag ⁡ ( 𝑥 1 , 𝑡 3 , … , 𝑥 𝑟 , 𝑡 3 ) 𝑉 𝑇 .

又因为

p t ( X T t X t )=Vdiag( p t ( x 2 1,t ),…, p t ( x 2 r,t )) V T (23) (23) 𝑝 𝑡 ( 𝑋 𝑡 𝑇 𝑋 𝑡 ) = 𝑉 diag ⁡ ( 𝑝 𝑡 ( 𝑥 1 , 𝑡 2 ) , … , 𝑝 𝑡 ( 𝑥 𝑟 , 𝑡 2 ) ) 𝑉 𝑇

在由 V 𝑉 张成的右奇异子空间上成立，所以

X t+1 = X t +(M− X t X T t X t ) p t ( X T t X t ) =Udiag ( x i,t +( σ i − x 3 i,t ) p t ( x 2 i,t )) r i=1 V T . (24) (25) (24) 𝑋 𝑡 + 1 = 𝑋 𝑡 + ( 𝑀 − 𝑋 𝑡 𝑋 𝑡 𝑇 𝑋 𝑡 ) 𝑝 𝑡 ( 𝑋 𝑡 𝑇 𝑋 𝑡 ) (25) = 𝑈 diag ⁡ ( 𝑥 𝑖 , 𝑡 + ( 𝜎 𝑖 − 𝑥 𝑖 , 𝑡 3 ) 𝑝 𝑡 ( 𝑥 𝑖 , 𝑡 2 ) ) 𝑖 = 1 𝑟 𝑉 𝑇 .

因此，从 X 0 =M 𝑋 0 = 𝑀 出发，迭代始终保持相同的左右奇异向量。整个矩阵迭代可以分解为一组独立的标量奇异值迭代

x i,t+1 = x i,t +( σ i − x 3 i,t ) p t ( x 2 i,t ). (26) (26) 𝑥 𝑖 , 𝑡 + 1 = 𝑥 𝑖 , 𝑡 + ( 𝜎 𝑖 − 𝑥 𝑖 , 𝑡 3 ) 𝑝 𝑡 ( 𝑥 𝑖 , 𝑡 2 ) .

令目标矩阵为

X ∗ ={M } 1/3 =U Σ 1/3 V T . (27) (27) 𝑋 ∗ = { 𝑀 } 1 / 3 = 𝑈 Σ 1 / 3 𝑉 𝑇 .

则

X ∗ X T ∗ X ∗ =U Σ 1/3 Σ 2/3 V T =UΣ V T =M. (28) (28) 𝑋 ∗ 𝑋 ∗ 𝑇 𝑋 ∗ = 𝑈 Σ 1 / 3 Σ 2 / 3 𝑉 𝑇 = 𝑈 Σ 𝑉 𝑇 = 𝑀 .

于是

M− X ∗ X T ∗ X ∗ =0. (29) (29) 𝑀 − 𝑋 ∗ 𝑋 ∗ 𝑇 𝑋 ∗ = 0.

代入式 (???) (???) ，得到

X t+1 = X ∗ . (30) (30) 𝑋 𝑡 + 1 = 𝑋 ∗ .

因此

X ∗ ={M } 1/3 是迭代的精确不动点。 (31) (31) 𝑋 ∗ = { 𝑀 } 1 / 3 是迭代的精确不动点。

这一点是该结构相对于一般拟合型多项式迭代的重要优势。无论 p t 𝑝 𝑡 如何选择，只要迭代达到目标，残差自动为零，后续迭代不会偏离目标。

当 d=0 𝑑 = 0 时， p 𝑝 是常数。由上述构造可得

p 0 (s)= 2 3(β+μ) . (32) (32) 𝑝 0 ( 𝑠 ) = 2 3 ( 𝛽 + 𝜇 ) .

于是

X t+1 = X t + 2 3(β+μ) (M− X t X T t X t ). (33) (33) 𝑋 𝑡 + 1 = 𝑋 𝑡 + 2 3 ( 𝛽 + 𝜇 ) ( 𝑀 − 𝑋 𝑡 𝑋 𝑡 𝑇 𝑋 𝑡 ) .

局部最坏收敛因子为

ρ 0 = β−μ β+μ . (34) (34) 𝜌 0 = 𝛽 − 𝜇 𝛽 + 𝜇 .

这就是常数步长格式在谱区间 [μ,β] [ 𝜇 , 𝛽 ] 上的最优折中选择。

当 d=1 𝑑 = 1 时，可得

p 1 (s)= 8(β+μ−s) 3( β 2 +6βμ+ μ 2 ) . (35) (35) 𝑝 1 ( 𝑠 ) = 8 ( 𝛽 + 𝜇 − 𝑠 ) 3 ( 𝛽 2 + 6 𝛽 𝜇 + 𝜇 2 ) .

对应的矩阵迭代为

X t+1 = X t + 8 3( β 2 +6βμ+ μ 2 ) (M− X t X T t X t )((β+μ)I− X T t X t ) . (36) (36) 𝑋 𝑡 + 1 = 𝑋 𝑡 + 8 3 ( 𝛽 2 + 6 𝛽 𝜇 + 𝜇 2 ) ( 𝑀 − 𝑋 𝑡 𝑋 𝑡 𝑇 𝑋 𝑡 ) ( ( 𝛽 + 𝜇 ) 𝐼 − 𝑋 𝑡 𝑇 𝑋 𝑡 ) .

在奇异值层面，它对应

x i,t+1 = x i,t + 8(β+μ− x 2 i,t ) 3( β 2 +6βμ+ μ 2 ) ( σ i − x 3 i,t ). (37) (37) 𝑥 𝑖 , 𝑡 + 1 = 𝑥 𝑖 , 𝑡 + 8 ( 𝛽 + 𝜇 − 𝑥 𝑖 , 𝑡 2 ) 3 ( 𝛽 2 + 6 𝛽 𝜇 + 𝜇 2 ) ( 𝜎 𝑖 − 𝑥 𝑖 , 𝑡 3 ) .

该格式不含 σ −1 i 𝜎 𝑖 − 1 ，也不含 x −1 i,t 𝑥 𝑖 , 𝑡 − 1 。矩阵上也不含逆矩阵。

在 s∈[μ,β] 𝑠 ∈ [ 𝜇 , 𝛽 ] 上，它的局部最坏收敛因子为

ρ 1 = (β−μ ) 2 β 2 +6βμ+ μ 2 . (38) (38) 𝜌 1 = ( 𝛽 − 𝜇 ) 2 𝛽 2 + 6 𝛽 𝜇 + 𝜇 2 .

相比常数步长的

ρ 0 = β−μ β+μ , (39) (39) 𝜌 0 = 𝛽 − 𝜇 𝛽 + 𝜇 ,

一次多项式预条件通常能给出更小的局部收敛因子，但它每步需要更多矩阵乘法。

回复评论 苏剑林 发表于 June 29th, 2026

哇，感谢你详尽的分析，对高阶迭代的参数选择给出了可行的分析思路。

不过，总的来说，我感觉这种理论分析还是很有局限性。首先，假设我们只考虑三阶迭代：

x t+1 = c t+1,1 x t + c t+1,2 m+ c t+1,3 x 3 t 𝑥 𝑡 + 1 = 𝑐 𝑡 + 1 , 1 𝑥 𝑡 + 𝑐 𝑡 + 1 , 2 𝑚 + 𝑐 𝑡 + 1 , 3 𝑥 𝑡 3

实测发现， x t+1 = x t +2/3×(m− x 3 t ) 𝑥 𝑡 + 1 = 𝑥 𝑡 + 2 / 3 × ( 𝑚 − 𝑥 𝑡 3 ) 的效果并不如贪心解。叠更高阶自然是可以的，但我们很难相信（虽然没有证据），同样的参数化和计算量下，理论分析得出的解，会比同样阶次的贪心解好。

而且贪心解的优势还在于它可以随意暴力构造，不追求可解释性，我们可以叠加很多小计算量但理论上很难解释的项，来进一步提高精度。从标量的视角看，这些项并不优雅，可能单纯就是把立方根像查表一样背下来，从效率上看可能还得不偿失，但对于矩阵立方根来说，它就有巨大优势，因为矩阵立方根没法查表。

又比如要求五次方根，那么不动点迭代就至少要考虑 x t+1 = x t +c(m− x 5 t ) 𝑥 𝑡 + 1 = 𝑥 𝑡 + 𝑐 ( 𝑚 − 𝑥 𝑡 5 ) 了，但贪心解依然可以考虑三阶迭代。

回复评论

你的大名

电子邮箱

个人网站（选填）

1. 可以使用LaTeX代码，点击“预览效果”可查看效果；
2. 可以通过点击评论楼层编号来引用该楼层；
3. 网站可能会有点卡，如非确认评论失败，请不要重复点击提交。

内容速览 两种函数 现有方法 暴力美学 贪心策略 逐一求解 参考实现 效果比较 一些结果 文章小结 智能搜索 支持整句搜索！网站自动使用结巴分词进行分词，并结合ngrams排序算法给出合理的搜索结果。 热门标签 生成模型 attention 优化 模型 语言模型 梯度 矩阵 优化器 概率 网站 转载 微分方程 分析 天象 深度学习 积分 python 几何 扩散 力学 无监督 节日 损失函数 生活 文本生成 随机文章 三角函数幂的定积分 【备忘】电脑远程控制手机的解决方案 更别致的词向量模型(三)：描述相关的模型 不成功的尝试：将多标签交叉熵推广到“n个m分类”上去 2012年天象 SVD分解(一)：自编码器与人工智能 几个有关集合势的“简单”证明 【NASA每日一图】不规则的NGC 55 关于中国人获得诺贝尔奖的情况 基于流式幂迭代的Muon实现：3. 雕琢 最近评论 zcj5918: 苏神你好，是不是事实上，在，$\log q_{\phi}(\boldsymbol{z}|\bo... linyu: 能稳定复现，方法比较简单：将safetensors模型文件的前三层tid2eid权重以一个概率... 江西老表: 既然要注册，就可能要走它的网站吧？手机和电脑会不会变肉鸡？ 苏剑林: 确实，已经更正，感谢指出 苏剑林: 确实很神奇，能稳定复现吗？随机设置的操作有什么规律呢？对模型效果的影响很小，指的是能正常说话，... xuan1918: (13) 是typo？应该没有积分吧 chenz: 抱歉打扰苏老师了，翻了前面的几页评论，明白您的意思了，确实这个分母对最终的最优解没有影响 chenz: 苏老师好，我的理解是：公式 (12) ( 12 ) 到公式 (15) ( 15 ) 推导出的Loss是让$s_\theta... linyu: 注：改完Deepseek-V4-Flash的tid2eid后，不需要任何训练。 linyu: 苏神您好，我发现一个非常神奇的现象。把Deepseek-V4-Flash的tid2eid完全或... 友情链接 Cool Papers 数学研发 Seatop Xiaoxia 积分表-网络版 丝路博傲 数学之家 有趣天文奇观 TwistedW godweiyang AI柠檬 王登科-DK博客 ESON 枫之羽 coding-zuo 博科园 孔皮皮的博客 运鹏的博客 jiming.site OmegaXYZ EAI猩球 文举的博客 申请链接

本站采用创作共用版权协议，要求署名、非商业用途和保持一致。转载本站内容必须也遵循“署名-非商业用途-保持一致”的创作共用协议。 © 2009-2026 Scientific Spaces. All rights reserved. Theme by laogui. Powered by Typecho. 备案号: 粤ICP备09093259号-1/2。