有些话得先坦白我第一次搜到的关键词其实是“流行优化”当时还以为是某个新出的网红优化器脑子里全是“流行、潮流、加持”。点进去才发现是“流形”两个字黎曼流形优化。那阵子我手头一个项目正被正交约束折磨——嵌入层每步更新后都要强行投影回正交矩阵用SVD修慢、不稳训练曲线还在原地抖。知道这世界上有“沿着曲面本身走”的优化方法之后我花了两周从数学书啃到手写代码这篇记录就是完整复盘。适合刚接触流形优化、或者正被各种矩阵约束搞到崩溃的同学参考我会尽量把概念说成“人话”代码也能直接抄。1. 为什么优化要跑到流形上约束才是真主角1.1 约束条件带来的麻烦远比你想象的大大多数人对“约束优化”的印象还停留在x 0或者a x b这种盒式约束上。可一旦约束变成矩阵层面的等式事情就开始离谱了。举个最常见的例子要求参数矩阵 (W) 满足 (W^T W I)列正交这在PCA、正交RNN、Lipschitz约束网络、去相关特征里到处都是又比如要求协方差矩阵是正定对称的SPD矩阵要求低秩分解里的因子保持低秩。这些约束集合根本不是向量空间也不是凸集它们本质上是一张“弯弯曲曲的曲面”。我见过不少人的第一反应是先正常梯度下降一步再把结果暴力投影回约束集合。投影到正交矩阵就做SVD投影到单位范数就做归一化。听起来很直接但实际一跑就翻车。原因有几个投影本身可能很贵。每步都做一次SVD计算量直接和优化主循环一个量级训练速度肉眼可见往下掉。这些集合不是凸的投影之后不能保证跳到“合理”的位置甚至可能在两个不连通的子集里来回横跳。更隐蔽的问题是欧氏空间里“距离最近”的投影在流形上不一定是最好的方向。你以为走的是捷径实际是在曲面外面穿行每一步都在犯错。1.2 与其把约束当外挂不如把约束当空间本身流形优化的核心思路很反直觉既然解必须落在某个曲面上那就干脆把优化过程定义在曲面内部让每一步都被迫“在面上走”而不是在面外乱跑。拿三维空间里的球面举个例子。球面上任意一点它的“可行方向”不是所有三维向量而是那些与径向垂直的向量——也就是球面上过该点的切线方向。你可以沿着球面走但不能穿透球面。如果我们直接用三维梯度去更新下一步必然飞到球面外于是你必须再投影回球面。流形优化换个玩法先把梯度投影到该点的切空间再在切空间里走一步然后用一个“缩回”retraction映射回球面。仔细想想这不就是在球面坐标系里做优化吗对更深层的矩阵流形也一样。正交矩阵集合Stiefel流形、单位模长向量集合球流形、正定对称矩阵集合SPD流形、固定秩矩阵集合都有各自的“切空间”和“缩回方式”。一旦你把这套东西配齐优化算法反而变得比“梯度下降暴力投影”更干净。不用拉格朗日乘子不用频繁做昂贵的投影收敛行为也稳定得多。所以我一直觉得流形优化的本质不是什么高深数学就是一句话把空间的形状真正尊重起来。2. 切空间、黎曼梯度与retraction三个绕不开的概念2.1 切空间曲面上的“局部线性化”为什么要先讲切空间因为它就是流形上“能走的方向”的集合。流形在局部看起来像欧氏空间但每个点的“平坦方向”都不一样。切空间就是这些方向构成的向量空间维度等于流形的维度。以最常见的三个流形为例它们各自的切空间长这样球面 (S^{n-1})(x \in \mathbb{R}^n, |x|1)。切向量 (v) 必须满足 (x^T v 0)即所有与 (x) 垂直的向量。Stiefel流形 (\mathrm{St}(p,n))(X \in \mathbb{R}^{n \times p}, X^T X I_p)。切向量 (V) 必须满足 (X^T V V^T X 0)也就是 (X^T V) 是个反对称矩阵。这一点在写代码时特别容易出错。SPD流形 (\mathrm{Sym}_n^)正定对称矩阵集合。切空间就是所有对称矩阵维度是 (n(n1)/2)。有了切空间优化问题的每步更新就有一个天然的“活动范围”。在这个范围内找方向就不会轻易跑出流形。很多初学者的第一个坑就是拿欧氏空间梯度的分量直接当切向量用结果完全不符合流形的几何关系。2.2 黎曼度量与黎曼梯度同一个函数不同的“陡峭”有了方向还不够还要知道“哪条路最陡”。在欧氏空间里负梯度就是最陡下降方向。但流形上每个点的切空间自带一个内积定义这个内积就是黎曼度量。不同度量下同一个函数的“最陡方向”可以完全不一样。怎么理解你在平地上看一座山最陡的方向由地面坐标决定但如果换一套度量比如在某个点附近放大了特定方向的距离那“最陡”就会变。黎曼流形优化的关键一步就是把欧氏梯度投影到切空间再按度量修正得到黎曼梯度。大多数情况下如果流形嵌入在欧氏空间里且用的度量是自然诱导的度量比如Frobenius内积黎曼梯度就是“欧氏梯度的切空间正交投影”。球面上很简单[ \mathrm{grad}{\text{man}} f(x) \mathrm{grad}{\text{euclid}} f(x) - (x^T \mathrm{grad}_{\text{euclid}} f(x)) x ]Stiefel流形稍微复杂一点设欧氏梯度为 (G)黎曼梯度公式是[ \mathrm{grad}_{\text{man}} G - X \cdot \mathrm{sym}(X^T G), \quad \mathrm{sym}(M) \frac{M M^T}{2} ]这个 “sym” 的来历不是拍脑袋它保证了投影后的方向仍然满足切空间的反对称条件。实际写代码时很多人漏掉后面那项的一半结果优化方向噪声极大收敛奇慢。2.3 指数映射、retraction 与向量传输方向和步长都有了怎么把“走一步”变成合法的流形上的点数学上最优雅的方式是指数映射沿测地线走。但指数映射多数情况下算起来很贵于是工程上几乎都用 retraction缩回/回缩映射。它的要求很简单在切点处等于原位置且一阶导等于切向量本身。满足这两条就能保证收敛性不被破坏。不同流形有不同的retraction我平常用得最多的是下面这几个流形Retraction 公式备注球面(R_x(v) \frac{xv}{|xv|})一维归一化最简单Stiefel(R_X(V) \mathrm{qf}(XV))对 (XV) 做QR取Q因子Stiefel(R_X(V) (XV)((XV)^T (XV))^{-1/2})Polar型数值更稳SPD流形(R_X(V) X V)步长小可近似但要小心正定性SPD流形(R_X(V) X^{1/2} \exp(X^{-1/2} V X^{-1/2}) X^{1/2})仿射不变度量下的精确指数映射除了 retraction还有一个“向量传输”vector transport概念。动量法、共轭梯度、BFGS需要把上一步的方向搬到当前点但流形不是平直的向量不能直接搬。向量传输就是平行移动的廉价替代品。实现也很简单最常见做法是把要传输的向量先按某种规则搬运到新点附近再做一次切空间投影。新手阶段如果只跑梯度下降可以先不管它但一旦上动量就躲不开了。3. 从RGD到Pymanopt流形优化算法与工具选型3.1 黎曼梯度下降把欧氏套路“翻译”到流形上黎曼梯度下降RGD的更新规则看起来和普通梯度下降几乎一样x_{k1} retract(x_k, -step * grad_man(x_k))差别就在两个词梯度必须是黎曼梯度位移必须用retraction而不是直接相加。就这么简单。换到代码里一次迭代的核心步骤只有四步在当前点计算欧氏梯度 (G)把 (G) 投影到切空间得到黎曼梯度取负方向乘步长得到切向量 (V)调用retraction把 (xV) 映射回流形。如果步长不固定想用线搜索也只需要把目标函数看成一维函数 (\phi(t) f(\mathrm{retract}(x, t \cdot V)))然后做标准的Armijo回溯。这比在欧氏空间做线搜索还要省事因为retraction把每一步都“粘”在流形上函数值有正常下降的几何保证。实践里我很推荐先手写一遍RGD再去用库。因为这个算法总共不到三十行写一遍能逼着你把切空间、黎曼梯度、retraction三件事彻底理清。我见过直接用库的人遇到问题只会换参数连“是不是梯度投影错了”都判断不了。3.2 二阶方法与自适应优化动量、BFGS、Trust RegionRGD本质一阶方法慢在收敛尾巴。想快可以上二阶思想但每个都要在流形上多加几件装备动量法需要在每次更新后把历史动量用向量传输搬到新切空间。若是长跑训练动量带来的加速很明显但传输写错还不如不加。共轭梯度CG需要计算当前切向量与前一步搜索方向的“传输版本”的夹角同样依赖向量传输。黎曼牛顿法 / Trust Region需要计算流形上的Hessian工程实现复杂但收敛速度是二次的适合中小规模、高精度的优化问题。Manopt库里的trustregions是这类算法的代表作。自适应优化RAdam等Geoopt库实现了流形版本的Adam梯度方向和矩统计都在切空间中进行矩的移动同样需要向量传输或近似处理。我的使用经验是小规模科学计算问题直接上Trust Region或黎曼牛顿收敛快且不挑步长大规模深度学习训练多半用RGD配合小步长或RAdam中间地带的矩阵补全、度量学习这类问题用RGD线搜索往往性价比最高。3.3 现成工具怎么选Pymanopt、Manopt、Geoopt自己实现是学原理真正做项目我更推荐站在巨人肩膀上。三套常用工具各有脾气工具语言定位一句话点评PymanoptPython科学计算、快速原型配合自动微分写代价函数即可梯度自动帮你算ManoptMATLAB / Julia经典库算法全文档丰富适合对照论文验证GeooptPython PyTorch深度学习嵌入能直接在PyTorch优化器里换成RAdam等Pymanopt 是我最常用的学习伴侣。它要求你定义流形和代价函数然后通过装饰器自动生成梯度。注意Pymanopt的自动梯度生成结果与欧氏梯度一致内部再用流形的egrad2rgrad转成黎曼梯度。所以即使数学推导不在行也能先把算法跑通再用数值验证自己的理解。Geoopt 则在神经网络里方便它的Stiefel流形封装了参数的正交约束前向传播不用你手动投影。不过工具永远替代不了理解。我在章节4里会给出一个完全不用第三方流形库、只用numpy和scipy跑通的正交PCA例子那才是把原理焊在脑子里最好的方式。4. 手写Stiefel流形上的PCA一份可直接跑的代码4.1 问题建模为什么PCA会变成流形优化假设我们有数据协方差矩阵 (\Sigma)想找到前 (p) 个主成分方向组成的矩阵 (W)那么优化目标是[ \min_{W^T W I_p} ; f(W) -\mathrm{tr}(W^T \Sigma W) ]这里 (W) 落在Stiefel流形 (\mathrm{St}(p,n)) 上。目标函数就是让投影后的方差尽量大约束保证列正交。如果 (pn)解就是特征向量矩阵的某种旋转(p1) 时退化为球面上的瑞利商最大化。先算欧氏梯度。对 (f(W) -\mathrm{tr}(W^T \Sigma W)) 求方向导数得到[ G -2 \Sigma W ]这里假设 (\Sigma) 对称求梯度时不用再对称化。接下来的任务就变成把 (G) 投影到Stiefel流形的切空间选一个retraction循环更新。4.2 完整实现切空间投影、QR回缩与线搜索下面这份代码能直接运行只依赖 numpy 和 scipyimport numpy as np from scipy.linalg import qr def euclidean_grad(W, Sigma): # 对 f(W) -tr(W^T Sigma W) 求欧氏梯度 return -2.0 * Sigma W def riemannian_grad(W, G): # 把欧氏梯度投影到 Stiefel 流形切空间 # 切空间条件: W^T V V^T W 0 A W.T G symA (A A.T) / 2.0 return G - W symA def retract_qr(X): # QR 回缩取 QR 分解的 Q 因子 Q, R qr(X, modeeconomic) # 可选对齐列符号避免 QR 输出符号跳变 sign np.sign(np.diag(R)) sign[np.abs(sign) 1e-12] 1.0 return Q np.diag(sign) def objective(W, Sigma): return -np.trace(W.T Sigma W) # 构造一个对称半正定的协方差矩阵 np.random.seed(42) n, p 10, 3 A np.random.randn(n, n) Sigma A.T A np.eye(n) * 1e-6 # 随机初始化正交矩阵 W0, _ qr(np.random.randn(n, p)) W W0.copy() max_iter 500 tol 1e-10 step 0.1 prev_loss objective(W, Sigma) for it in range(max_iter): G euclidean_grad(W, Sigma) grad riemannian_grad(W, G) gnorm np.linalg.norm(grad, ordfro) if gnorm tol: break # 固定步长更新也可以用回溯线搜索 V -step * grad W_next retract_qr(W V) # 简单的回溯线搜索Armijo loss objective(W_next, Sigma) armijo_c 1e-4 while loss prev_loss armijo_c * step * np.sum(grad * V): step * 0.5 V -step * grad W_next retract_qr(W V) loss objective(W_next, Sigma) if step 1e-12: break W W_next prev_loss loss if it % 50 0: print(fiter {it:4d} | loss {loss:.6f} | grad_norm {gnorm:.4e}) # 用标准特征分解核对 eigvals, eigvecs np.linalg.eigh(Sigma) true_W eigvecs[:, -p:] # 最大的 p 个特征向量 # 比较子空间投影矩阵 W W^T 的差异 proj_mine W W.T proj_true true_W true_W.T print(子空间误差:, np.linalg.norm(proj_mine - proj_true, ordfro))几个容易写错的地方我特意标一下riemannian_grad里的symA必须是对称化后的 (X^T G)漏掉(A A.T)/2.0里的转置项方向就会带噪声。QR回缩一般直接用qr(X, modeeconomic)的Q。我在后面补了符号对齐是为了避免QR每轮列方向符号随意翻转导致目标函数在相邻两步之间出现不自然的跳变。实际训练中如果发现损失曲线有锯齿优先检查这一处。线搜索里用的是 Frobenius 内积grad, V因为流形度量是欧氏内积诱导的。如果换到SPD流形上的仿射不变度量内积公式也要跟着换。4.3 跑通之后你会发现三个有趣现象第一个现象用固定步长也能收敛但步长太大时损失曲线会在最后几步出现小的振荡。这不是算法错了而是retraction只保证一阶近似大方向下QR回缩与理想测地线的误差被放大了。线搜索可以有效压住这个振荡。第二个现象最终得到的 (W) 和标准PCA的 (W) 往往差一个正交旋转。因为目标函数只在乎子空间不在乎基的顺序所以不要直接比较 (W) 的每一列要比 (WW^T) 这个投影矩阵。做流形优化时这种“排列/旋转不确定性”会频繁出现在各种问题里提前有心理准备能省掉大量排查时间。第三个现象最开始的几步目标函数下降特别快后面梯度范数会缓慢爬向0。这是流形优化的典型收敛模式和欧氏空间一样所以别指望曲线像直线一样掉到0。5. 踩坑实录收敛震荡、正交性漂移与库的脾气5.1 梯度/retraction写错用流形上的差分检验我犯过最愚蠢的错误是把切空间投影里的sym写成了0.5 * A而不是0.5 * (A A.T)。当时损失曲线虽然也在降但方向一直不对收敛慢得离谱。后来学到一个很有用的自查方法在流形上做有限差分。取一个随机切向量 (V)数值上计算t 1e-6 num_grad (objective(retract_qr(W t * V), Sigma) - objective(retract_qr(W - t * V), Sigma)) / (2 * t) ana_grad np.sum(riemannian_grad(W, G) * V)理想情况下两者应该非常接近。如果差得远要么切空间投影错了要么retraction的一阶性质不满足。这个小测试只需要几行代码却能帮你定位绝大多数“玄学”问题建议写入你自己的测试模块里。5.2 步长与收敛别再拿欧氏最优步长来赌很多从欧氏优化转过来的朋友会对固定步长特别执着。确实在欧氏二次型上最优步长有解析公式但流形上每一点曲率不同同样步长在不同位置的实际“前进距离”不一样。所以我建议小规模问题直接上回溯线搜索不要手动调步长。大规模深度学习场景若想要训练稳定步长宁可小一些也别指望一个步长走天下。如果必须用自适应步长注意RAdam这类流形自适应优化器内部有额外的切空间向量传输改超参数时要看训练曲线的平滑度而不是只看最终损失。5.3 库和框架的暗坑QR回缩、GPU与SPD流形用PyTorch做批量Stiefel优化时很多人会踩torch.linalg.qr的坑。批量QR默认返回的是“完整”分解内存开销大正确做法是指定modereduced或者直接改用极分解型retraction先torch.linalg.svd再取U Vt。SVD在GPU上通常比批量QR稳定反传也相对可控。SPD流形是另一个重灾区。如果你直接在欧氏空间里做X V很可能下一步就失去正定性然后出现NaN。用我前面给过的指数映射retraction虽然计算贵一点但至少每一步都合法。另外SPD流形上不同度量差别很大欧氏度量下做均值、仿射不变度量下做均值、对数欧氏度量下做均值结果都不一样。这不是数值误差是几何定义本来就不一样项目文档里一定要写清楚到底用了哪个度量。6. 自学路线与资料清单照着学不迷路6.1 数学准备够用就好别一头扎进微分几何流形优化最大的学习门槛是“心理门槛”。我最初也以为要先把光滑流形、纤维丛、测地线全啃完才能动手后来发现做落地应用根本不需要那么重。你只需要几块知识流形的定义、切空间、黎曼度量、指数映射与retraction、向量传输。至于更深的内容用到再补。如果时间有限我建议按这个顺序读Absil、Mahony、Sepulchre 的《Optimization Algorithms on Matrix Manifolds》2008年的书例子多适合当查字典。Boumal 的《An Introduction to Optimization on Smooth Manifolds》2023年新书更现代算法讲得更清楚网上有电子版和配套讲义。Lee 的《Introduction to Smooth Manifolds》不作为第一本只当几何概念的补充查询。不要一上来就背定义。拿第4章的代码跑通一遍再回头看书里的公式效率会高很多。6.2 由易到难的练习顺序与验收标准我给自己的练习清单是这样的每个都能在半天内完成验收标准也明确球面上的瑞利商最大化目标 (\max_{x^T x1} x^T A x)。手写RGD和特征分解结果对比重点练切空间投影和归一化retraction。Stiefel流形上的PCA就是本文的代码。验收标准是投影矩阵 (WW^T) 与标准PCA误差在 (10^{-6}) 级别。Grassmann流形上的子空间估计重点体验Grassmann流形与Stiefel流形的差别很多子空间类问题用它更自然。SPD流形上的黎曼均值给一堆正定矩阵用仿射不变度量求黎曼均值和欧氏均值做对比理解“度量变了结论就变”。每个练习都配上第5.1节的数值梯度校验。这样每一个流形你都能形成“定义、切空间、retraction”三行速查表后续套算法就是换皮。6.3 最后一点个人心得这套东西真正难的不是公式是思维的转向。以前我遇到正交约束条件反射是“先无约束更新再投影”现在我会先问自己这个约束集合是不是流形切空间里能不能找到方向有没有便宜的retraction如果三个问题都能答上来剩下的就是把熟悉的算法翻译到流形上。要说有什么最想提醒后来人的就一句话不要迷信库。Pymanopt、Geoopt写得再漂亮也只能帮你算数值不能帮你建立直觉。我第三次重写Stiefel代码时才真正理解为什么切空间条件长那样也才明白为什么流形优化比“SVD投影”稳那么多。希望你不用等第三次才懂。