如何计算R中矩阵的幂

Question

我正在尝试计算以下矩阵的-0.5幂:

S <- matrix(c(0.088150041, 0.001017491 , 0.001017491, 0.084634294),nrow=2)

在Matlab中,结果是(S^(-0.5)):

S^(-0.5)
ans =
 3.3683   -0.0200
-0.0200    3.4376

Answer 1

> library(expm)
> solve(sqrtm(S))
            [,1]        [,2]
[1,]  3.36830328 -0.02004191
[2,] -0.02004191  3.43755429

Answer 2

一段时间后，出现了以下解决方案：

"%^%" <- function(S, power) 
   with(eigen(S), vectors %*% (values^power * t(vectors))) 
S%^%(-0.5)

结果给出了预期的答案：

              [,1]        [,2]
  [1,]  3.36830328 -0.02004191
  [2,] -0.02004191  3.43755430

它们只是“本征”的结果。 (2认同)
警告！这是错误的，并且在一般情况下不起作用。根据线性代数课程，此解决方案适用于将 `t(vectors)` 替换为 `solve(vectors)`，但前提是 S 可对角化。它也适用于其他一些情况，例如，如果 S 是对称的，因为矩阵 `vectors` 是正交的，即 `t(vectors) = solve(vectors)`（[参见此处]（https://en.wikipedia） .org/wiki/Diagonalizable_matrix#Diagonalizable_matrices))。这就是为什么它适用于原始问题给出的矩阵，这是对称的。但是再一次，对于随机矩阵，这失败了。见上一条评论。 (2认同)

Answer 3

矩阵的平方根不一定是唯一的（大多数实数至少有 2 个平方根，因此不仅仅是矩阵）。有多种算法可用于生成矩阵的平方根。其他人已经展示了使用 expm 和特征值的方法，但 Cholesky 分解是另一种可能性（请参阅chol函数）。