形式陈述
矩阵方程 e X = A 可能有许多解。怎样选出一个可连续计算的对数,并说明何时不能这样选?
设 A ∈ C n × n 的谱不与闭负实轴 ( − ∞ , 0 ] 相交。把标量主值对数作用为 primary 矩阵函数 公理库 主矩阵函数与 Hermite 插值 Primary matrix function · Matrix function by Hermite interpolation 用极小多项式规定所需的函数值与导数,经 Hermite 插值定义一般方阵的函数,并区分 primary 与 principal 两种不同的分支要求。 ,得到 主矩阵对数 Log A 。它满足
e Log A = A , σ ( Log A ) ⊂ { z : − π < Im z < π } . 它是谱位于这个开条带内的唯一矩阵对数。如果 A 为实矩阵,Log A 也为实矩阵。
存在性来自标量主值的解析函数演算。若另一矩阵 X 的谱在条带内且 e X = A ,标量恒等式 Log ( e z ) = z 在该条带解析成立,代入矩阵得 Log ( e X ) = X ,故 X = Log A 。实矩阵结论则由复共轭与唯一性得到。
逆缩放平方的计算路线
先取若干次主平方根
B = A 1 / 2 s , X = B − I . 主分支保证 Log A = 2 s Log B 。当 X 足够小时,使用积分表示
Log ( I + X ) = ∫ 0 1 X ( I + t X ) − 1 d t . 以 m ≥ 1 点 Gauss–Legendre 求积 公理库 Gaussian 求积 Gaussian quadrature · Gauss quadrature 以正交多项式零点选择节点,使 n 个正权节点对最高 2n−1 次多项式精确积分。 的节点 t j 、权重 w j 近似,得到一个 Padé 有理函数 公理库 有理函数 Rational function · 有理式 由两个多项式之商表示,并按交叉相乘关系识别相等表示的函数域元素。 :
Log A ≈ 2 s ∑ j = 1 m w j X ( I + t j X ) − 1 . 每项通过多右端线性求解计算。实践中先作 Schur 分解 公理库 Schur 分解 Schur decomposition · Schur triangularization · Real Schur form 用酉或正交相似变换把一般矩阵化为上三角或实准上三角形,作为浮点全谱计算的稳定结构目标。 ,对三角矩阵重复开方和求解,再映回原坐标,能利用三角结构。
一个明确但保守的截断界
若在相容范数下 r = ‖ X ‖ < 1 ,展开积分中的 resolvent。m 点 Gauss–Legendre 对次数不超过 2 m − 1 的多项式精确,所以前 2 m 项矩阵幂的系数一致。利用权重为正且总和为一,可得
‖ Log A − 2 s ∑ j = 1 m w j X ( I + t j X ) − 1 ‖ ≤ 2 s 2 r 2 m + 1 1 − r . 这是精确平方根、精确求解前提下的逼近截断界,不包含浮点开方误差。高效实现通常用更紧的幂范数后向误差界选择 s , m 。
直觉
指数算法先缩小输入再反复平方;对数算法先对输入反复开方,使它靠近单位矩阵,最后把一个小对数放大 2 s 倍。每次必须沿同一主分支前进,否则“开方后再乘二”可能跳到另一个对数。
图片加载失败 图的上部用复数 e i θ 标记旋转参数,并标出对应的主矩阵对数;负实轴是分支边界。中部展示正谱矩阵反复开方,特征值逐渐靠近 1 。靠近单位阵便于逼近,也可能使 B − I 的直接相减丢失精度,缩放次数需要有节制。
例子与边界
旋转跨过 π 时会发生什么
记 J = ( 0 − 1 1 0 ) ,R ( θ ) = e θ J 。当 − π < θ < π 时,Log R ( θ ) = θ J 。当 π < θ < 3 π 时,主值角要减去 2 π ,所以主对数变成 ( θ − 2 π ) J 。
令 θ 从两侧趋于 π ,两个极限分别是 π J 和 − π J 。在 θ = π ,矩阵等于 − I ,其谱落在切线上,主对数不定义;但它确实有实对数 π J 。“没有主对数”不等于“没有任何对数”。
类似地,e 2 π J = I ,而 Log I = 0 。因此 Log ( e X ) = X 必须附带谱位于主值条带的条件。
两次开方,再做三点有理逼近
取
A = ( 4 6 0 9 ) , Log A = ( log 4 6 5 log ( 9 / 4 ) 0 log 9 ) . 第一次主平方根是 B 1 = ( 2 6 / 5 0 3 ) 。第二次为
B 2 = ( 2 6 / 5 2 + 3 0 3 ) ≈ ( 1.41421356 0.38140469 0 1.73205081 ) . 令 X = B 2 − I 。三点 Gauss–Legendre 节点为 ( 1 − 3 / 5 ) / 2 , 1 / 2 , ( 1 + 3 / 5 ) / 2 ,对应权重 5 / 18 , 4 / 9 , 5 / 18 。三次求解后乘以 4 ,得到
L ≈ ( 1.38629352 0.97309258 0 2.19720400 ) . 它与精确对数的 1 -范数误差约为 4.42 × 10 − 5 。这里 ‖ B 2 − I ‖ 1 ≈ 1.11346 > 1 ,不能套用上面的范数截断界;例子只是执行示范。
继续到四次开方,得到 r ≈ 0.21523664 。使用十点求积,保守截断界降至 4.00 × 10 − 13 以下;附带脚本还与八十位精度的主对数交叉核验。这个组合的求值核心是四次三角平方根与十次多右端求解。
接近单位阵时要防止相减损失
过多开方会让 B 极接近 I ,而开方阶段的小误差在 2 s ( B − I ) 中被放大。标量中可用 a − 1 = ( a − 1 ) / ( a + 1 ) 减少相消;成熟三角算法还会从原始数据重算对角和第一超对角的关键条目。上面的教学核展示流程和截断界,没有实现全部这些精度修正。
推论与应用
三角平方根由 R 2 = T 的递推得到:r i i = t i i ,按超对角推进时
( r i i + r j j ) r i j = t i j − ∑ i < k < j r i k r k j . 主根保证分母不为零。块版本对应 Sylvester 方程 公理库 Sylvester 方程与谱分离 Sylvester equation 刻画 AX−XB=C 对所有右端唯一可解的充要条件,并用非正规矩阵说明特征值间距不能代替 sep 敏感性。 ,也说明分支条件为何进入算法本身。稠密 Schur 路线的总成本约为 O ( ( 1 + s + m ) n 3 ) ,存储 O ( n 2 ) 。
对数常用于从一步传播矩阵恢复连续生成元,但多值性意味着必须明确所选分支。验算 e L ≈ A 只能验证“某个对数”;要接受主对数,还应检查 L 的谱虚部在 ( − π , π ) ,并注意输入谱接近切线时的敏感性。
参考资料
Nicholas J. Higham and Lijing Lin, Matrix Functions: A Short Course , 2013, §3.1.6:主对数的存在、谱条带及实矩阵性质。作者版本
Awad H. Al-Mohy and Nicholas J. Higham, “Improved Inverse Scaling and Squaring Algorithms for the Matrix Logarithm,” SIAM Journal on Scientific Computing 34(4), 2012, pp. C153–C169, §§1–4:Gauss–Legendre 部分分式、后向误差、相消修正和 Schur 算法。作者版本