How to Raise a square matrix to a negative half power in R?

The singular power operation can be performed on square matrices itself, where the specified power is applied to each element of the matrix. Base R has many methods and routines to compute power for any k > 1 , where k is the integer value. However, raising matrices to a non-integral power is a challenge, and a limited number of solutions are available. External packages can be invoked to perform power computation in a more reliable manner in R Programming Language.

Method 1: Using sqrtm method

The expm package in R is used to compute exponential powers, logarithmic powers and square roots of matrices in R. The package needs to be first installed into the working space by the following command execution : 

install.packages("expm")

The package has the method sqrtm() which is used for matrix square root computation of the square matrices in R. 

sqrtm (mat)

This is followed by the application of the solve() method in R which, in general, solves the equation a %*% x = b for x, where b can be either a vector or a matrix or a real number value <0. This method is available in base R itself. 

solve( mat , power)

In case the power is empty, the inverse of the matrix mat is computed, that is the matrix is raised to the power of -1. Hence, to summarize, we compute the half-power using the sqrtm() method and negate it using the solve() method 

Code:

R

# loading the required libraries
library ("expm")

# declaring a data frame in R
mat <- matrix(1 : 9, ncol = 3)
print ("Original Matrix")
print (mat)
print ("Power matrix")

# computing -0.5 power of matrix
pow <- solve(sqrtm(mat))
print (pow)

Output:

[1] "Original Matrix" 
     [,1] [,2] [,3] 
[1,]    1    4    7 
[2,]    2    5    8 
[3,]    3    6    9 
[1] "Power Matrix" 
          [,1]        [,2]        [,3]
[1,] 0- 6980476i 0+13960950i 0- 6980475i 
[2,] 0+13960950i 0-27921901i 0+13960951i 
[3,] 0- 6980475i 0+13960951i 0- 6980475i

Method 2: Using eigen vectors approach

^%^ operator can be used to compute the matrix power operation, where eigen value decomposition is performed for the specified matrix. Its corresponding vectors, as well as values, are obtained in the form of arrays. The transpose of the vectors along with the values vector is then used to return the customized function value to compute matrix power. However, this method is considered unsuitable for the working of random square matrices, as it is guided by many constraints, some of them being : 

• This method doesn’t work for a matrix that has no eigen value decomposition.
• This method doesn’t work for a matrix which is not diagonalized.
• The matrix should be preferably symmetric.

Code: 

R

# creating data for matrix
vec <- c(0.088150041, 0.001017491 ,
         0.001017491, 0.084634294)

# declaring matrix
mat <- matrix(vec , nrow = 2)
print ("Original Matrix")
print (mat)

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

power <- -0.5

# raising matrix to -0.5 power
pow <- mat%^%(power)
print ("Modified Matrix")
print (pow)

Output:

[1] "Original Matrix" 
        [,1]     [,2] 
[1,] 0.088150 0.001017 
[2,] 0.001017 0.084634 
[1] "Modified Matrix" 
[,1]     [,2] 
[1,]  3.36830 -0.02004 
[2,] -0.02004  3.43755