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 c a 2 1 4 b 2 -3 1 c 1 3 2

The identity matrix is I:

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

> I a b c a 1 0 0 b 0 1 0 c 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 = I temp = Tx[1,] Tx[1,] = Tx[3,] Tx[3,] = temp

> Tx a b c a 0 0 1 b 0 1 0 c 1 0 0 > Tx %*% m a b c a 1 3 2 b 2 -3 1 c 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 f a 2 1 4 1 0 0 b 2 -3 1 0 1 0 c 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 = I T2[2,2] = -1/3 M2 = T2 %*% M1

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

We add -1 x row 2 to row 1:

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

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

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

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

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

We multiply row 3 by 1/5:

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

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

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

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

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

Finally, we multiply row 2 by 1/3:

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

> M7 a b c d e f a 1 0 0 -0.6 0.6666667 0.8666667 b 0 1 0 -0.2 0.0000000 0.4000000 c 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] %*% m round(result)

> round(result) a b c a 1 0 0 b 0 1 0 c 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 f a 1 4.440892e-16 0 -0.6 0.6666667 0.8666667 b 0 1.000000e+00 0 -0.2 0.0000000 0.4000000 c 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.87 b -0.2 0.00 0.40 c 0.6 -0.33 -0.53 > round(Tx %*% M, 2) a b c d e f a 1 0 0 -0.6 0.67 0.87 b 0 1 0 -0.2 0.00 0.40 c 0 0 1 0.6 -0.33 -0.53