## Tuesday, July 14, 2009

### Matrix manipulation

Here is a little more "homework" related to matrices. Suppose we construct a simple 3 x 3 matrix using R.

 `v = c(2,1,4,2,-3,1,1,3,2)m = matrix(v,nrow=3,byrow=T)rownames(m) = c('a','b','c')colnames(m) = rownames(m)`

 `> m a b ca 2 1 4b 2 -3 1c 1 3 2`

The identity matrix is I:

 `I = diag(3)rownames(I) = rownames(m)colnames(I) = rownames(m)`

 `> I a b ca 1 0 0b 0 1 0c 0 0 1`

If we switch rows 1 and 3 in I and then do I x m, we switch the corresponding rows in m:

 `Tx = Itemp = Tx[1,]Tx[1,] = Tx[3,]Tx[3,] = temp`

 `> Tx a b ca 0 0 1b 0 1 0c 1 0 0> Tx %*% m a b ca 1 3 2b 2 -3 1c 2 1 4`

Our real problem today is to find the inverse of M, such that M-1 x M = I. It turns out that if we do a series of transformations of M that end by converting it to I, and then do the same series of transformations on I, we generate M-1!

If T x M = I, then T x I = M-1.

To make things easier to follow, we'll bind our matrix and the identity matrix together:

 `M = cbind(m,I)colnames(M) = c('a','b','c','d','e','f')`

 `> M a b c d e fa 2 1 4 1 0 0b 2 -3 1 0 1 0c 1 3 2 0 0 1`

The first thing we'll do is add -2 x row 3 to both row 1 and row 2. At the same time, we switch row 3 and row 1.

 `T1 = matrix(c(0,0,1,0,1,-2,1,0,-2), nrow=3,byrow=T)M1 = T1 %*% M`

 `> T1 [,1] [,2] [,3][1,] 0 0 1[2,] 0 1 -2[3,] 1 0 -2> M1 a b c d e f[1,] 1 3 2 0 0 1[2,] 0 -9 -3 0 1 -2[3,] 0 -5 0 1 0 -2`

At this point, if we had a system of simultaneous equations, we would have solved the system. The third row allows us to solve for y, plugging y into the second row and we can solve for z, and then into the first and we solve for x. But here, we continue to transform until we get the identity matrix I.

Next, we multiply row 2 by -1/3:

 `T2 = IT2[2,2] = -1/3M2 = T2 %*% M1`

 `> M2 a b c d e fa 1 3 2 0 0.0000000 1.0000000b 0 3 1 0 -0.3333333 0.6666667c 0 -5 0 1 0.0000000 -2.0000000`

We add -1 x row 2 to row 1:

 `T3 = IT3[1,2] = -1M3 = T3 %*% M2`

 `> M3 a b c d e fa 1 0 1 0 0.3333333 0.3333333b 0 3 1 0 -0.3333333 0.6666667c 0 -5 0 1 0.0000000 -2.0000000`

We multiply row 2 by 5 and add it to three times row 3:

 `T4 = IT4[3,2] = 5T4[3,3] = 3M4 = T4 %*% M3`

 `> M4 a b c d e fa 1 0 1 0 0.3333333 0.3333333b 0 3 1 0 -0.3333333 0.6666667c 0 0 5 3 -1.6666667 -2.6666667`

We multiply row 3 by 1/5:

 `T5 = IT5[3,3] = 1/5M5 = T5 %*% M4`

 `> M5 a b c d e fa 1 0 1 0.0 0.3333333 0.3333333b 0 3 1 0.0 -0.3333333 0.6666667c 0 0 1 0.6 -0.3333333 -0.5333333`

And then subtract row 3 from row 1 and row 2:

 `T6 = IT6[1,3] = -1T6[2,3] = -1M6 = T6 %*% M5`

 `> M6 a b c d e fa 1 0 0 -0.6 0.6666667 0.8666667b 0 3 0 -0.6 0.0000000 1.2000000c 0 0 1 0.6 -0.3333333 -0.5333333`

Finally, we multiply row 2 by 1/3:

 `T7 = IT7[2,2] = 1/3M7 = T7 %*% M6`

 `> M7 a b c d e fa 1 0 0 -0.6 0.6666667 0.8666667b 0 1 0 -0.2 0.0000000 0.4000000c 0 0 1 0.6 -0.3333333 -0.5333333`

The identity matrix has been transformed into M-1, as you can check:

 `result = M7[,4:6] %*% mround(result)`

 `> round(result) a b ca 1 0 0b 0 1 0c 0 0 1`

We could also do the transformations as one long chain:

 `T7 %*% T6 %*% T5 %*% T4 %*% T3 %*% T2 %*% T1 %*% M`

 `> T7 %*% T6 %*% T5 %*% T4 %*% T3 %*% T2 %*% T1 %*% M a b c d e fa 1 4.440892e-16 0 -0.6 0.6666667 0.8666667b 0 1.000000e+00 0 -0.2 0.0000000 0.4000000c 0 -4.440892e-16 1 0.6 -0.3333333 -0.5333333`

or multiply the transformation matrices together and keep them:

 `Tx = T7 %*% T6 %*% T5 %*% T4 %*% T3 %*% T2 %*% T1`

 `> round(Tx,2) [,1] [,2] [,3]a -0.6 0.67 0.87b -0.2 0.00 0.40c 0.6 -0.33 -0.53> round(Tx %*% M, 2) a b c d e fa 1 0 0 -0.6 0.67 0.87b 0 1 0 -0.2 0.00 0.40c 0 0 1 0.6 -0.33 -0.53`