2×2矩阵指数公式的鲁棒计算方法咨询 #

首先先明确你提到的Al-Mohy和Higham 2009年论文中的2×2矩阵指数公式：

对于矩阵
$$A = \begin{bmatrix}
a & b \
c & d
\end{bmatrix}$$
其矩阵指数为：
$$e^A = e^{(a+d) / 2}\begin{bmatrix}
\cosh(\mu / 2) + \frac{1}{2}(a-d)\text{sinch}(\mu/2) & b \cdot \text{sinch}(\mu/2) \
c \cdot \text{sinch}(\mu/2) & \cosh(\mu / 2) + \frac{1}{2}(a-d)\text{sinch}(\mu/2)
\end{bmatrix}$$
其中 $\mu = \sqrt{(a-d)^2 + 4bc}$

针对你遇到的$\mu$过大导致$\cosh$溢出的问题，结合数值计算的鲁棒性要求，这里给你几个针对性的解决方案，不用依赖高次矩阵幂：

1. 大$\mu$场景下的公式渐近变形 #

当$|\mu/2|$足够大时（比如大于10，此时$\cosh(\mu/2)$已经会快速趋近于$e^{\mu/2}/2$），我们可以用双曲函数的渐近近似来改写原公式，避免直接计算大参数的$\cosh$：

• $\cosh(x) \approx \frac{e^x}{2}$，$\sinh(x) \approx \frac{e^x}{2}$
• $\text{sinch}(x) = \frac{\sinh(x)}{x} \approx \frac{e^x}{2x}$

把这些近似代入原式，提取公共的指数因子后，得到：
$$e^A \approx \frac{e^{(a+d+\mu)/2}}{2}\begin{bmatrix}
1 + \frac{a-d}{\mu} & \frac{2b}{\mu} \
\frac{2c}{\mu} & 1 + \frac{a-d}{\mu}
\end{bmatrix}$$
这个形式完全避开了大参数的双曲函数计算，直接通过指数的线性组合来表达，不会触发溢出。你可以设置一个阈值（比如$|\mu|>20$），当$\mu$超过阈值时用这个近似，否则用原公式。

2. 改进版缩放平方策略，避免高次幂 #

你之前用Frobenius范数缩放的思路方向是对的，但不用直接计算矩阵的$|A|$次幂，而是用二进制分解的平方操作——这也是Al-Mohy和Higham论文的核心思路：

• 选择一个整数$k$，使得$|A/2^{k|$足够小（比如小于1，这样计算$e}{A/2^{k}$时，$\mu/2}{k+1}$会很小，双曲函数不会溢出）
• 计算$B = e^{A/2k}$（此时用原公式完全安全）
• 对$B$重复平方$k$次，得到$e^A = B^{2k}$

这里的平方操作只是矩阵自身相乘，效率非常高，而且$k$的取值通常很小（比如当$|A|=1000$时，$k=10$就足够，因为$1000/2^{10}≈0.976$），完全不会有高次幂的计算负担。

3. 特征值分解辅助计算（分场景使用） #

如果你的矩阵$A$是可对角化的，也可以用特征值分解来计算：

• $A$的特征值为$\lambda_1 = \frac{(a+d)+\mu}{2}$，$\lambda_2 = \frac{(a+d)-\mu}{2}$
• 找到对应的特征向量组成可逆矩阵$P$，那么$e^A = P \cdot \text{diag}(e^{\lambda_1}, e^{\lambda_2}) \cdot P^{-1}$

不过要注意，当$\mu$非常小时（矩阵接近不可对角化），这个方法的数值稳定性会下降，此时原公式的$\text{sinch}$项反而更鲁棒，所以建议只在$\mu$较大时结合这个方法使用。

总结最优流程 #

你可以把这些方法结合起来，形成一个鲁棒的计算流程：

• 计算矩阵的迹$t = a+d$和$\mu = \sqrt{(a-d)^2 + 4bc}$
• 判断$|\mu|$的大小：
  • 若$|\mu|$过大（超过阈值）：用渐近近似公式计算
  • 若$|\mu|$适中：直接用原公式计算
  • 若矩阵的范数过大但$\mu$适中：用缩放平方策略，先计算小矩阵的指数再重复平方

这样既解决了溢出问题，又保证了数值稳定性，还不用做复杂的高次矩阵幂操作。

备注：内容来源于stack exchange，提问作者Ghanem