Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

# Introduction

“There is no unique definition of a square root of a matrix. While there are a few ways to calculate it in closed form, the results differ”. It was from these notes that I read in my multivariate statistics class that inspired me to explore the the geometric interpretations of varying results. In this blog, I am going to explore the geometric differences between the calculated square root of a symmetric positive-definite matrix via the singular value decomposition (SVD) and with the Cholesky decomposition.

# Some Math

## Using the SVD

Using the SVD of a matrix $A = UDV'$ , we can calculate the square root $A^{\frac{1}{2}} = UD^{\frac{1}{2}}V'$. Where $D$ is the diagonal matrix of the eigenvalues, and $U$ and $V$, for our context, are matrices containing eigenvectors calculated (i.e. $U=V$).

As a result, the code to find the $A^{\frac{1}{2}}$ using SVD is:

# Our symmetric matrix

A<- matrix(c(5,-1,-1,5),byrow=TRUE,ncol=2)

A

##      [,1] [,2]
## [1,]    5   -1
## [2,]   -1    5

# SVD of A

sVD<-svd(A)

sVD

## $d ##  6 4 ## ##$u
##            [,1]      [,2]
## [1,] -0.7071068 0.7071068
## [2,]  0.7071068 0.7071068
##
## $v ## [,1] [,2] ## [1,] -0.7071068 0.7071068 ## [2,] 0.7071068 0.7071068 # The square root of A using SVD svdResult<-sVD$u %*% diag(sqrt(sVD$d),ncol=2) %*% t(sVD$v)

svdResult

##            [,1]       [,2]
## [1,]  2.2247449 -0.2247449
## [2,] -0.2247449  2.2247449


## Using the Cholesky decompositon

The Cholesky decompositon of a real symmetric positive definite matrix is defined as $A = LL'$, where $L$ is a lower triangular matrix. As a result $A^{\frac{1}{2}} = L$.

The code for this is:

# Transposing to get the lower triangular matrix
choleskyResult<-t(chol(A))

choleskyResult

##            [,1]    [,2]
## [1,]  2.2360680 0.00000
## [2,] -0.4472136 2.19089


As observed, we note that the square-root matrices created by using the SVD differ from that using the Cholesky decomposition. Lets see how they differ geometrically.

# Visualizing the results.

Using the matlib package, we are able to visualize the above matrices. The code for this was inspired from this StackOverflow post, so check out the answers there for context.

library(matlib)

# Setting our limits
xlim <- c(-3,5)
ylim <- c(-2,5)

# Making a blank plot
plot(xlim, ylim, type="n", xlab="x", ylab="y", asp=1)
grid()

vectors(A[,1], labels="", pos.lab=4, frac.lab=.5, col="Red")
vectors(A[,2], labels="", pos.lab=4, frac.lab=.5,col="Red")
vectors(choleskyResult[,1], labels="", pos.lab=4, frac.lab=.5, col="black")
vectors(choleskyResult[,2], labels="", pos.lab=4, frac.lab=.5,col="black")
vectors(svdResult[,1], labels="", pos.lab=4, frac.lab=.5,col="blue")
vectors(svdResult[,2], labels="", pos.lab=4, frac.lab=.5,col="blue")
legend(x=-6,
y= 0,
legend=c("A Matrix",
expression(paste("A"^frac(1,2)," Using SVD")),
expression(paste("A"^frac(1,2),"Using Cholesky"))),
col=c("red", "black","blue"),
lty=1:2, cex=0.8,
box.lty=0)
title(main=expression(paste("Calculating ","A"^frac(1,2))))



The vectors generated for $A^{\frac{1}{2}}$ appear to have the same lengths but differ in rotation.

xlim <- c(-1,3)
ylim <- c(-1,3)

plot(xlim, ylim, type="n", xlab="x", ylab="y", asp=1)
grid()
vectors(choleskyResult[,1], labels="", pos.lab=4, frac.lab=.5, col="black")
vectors(choleskyResult[,2], labels="", pos.lab=4, frac.lab=.5,col="black")
vectors(svdResult[,1], labels="", pos.lab=4, frac.lab=.5,col="blue")
vectors(svdResult[,2], labels="", pos.lab=4, frac.lab=.5,col="blue")

legend(x=-3,
y= 0,
legend=c(expression(paste("A"^frac(1,2)," Using SVD")),
expression(paste("A"^frac(1,2),"Using Cholesky"))),
col=c("black","blue"),
lty=1:2, cex=0.8,
box.lty=0)
title(main=expression(paste("SVD vs Cholesky for Calculating ","A"^frac(1,2))))



Interestingly enough, they both share the same distances between their rows, which seems to me that the calculation of the square root of a symmetric matrix using the SVD vs the Cholesky decomposition only differs in rotation.

dist(choleskyResult)

##          1
## 2 3.464102

dist(svdResult)

##          1
## 2 3.464102


# Conclusion

I really have only began to scratch the surface of this, but these are my observations along the way. Do you know anything else about this? Let me know!

## Want to see more of my content?

Be sure to subscribe and never miss an update!