Exponentiation of a matrix (including pseudoinverse)

The following function "exp.mat" allows for the exponentiation of a matrix (i.e. calculation of a matrix to a given power). The function follows three steps:
1) Singular Value Decomposition (SVD) of the matrix
2) Exponentiation of the singular values
3) Re-calculation of the matrix with the new singular values

The most common case where the method is applied is in the calculation of a "Moore–Penrose pseudoinverse", or a matrix raised to the power of -1. In this case, the function returns the same result as the "pseudoinverse" function from the corpcor package (although maybe not as nicely implemented). This should also be analogous to the pinv function in Matlab.

Less common is the need to calculate other powers of matrices. For example, I use this function to calculate the square root of a matrix (and square root of the inverse of a matrix) within Canonical Correlation Analysis (CCA). I hope to write a post shortly on the use of CCA in identifying coupled patterns in climate data.

Below is the code for the exp.mat function as well as some demonstrations of its use in R.


the exp.mat function:

#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)
}

Created by Pretty R at inside-R.org



some demonstrations:

#required functions (from "www.menugget.blogspot.com")
source("exp.mat.R")
 
#example matrix from Wilks (2006)
A <- matrix(c(185.47,110.84,110.84,77.58),2,2)
A
solve(A) #inverse
exp.mat(A, -1) # pseudoinverse
exp.mat(exp.mat(A, -1), -1) #inverse of the inverse -return to original A matrix
exp.mat(A, 0.5) # square root of a matrix
exp.mat(A, -0.5) # square root of its inverse
exp.mat(exp.mat(A, -1), 0.5) # square root of its inverse (same as above)
 
#Compare to pseudoinverse function or corpcor package
library("corpcor")
solve(A)
pseudoinverse(A)
exp.mat(A, -1)
 
#Pseudoinversion of a non-square matrix
set.seed(1)
D <- matrix(round(runif(24, min=1, max=100)), 4, 6)
D
pseudoinverse(D)
exp.mat(D, -1)
pseudoinverse(t(D))
exp.mat(t(D), -1)
 
#Pseudoinversion of a square matrix
set.seed(1)
D <- matrix(round(runif(25, min=1, max=100)), 5, 5)
solve(D)
pseudoinverse(D)
exp.mat(D, -1)
solve(t(D))
pseudoinverse(t(D))
exp.mat(t(D), -1)
 
 
###Examples from corpcor manual
# a singular matrix
m = rbind(
c(1,2),
c(1,2)
)
 
# not possible to invert exactly
solve(m)
pseudoinverse(m)
exp.mat(m, -1) #the exp.mat comparison
p <- pseudoinverse(m)
 
# characteristics of the pseudoinverse
zapsmall( m %*% p %*% m )  ==  zapsmall( m )
zapsmall( p %*% m %*% p )  ==  zapsmall( p )
zapsmall( p %*% m )  ==  zapsmall( t(p %*% m ) )
zapsmall( m %*% p )  ==  zapsmall( t(m %*% p ) )
 
# example with an invertable matrix
m2 = rbind(
c(1,1),
c(1,0)
)
zapsmall( solve(m2) ) == zapsmall( pseudoinverse(m2) )
zapsmall( solve(m2) ) == zapsmall( exp.mat(m2,-1) ) #the exp.mat comparison

Created by Pretty R at inside-R.org



Reference
Wilks, D.S., 2006. Statistical Methods in the Atmospheric Sciences. Second Edition. Academic Press.