如何在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

- WAF

Answer 1

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

- caracal

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

- WAF

---

你应该解释一下 vectors 和 values 来自哪里，否则几周（或几小时）后没人会记得 :-0 - Carl Witthoft

它们只是“eigen”的结果。 - Sacha Epskamp

这对于没有特征值分解的矩阵有效吗？https://en.wikipedia.org/wiki/Diagonalizable_matrix#Matrices_that_are_not_diagonalizable? - Basj

等等，这个 %^% 连平方都不行：使用 A = matrix(c(0,1,2,2),2,2)，A %*% A 的结果是正确的，但是 A %^% 2 的结果却不同。 - Basj

警告！这是错误的，通常不起作用。根据线性代数课程，此解决方案将使用t(vectors)替换为solve(vectors)才能正常工作，但仅当S可对角化时才有效。它在某些其他情况下也可以工作，例如如果S是对称的，因为那么矩阵vectors是正交的，即t(vectors) = solve(vectors)（见此处）。这就是为什么它适用于原始问题给出的矩阵，该矩阵是对称的。但再次强调，对于随机矩阵，这种方法会失败。请参见上一个评论。 - Basj

Answer 3

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

- Greg Snow

---

确实。我找到了这个（显然的）页面，http://en.wikipedia.org/wiki/Square_root_of_a_matrix，并惊讶/恐惧地发现2x2单位矩阵有无限多个由毕达哥拉斯三元数组成的平方根！ - Carl Witthoft

Answer 4

为了将此答案扩展到平方根以外，以下函数exp.mat()推广了矩阵的“Moore–Penrose pseudoinverse”，并允许通过奇异值分解（SVD）来计算矩阵的指数（即使对于非正方形矩阵也有效，尽管我不知道什么时候需要这样做）。

#The exp.mat function performs can calculate the pseudoinverse of a matrix (EXP=-1)
#and other exponents of matrices, such as square roots (EXP=0.5) or square root of 
#its inverse (EXP=-0.5). 
#The function arguments are a matrix (MAT), an exponent (EXP), and a tolerance
#level for non-zero singular values.
exp.mat<-function(MAT, EXP, tol=NULL){
 MAT <- as.matrix(MAT)
 matdim <- dim(MAT)
 if(is.null(tol)){
  tol=min(1e-7, .Machine$double.eps*max(matdim)*max(MAT))
 }
 if(matdim[1]>=matdim[2]){ 
  svd1 <- svd(MAT)
  keep <- which(svd1$d > tol)
  res <- t(svd1$u[,keep]%*%diag(svd1$d[keep]^EXP, nrow=length(keep))%*%t(svd1$v[,keep]))
 }
 if(matdim[1]<matdim[2]){ 
  svd1 <- svd(t(MAT))
  keep <- which(svd1$d > tol)
  res <- svd1$u[,keep]%*%diag(svd1$d[keep]^EXP, nrow=length(keep))%*%t(svd1$v[,keep])
 }
 return(res)
}

示例

S <- matrix(c(0.088150041, 0.001017491 , 0.001017491, 0.084634294),nrow=2)
exp.mat(S, -0.5)
#            [,1]        [,2]
#[1,]  3.36830328 -0.02004191
#[2,] -0.02004191  3.43755429

其他示例可以在这里找到。

- Marc in the box