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

**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

**Reference**

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

## No comments:

## Post a Comment