Skip to content

算法Algorithm

极分解与 Newton 正交化

Polar decomposition · Newton iteration for the polar factor

证明可逆实矩阵的极分解及最近正交矩阵性质,以奇异值标量迭代解释 Newton 正交化,并区分正交因子与受限旋转。

形式陈述 ​

一个可逆实矩阵既改变长度,也改变方向。如何把这两种作用拆开,并找出离它最近的正交矩阵?

对非奇异 A∈Rn×n,存在唯一分解

A=UH,UTU=I,H=HT≻0.

它称为 极分解,其中 H 是对称正定伸缩,U 是正交因子。具体地,H=(ATA)1/2、U=AH−1。若 奇异值分解为 A=PΣQT,则

U=PQT,H=QΣQT.

唯一性来自 H2=ATA 的正定平方根唯一,再由 U=AH−1 确定方向因子。

最近正交矩阵 ​

U 是 Frobenius 范数下唯一的最近正交矩阵:

‖A−U‖F=minZTZ=I‖A−Z‖F,‖A−U‖F2=∑i(σi−1)2.

证明时令 W=PTZQ,则 W 正交,并且

‖A−Z‖F2=‖Σ‖F2+n−2tr(ΣW).

每个 wii≤1,而所有 σi>0,所以迹在 W=I 时达到最大。达到最大要求每个 wii=1,正交性进一步迫使 W=I,故最近点唯一。

Newton 正交化 ​

不显式求 SVD 时,可运行

X0=A,Xk+1=12(Xk+Xk−T).

实现为解 XkTYk=I,再令 Xk+1=(Xk+Yk)/2。这是一种 固定点迭代;非奇异输入下,精确运算中的所有迭代都有定义,并收敛到 U。

证明保持相同的左右奇异向量:Xk=PΣkQT,每个正奇异值独立更新为

σk+1=12(σk+σk−1),σk+1−1=(σk−1)22σk.

第一轮后 σk≥1,随后单调趋于 1,接近目标时误差满足二次上界;若初值已正交,轨道从一开始就保持不变。因此 Xk→PQT。

直觉

H 沿一组正交轴拉伸或压缩,U 再整体改变方向。Newton 更新保留这组左右方向,专门把每个奇异值推向一:过大的奇异值与其倒数取平均,过小的奇异值也被倒数拉回。它不是逐列单独归一化。

图中原始单位圆被送成椭圆;正交因子只改变方向,不改变长度。下方两条奇异值轨道在第一轮后都落到 1 上方,再向同一个目标靠拢。

例子与边界

一个剪切与伸缩混合的矩阵 ​

取

A=(2101).

它的极分解为

U=110(31−13),H=110(6224).

直接检查 UTU=I、UH=A,且 H 的特征值为正。前两轮 Newton 更新是

X1=(5/41/2−1/41),X2=(87/8815/44−27/8821/22).

奇异值从约 (2.28825,0.874032) 变为 (1.36263,1.00908),再变为 (1.04825,1.0000408);正交性残差 ‖XkTXk−I‖F 依次约为 4.24264,0.856957,0.0988337。

这个 A 已是正对角上三角矩阵,其常规 QR 分解可取 Q=I。但

‖A−I‖F=2,‖A−U‖F=8−210≈1.29439<2.

QR 给出的正交基解决列空间表示问题,极因子解决最近正交矩阵问题,两个目标产生不同答案。

正交因子可能是反射 ​

取 A=diag(−2,1)。其极因子是 U=diag(−1,1),行列式为 −1,距离为 1。如果应用要求行列式必须为 +1,可行集合就变成旋转群。

对二维旋转 R(θ),tr(ATR(θ))=−cos⁡θ,最大值在 θ=π 取得。故最近旋转是 −I,距离 5。直接返回普通极因子,会违反这个额外的方向保持约束。

秩亏时改变了哪些结论 ​

当 A 秩亏时,H=(ATA)1/2 仍唯一,但完整正交因子在零奇异值方向上不唯一,Newton 更新中的逆也不存在。矩形满列秩矩阵有相应的薄极分解,可先作 QR 再处理方形因子;本页迭代及唯一最近正交矩阵证明的主接口限定为非奇异方阵。

推论与应用

对于精确 Newton 轨道,第一轮以后所有奇异值至少为一,所以

‖Xk−U‖F2=∑i(σi−1)2≤14∑i(σi2−1)2=14‖XkTXk−I‖F2.

这给出轨道上正交性残差到因子误差的界。对任意一个外部给出的近正交矩阵,仅有小正交性残差当然不能证明它接近正确的 U;还要检查 H=UTA 的对称正定性与重构关系。

每轮需要一次多右端线性求解,稠密工作为 O(n3),存储 O(n2)。极分解可用于将数值演化后的矩阵拉回正交约束、刚体配准和正交近似,但必须先明确是否允许反射,以及接近秩亏时能接受怎样的敏感性。

单元任务:把函数值、导数与分支分别验清 ​

对 Tε=(12101+ε3004),取 ε=0.1,10−6,10−12,0,计算指数、主平方根和主对数。用 ε=0 的精确公式核对极限,再计算非交换方向 E=e2e1T 的指数导数。比较完整指数与只求 eTεe3,记录真正执行的乘法、求解和矩阵向量乘次数。

验收时需要分别给出函数值误差、平方或指数重构残差、主分支条件及导数交叉核验。另用四节点扩散例检查 Krylov 缺陷积分界,不能把终点缺陷当作误差。完整推导、数值和接受标准见单元题解;可运行核验脚本依赖 NumPy、SciPy 和 mpmath,并给出对应版本及八十位参照计算。

参考资料
  • Nicholas J. Higham, “Computing the Polar Decomposition—with Applications,” SIAM Journal on Scientific and Statistical Computing 7(4), 1986, pp. 1160–1174, §§2.1–2.2、3.1–3.2:极分解、最近正交因子及 Newton 奇异值迭代。作者版本
  • Nicholas J. Higham and Lijing Lin, Matrix Functions: A Short Course, 2013, §6.9:矩阵迭代与稳定实现的区分。作者版本
关系图谱7 个相邻概念 · 3 类关系

拖动节点调整位置。

显示关系

显示:依赖

  1. 前置三跳
  2. 前置二跳
  3. 前置一跳
  4. 当前条目
  5. 后续一跳
  6. 后续二跳
  7. 后续三跳
文字版关系按与当前条目的最短距离分组
分类位置

上位 / 更一般

下位 / 直接特例

暂未标注直接特例。

类型化关系

使用的工具