PmWiki | RExcursions / SVD

First let us create a sample dataset with which to work.

  
 rect = vector("list",2)
 rect[[1]] = 1:40
 rect[[1]] = c(rect[[1]],rect[[1]])
 rect[[1]] = c(rect[[1]],rect[[1]])

 for(i in 1:4){
   for(j in 1:40){
     rect[[2]][((i-1)*40)+j] = i
   }
 }

 plot(rect[[1]],rect[[2]],xlim = c(0,40), ylim=c(-20,20))

What we should now have is a set of points whose X varies from 1:40 and whose Y varies from 1:4. A 4x40 rectangle of points. We can simplify the datastructure a bit by turning it into a matrix.

 rectMat = matrix(0,160,2)
 rectMat[,1] = rect[[1]]
 rectMat[,2] = rect[[2]]

Or alternatively

 rectMat = matrix(c(rect[[1]],rect[[2]]),160,2)

We can center this point space on the coordinate space by creating a simple function

 center <- function(matrix){
   ans = matrix
   ans[,1] = matrix[,1] - mean(matrix[,1])
   ans[,2] = matrix[,2] - mean(matrix[,2])

  ans
 }

We can now center on the fly as in this plot statement. A more generic n dimensional center function could easily be constructed with a loop. A pre-written function probably also exists, though it is a tossup whether or not finding it would be faster than implementing it for simple problems. In particular, if you need to fully understand what is happening, it is often useful to implement one's self. If you want to ensure no errors exist (or at least hopefully so) it is better to search and find existing functions. How one performs operations in R depends largely on the goal of the R interaction.

 plot(center(rectMat))

 theta = pi/4
 rotationMat = matrix(c(cos(theta),sin(theta),-sin(theta),cos(theta)),2,2)
 plot(center(rectMat) %*% rotationMat)

Since this is what we will be working with for now we can simply assign these values and eliminate needing to constantly use functions.

 rectMat = center(rectMat)
 rectRot = rectMat %*% rotationMat

Now that we have this point space, the interesting bit can begin. We start by running a Singular Value Decomposition on our space.

 svdAns = svd(rectRot)

The singular value decomposition has 3 parts. In R the functions defines this as X = UDV' though the notation varies somewhat depending on the source. If we examine this particular example, we can see that V is the rotation applied to X to put it into the orthogonal axis of our coordinate system. V' is precisely the rotation we applied originally to rotate the original rectangle.

> svdAns$v
           [,1]      [,2]
[1,]  0.7071068 0.7071068
[2,] -0.7071068 0.7071068
> rotationMat
          [,1]       [,2]
[1,] 0.7071068 -0.7071068
[2,] 0.7071068  0.7071068
> t(svdAns$v)
          [,1]       [,2]
[1,] 0.7071068 -0.7071068
[2,] 0.7071068  0.7071068

D is a vector which represents the diagonal entries of a diagonal matrix containing the Singular Values. X and D must have the same rank so at most p (where our original matrix was nxp) non-zero values on the diagonal. D can be used to interpret the amount of variance which exists in orthogonal directions. In our case we are working with 2 dimensional points (p = 2) and variance exists in both directions. If D had a rank of 1 it would indicate that, despite p = 2 in our original data, the data can be transformed to have all of it's variance lie along a single dimension. SVD should recover this as it tries to decompose the original space into UDV' where V and U are orthogonal.

If we look at D in this example

 > svdAns$d
 [1] 146.01370  14.14214

We can use this to recover our original un-rotated data.

 > sum(rectMat - (svdAns$u %*% diag(svdAns$d)))
 [1] 1.154632e-14

The answer is very near 0. Since SVD is being calculated using an iterative solver, the answer doesn't quite come to zero, however if you look at the matrix resulting from this subtraction and you round each to 0, you find that the two matrices differ only very slightly as a result of this long floating-point result from the solver. Since D is used to scale U to the full variance of the original space, D can be seen as a relative difference in the amount of variance in the orthogonal dimensions present in the data. If both entries in D were identical, the original data is equally spread in both of it's dimensions. If there is a large difference between the two, as is in this example, the data in much more spread-out in one dimension than the other.