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 , we can calculate the square root . Where is the diagonal matrix of the eigenvalues, and and , for our context, are matrices containing eigenvectors calculated (i.e. ).

As a result, the code to find the 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 , where is a lower triangular matrix. As a result .

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