RObservations #14: Comparing the Calculated Square Roots of Symmetric Postive Matrices

[This article was first published on r – bensstats, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
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
## [1] 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()

# Adding our vectors
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!

To leave a comment for the author, please follow the link and comment on their blog: r – bensstats.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)