Computing the Inverse of a Matrix via Row Reduction

Some Preliminary Functions

You don't need to see these, but they will generate the random matrices I'll use in my computations.

In[39]:=

Rand[x_, y_] := Random[] ; RandomMatrix[n_] := Array[Rand, {n, n}] ; RandomMatrix[n_, m_] := Array[Rand, {n, m}] ;

Some examples with 3x3 Matrices

I want to implement the procedure for calculating the inverse of an operator which we talked about in class.  We start by selecting a random 3x3 matrix which we'll call A.  The procedure I've written to generate a random matrix gives an n x n matrix with entries selected randomly between 0 and 1.  

In[40]:=

A = RandomMatrix[3]

Out[40]=

( {{0.170254, 0.646524, 0.0262419}, {0.165057, 0.902071, 0.976832}, {0.0732227, 0.00987297, 0.250392}} )

Now we augment the matrix A with the identity matrix.  It takes some work to get Mathematica to do this, but you shouldn't worry about that.  

In[41]:=

AugmentedA = Transpose[Join[Transpose[A], IdentityMatrix[3]]]

Out[41]=

( {{0.170254, 0.646524, 0.0262419, 1, 0, 0}, {0.165057, 0.902071, 0.976832, 0, 1, 0}, {0.0732227, 0.00987297, 0.250392, 0, 0, 1}} )

Finally,we perform elementary row operations until our matrix is in reduced row echelon form. RowReduce[ ] is the function Mathematica uses to put a matrix in reduced row echelon form.

In[42]:=

AugmentedAInv = RowReduce[AugmentedA]

Out[42]=

( {{1, 0., 0., 3.95685, -2.95766, 11.1238}, {0, 1, 0., 0.552596, 0.74495, -2.96412}, {0, 0, 1, -1.1789, 0.835541, 0.857667}} )

Since the left hand side is the identity matrix,the right hand side is the inverse of our matrix.

In[43]:=

AInv = Transpose[Drop[Transpose[AugmentedAInv], 3]]

Out[43]=

( {{3.95685, -2.95766, 11.1238}, {0.552596, 0.74495, -2.96412}, {-1.1789, 0.835541, 0.857667}} )

Later on when we talk about matrix multiplication we'll see that the inverse we've just calculated is the multiplicative inverse of A.  That is to say, when I multiply A and AInv I should get the identity matrix.  Let's check that here.

In[44]:=

A . AInv

Out[44]=

( {{1., -5.89806*10^-17, -1.45717*10^-16}, {-2.22045*10^-16, 1., 4.44089*10^-16}, {2.22045*10^-16, -1.66533*10^-16, 1.}} )

In[45]:=

AInv . A

Out[45]=

( {{1., -3.88578*10^-16, -2.22045*10^-15}, {-5.55112*10^-17, 1., 5.55112*10^-16}, {1.38778*10^-17, -5.20417*10^-18, 1.}} )

The numbers off the diagonal which look nonzero are just roundoff error.  You can see that they are essentially zero.
We can do this again with another randomly selected 3x3 matrix.  (We'll just run the commands over)

In[46]:=

B = RandomMatrix[3]

Out[46]=

( {{0.693553, 0.69199, 0.430368}, {0.432805, 0.5849, 0.398616}, {0.741279, 0.439683, 0.0465146}} )

Now we augment.

In[47]:=

AugmentedB = Transpose[Join[Transpose[B], IdentityMatrix[3]]]

Out[47]=

( {{0.693553, 0.69199, 0.430368, 1, 0, 0}, {0.432805, 0.5849, 0.398616, 0, 1, 0}, {0.741279, 0.439683, 0.0465146, 0, 0, 1}} )

Finally,we perform elementary row operations until our matrix is in reduced row echelon form.

In[48]:=

AugmentedBInv = RowReduce[AugmentedB]

Out[48]=

( {{1, 0., 0., 8.79071, -9.32385, -1.43186}, {0, 1, 0., -16.3487, 17.026, 5.35522}, {0, 0, 1, 14.4441, -12.3505, -6.30318}} )

Since the left hand side is the identity matrix,the right hand side is the inverse of our matrix.

In[49]:=

BInv = Transpose[Drop[Transpose[AugmentedBInv], 3]]

Out[49]=

( {{8.79071, -9.32385, -1.43186}, {-16.3487, 17.026, 5.35522}, {14.4441, -12.3505, -6.30318}} )

Again, let's check that this is actually the inverse.

In[50]:=

B . BInv

Out[50]=

( {{1., -1.77636*10^-15, -8.88178*10^-16}, {-1.77636*10^-15, 1., -4.44089*10^-16}, {5.55112*10^-16, -1.66533*10^-15, 1.}} )

In[51]:=

BInv . B

Out[51]=

( {{1., -8.88178*10^-16, -4.44089*10^-16}, {4.44089*10^-16, 1., -9.71445*10^-16}, {-8.88178*10^-16, 0., 1.}} )

Aside from these roundoff errors, we get the identity matrix!

A bigger example - an 8x8 matrix

Let's do a larger computation just for fun.  Notice that without a computer, I wouldn't come close enough to an 8 x 8 matrix to touch it with a 10 foot pole, let alone row reduce it.

In[52]:=

F = RandomMatrix[8]

Out[52]=

Our procedure says we should augment F with the identity matrix.

In[53]:=

AugmentedF = Transpose[Join[Transpose[F], IdentityMatrix[8]]]

Out[53]=

Finally,we perform elementary row operations until our matrix is in reduced row echelon form.

In[54]:=

AugmentedFInv = RowReduce[AugmentedF]

Out[54]=

Since the left hand side is the identity matrix,the right hand side is the inverse of our matrix.

In[55]:=

FInv = Transpose[Drop[Transpose[AugmentedFInv], 8]]

Out[55]=

Let's check that this is actually the inverse

In[56]:=

F . FInv

Out[56]=

In[57]:=

FInv . F

Out[57]=

We did it!