主矩阵平方根与 Denman–Beavers 迭代

形式陈述

满足 $S^2=A$ 的矩阵可能很多。若 $A$ 的谱避开闭负实轴，标量主平方根通过矩阵函数演算给出唯一的主平方根 $S=A^{1/2}$ ，其谱位于开右半平面。Denman–Beavers 迭代同时逼近它和它的逆：

Y_0=A,\qquad Z_0=I,

Y_{k+1}=\frac12(Y_k+Z_k^{-1}),\qquad Z_{k+1}=\frac12(Z_k+Y_k^{-1}).

两次更新必须都读取旧的 $Y_k,Z_k$ 。 实现时分别解 $Z_kU_k=I$ 、 $Y_kV_k=I$ ，再同时赋值 $Y_{k+1}=(Y_k+U_k)/2$ 、 $Z_{k+1}=(Z_k+V_k)/2$ 。在上述谱条件下，精确运算中的迭代始终有定义，并满足 $Y_k\to S$ 、 $Z_k\to S^{-1}$ 。

为什么这对迭代收敛到主根

从初始化出发，所有迭代矩阵都是 $A$ 的有理函数，彼此交换，并保持 $Y_k=AZ_k$ 。因此第一条更新也可写成

Y_{k+1}=\frac12(Y_k+AY_k^{-1}).

成对状态 $(Y_k,Z_k)$ 每次应用同一个更新映射；上述谱条件保证初始化轨道始终留在有定义的状态内，因此这是不动点迭代的一种具体形式。为了证明它选择正确的根，定义 Cayley 型误差

W_k=(Y_k-S)(Y_k+S)^{-1}.

利用各矩阵交换和 $S^2=A$ ，直接代数化简得到

W_{k+1}=W_k^2,\qquad W_0=(S-I)(S+I)^{-1}.

若 $\operatorname{Re}\lambda(S)>0$ ，则 $|\frac{\lambda-1}{\lambda+1}|<1$ 。故 $\rho(W_0)<1$ ， $W_k=W_0^{2^k}\to0$ ，并由

Y_k=S(I+W_k)(I-W_k)^{-1}

得到 $Y_k\to S$ 。这些式子还说明各次所需逆矩阵存在。进入根的邻域后，误差满足二次上界；某些幂零结构甚至会有限步精确到达，因而不必具有恰为二的渐近阶。这个论证允许非对角化矩阵，因为谱半径小于一仍保证矩阵幂趋于零。

直觉

$Y$ 猜测平方根， $Z$ 猜测逆平方根。若一个方向被 $Y$ 放大过头，来自 $Z^{-1}$ 的信息会把它拉回；两条轨道在精确运算中始终保持同一个矩阵函数关系。

图中两条轨道同步前进。横向读数看平方残差是否下降，纵向关系看 $YZ$ 是否接近单位矩阵。两个量检查的是不同关系，都值得保留。

例子与边界

非对角矩阵的前两轮

取

A=\begin{pmatrix}4&6\\0&9\end{pmatrix},\qquad S=\begin{pmatrix}2&6/5\\0&3\end{pmatrix}.

$S^2$ 的超对角是 $2\cdot(6/5)+(6/5)\cdot3=6$ 。第一轮为

Y_1=\frac{A+I}{2}=\begin{pmatrix}5/2&3\\0&5\end{pmatrix},\qquad Z_1=\frac{I+A^{-1}}2=\begin{pmatrix}5/8&-1/12\\0&5/9\end{pmatrix}.

第二轮读取这两个旧矩阵，得到

Y_2=\begin{pmatrix}41/20&81/50\\0&17/5\end{pmatrix},\qquad Z_2=\begin{pmatrix}41/80&-97/600\\0&17/45\end{pmatrix}.

例如 $AZ_2$ 的超对角为 $4(-97/600)+6(17/45)=81/50$ ，再次验证 $Y_2=AZ_2$ 。两次残差的 Frobenius 范数如下：

$k$	$\|Y_k^2-A\|_F$	$\|Y_kZ_k-I\|_F$
0	$102.5280$	$10.4403$
1	$23.0936$	$2.3672$
2	$3.82071$	$0.402739$

继续迭代会迅速进入二次收敛阶段。脚本记录每一轮矩阵与残差，而不是只输出最终近似。

零残差不能识别分支

$-S$ 也满足 $(-S)^2=A$ ，平方残差完全为零；但其特征值为 $-2,-3$ ，不在开右半平面，因此不是主平方根。若同时令 $Z=-S^{-1}$ ，甚至 $YZ=I$ 也成立。必须额外检查主分支条件，不能靠这两个残差自动选根。

若 $A$ 的谱落在闭负实轴，上面的主分支和收敛定理不再适用。以 $A=-I$ 为例，第一轮 $Y_1=0$ ，下一轮求逆就失败；它虽在复数上有平方根，却不满足这里的算法假设。

精确收敛不是浮点稳定性的替身

在精确运算中，成对迭代可约化成单矩阵 Newton 形式；浮点下各矩阵不再严格交换，两种实现的误差传播可以不同。直接使用 $Y\leftarrow(Y+Y^{-1}A)/2$ 不能仅凭标量 Newton 的收敛性就宣称稳定。成对或乘积形式及缩放能改善表现，但求解病态、溢出和原问题敏感性仍需单独处理。

这正是前向与后向误差分析要补充的内容：残差衡量方程满足程度，谱分支选择目标，而条件性衡量残差能否可靠代表根的误差。

推论与应用

若近似根写成 $Y=S+E$ ，则

Y^2-A=SE+ES+E^2.

一阶误差由 Sylvester 算子 $E\mapsto SE+ES$ 控制。局部误差尺度约为残差除以 $\operatorname{sep}(S,-S)$ ；当该值很小时，平方残差小仍不足以保证根准确。这个关系也给出针对残差的 Newton 校正方程。

每轮成对更新需要两次 $n$ 右端的线性求解，稠密成本 $O(n^3)$ ，存储 $O(n^2)$ ；总成本还乘以实际迭代轮数。主平方根既是独立计算目标，也是逆缩放平方对数算法的关键步骤。

参考资料

Nicholas J. Higham and Lijing Lin, Matrix Functions: A Short Course, 2013, §6.9：平方根 Newton 迭代的数值风险、成对及乘积形式。作者版本
Nicholas J. Higham, “Functions of Matrices,” Handbook of Linear Algebra, 2nd ed., 2014, §17.9：主平方根与 Denman–Beavers 迭代。作者版本

形式陈述 ​

为什么这对迭代收敛到主根 ​

非对角矩阵的前两轮 ​

零残差不能识别分支 ​

精确收敛不是浮点稳定性的替身 ​

反向链接

形式陈述

为什么这对迭代收敛到主根

非对角矩阵的前两轮

零残差不能识别分支

精确收敛不是浮点稳定性的替身