R语言中如何正确计算对角矩阵的-1/2次幂？求解异常问题

我通过以下代码生成对角矩阵V：

X <- rnorm(n)
V <- matrix(diag(abs(X)), ncol = n)

需要计算V的-1/2次幂，尝试使用expm包的%^%运算符：

K <- V %^% (-1/2)

但得到的结果是全1的对角矩阵（如下方输出），而正确结果应该是每个对角元素v_i的(-1/2)次幂。

输出示例：

# 原对角矩阵V
[1,] 0.08378436 0.0000000 0.000000 0.0000000 0.0000000
[2,] 0.00000000 0.9829437 0.000000 0.0000000 0.0000000
[3,] 0.00000000 0.0000000 1.875067 0.0000000 0.0000000
[4,] 0.00000000 0.0000000 0.000000 0.1861447 0.0000000
[5,] 0.00000000 0.0000000 0.000000 0.0000000 0.6334857

# 错误的计算结果K
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    0    1    0    0    0
[3,]    0    0    1    0    0
[4,]    0    0    0    1    0
[5,]    0    0    0    0    1

对角矩阵的幂运算无需借助expm包的通用矩阵幂运算符，直接操作对角元素即可得到正确结果：

步骤1：提取对角矩阵的对角元素 #

用diag()函数提取V的对角元素：

diag_elements <- diag(V)

步骤2：计算每个对角元素的-1/2次幂 #

直接对提取出的元素做幂运算：

diag_powered <- diag_elements^(-1/2)

步骤3：重新构造对角矩阵 #

将计算后的元素放回对角矩阵：

K <- diag(diag_powered)

完整示例代码 #

set.seed(123) # 固定随机种子，确保结果可复现
n <- 5
X <- rnorm(n)
V <- diag(abs(X)) # 直接用diag()生成对角矩阵，原代码的matrix包裹是冗余操作

# 正确计算-1/2次幂
diag_elements <- diag(V)
diag_powered <- diag_elements^(-1/2)
K <- diag(diag_powered)

# 查看结果
print("原对角矩阵V：")
print(V)
print("正确的-1/2次幂矩阵K：")
print(K)

错误原因说明 #

expm包的%^%运算符是针对通用矩阵的指数运算，通过矩阵的特征值分解等方式计算。对于对角矩阵，当指数为分数时，可能因数值精度或计算逻辑的特殊处理，导致错误地将特征值的幂计算为1（比如极小值的浮点误差干扰）。而直接操作对角元素是对角矩阵幂运算的本质逻辑，既高效又准确。

内容的提问来源于stack exchange，提问作者Cole Hendrickson