Worksheet 5

In this worksheet we will compute the inverse of a non-singular matrix using the Gauss-Jordan Method.

Computing the inverse of a matrix is extremely usefull when you need to compute the solution of the linear system \begin{equation} Ax_j = b_j, \qquad \mbox{for } j = 1, ..., R. \end{equation} I.e. you need to solve the same system for many different right-hand sides.

If you know the inverse of $A$, denoted $A^{-1}$, then you can easily solve the system for any right-hand side by applying the inverse, $A^{-1}$, which results in \begin{equation} x_j = A^{-1}b_j. \end{equation}

Gauss-Jordan method

You know that a $n\times n$ matrix $A$ is invertible if and only if, it is row equivalent to an identity matrix, which in in this case is denoted $I_n$.

You know that if we apply the Gauss-Jordan method to the augmented matrix $[A, I_n]$, the resulting matrix will have the form $[I_n, A^{-1}]$. Where $A^{-1}$ is the inverse of $A$.

In this worksheet we will show you, step by step, how to apply the Gauss-Jordan method to the augmented matrix $[A, I_n]$, where $n = 4$, and $A$ is a $n \times n$ randomly generated matrix.

For the sake of presentation we will split the computation in two stages,

  • a downward pass that corresponds to the Gaussian Eliminiation,
  • and an upwards pass, to obtain the reduced row echelon form. In this stage, the elements above the diagonal will be eliminated and the diagonal elements will be normalized to one.

In [47]:
n = 4
A = rand(n,n) + n*eye(n) # we add a identity to be sure that the matrix is well conditioned


Out[47]:
4x4 Array{Float64,2}:
 4.1951     0.908264  0.855442   0.618706 
 0.258505   4.46066   0.0740031  0.318408 
 0.0332677  0.510741  4.26229    0.0569284
 0.103489   0.411466  0.808493   4.73789  

and we suppose a $n \times n$ Identity matrix


In [48]:
I_n = eye(n)


Out[48]:
4x4 Array{Float64,2}:
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0

In order to perform the Gaussian elimination we build the augmented matrix


In [49]:
M = hcat(A,I_n) # or in equivalent fashion M = [A b]


Out[49]:
4x8 Array{Float64,2}:
 4.1951     0.908264  0.855442   0.618706   1.0  0.0  0.0  0.0
 0.258505   4.46066   0.0740031  0.318408   0.0  1.0  0.0  0.0
 0.0332677  0.510741  4.26229    0.0569284  0.0  0.0  1.0  0.0
 0.103489   0.411466  0.808493   4.73789    0.0  0.0  0.0  1.0

We start the first step of the Gaussian elimination. We look for the pivot, i.e. the first non-zero in the first column


In [50]:
i = 1; #first column
find(M[i:end,i])


Out[50]:
4-element Array{Int64,1}:
 1
 2
 3
 4

The function find provides the indices of all the non zeros elements in the input vector. In this case all of them are non zero, so we don't need to perform a row swap. Then we can start introducing zeros in the first column, below the diagonal. In order to do this we perform a row operation, given by


In [51]:
M[2,:] = M[2,:] - M[2,1]/M[1,1]*M[1,:]
M


Out[51]:
4x8 Array{Float64,2}:
 4.1951     0.908264  0.855442   0.618706    1.0        0.0  0.0  0.0
 0.0        4.4047    0.0212902  0.280283   -0.0616207  1.0  0.0  0.0
 0.0332677  0.510741  4.26229    0.0569284   0.0        0.0  1.0  0.0
 0.103489   0.411466  0.808493   4.73789     0.0        0.0  0.0  1.0

As you can see we introduced a zero in (2,1). We can repeat this operation with a for for the rest of the column


In [52]:
for j = 1+1:n
    M[j,:] = M[j,:] - M[j,1]/M[1,1]*M[1,:]
end
M


Out[52]:
4x8 Array{Float64,2}:
 4.1951  0.908264  0.855442   0.618706   1.0         0.0  0.0  0.0
 0.0     4.4047    0.0212902  0.280283  -0.0616207   1.0  0.0  0.0
 0.0     0.503539  4.2555     0.052022  -0.00793013  0.0  1.0  0.0
 0.0     0.38906   0.78739    4.72263   -0.0246689   0.0  0.0  1.0

You can see that we introduced zeros in the whole first column under the diagonal.
Now, we can start with the second row.
We check for the first non-zero in the second column, under the diagonal.


In [53]:
indnz = find(M[2:n,2])[1] # you pick the first index


Out[53]:
1

indz is the first non-zero index in the 2nd column starting from the 2nd row (i.e. is a local index) we need to transorm it back to a global index


In [54]:
indj = 1 + indnz


Out[54]:
2

If indj is different from i, this means that $M_{i,i} = 0$ and we need to change the pivot using a row swap


In [55]:
if indj != 2
    buffer = M[2,:]
    M[2,:] = M[indj,:]
    M[indj,:]  = buffer
end

Now we are sure that we have a non zero pivot, so we can continue to introduce zeros


In [56]:
for j = 3:n
    M[j,:] = M[j,:] - M[j,2]/M[2,2]*M[2,:]
end
M


Out[56]:
4x8 Array{Float64,2}:
 4.1951  0.908264     0.855442   …   1.0           0.0        0.0  0.0
 0.0     4.4047       0.0212902     -0.0616207     1.0        0.0  0.0
 0.0     0.0          4.25307       -0.000885742  -0.114319   1.0  0.0
 0.0     5.55112e-17  0.78551       -0.0192261    -0.0883284  0.0  1.0

Now we can go to third column


In [57]:
indnz = find(M[3:n,3])[1] # you pick the first index
indj = 2 + indnz
if indj != 3
    buffer = M[3,:]
    M[3,:] = M[indj,:]
    M[indj,:] = buffer
end
for j = 4:n
    M[j,:] = M[j,:] - M[j,3]/M[3,3]*M[3,:]
end
M


Out[57]:
4x8 Array{Float64,2}:
 4.1951  0.908264     0.855442   0.618706   …   0.0         0.0       0.0
 0.0     4.4047       0.0212902  0.280283       1.0         0.0       0.0
 0.0     0.0          4.25307    0.0199805     -0.114319    1.0       0.0
 0.0     5.55112e-17  0.0        4.69418       -0.0672147  -0.184693  1.0

Now we can extract $U$ and $L^{-1}b = b1$ from the augmented matrix


In [58]:
U = M[1:n,1:n]
B = M[1:n,n+1:end]


Out[58]:
4x4 Array{Float64,2}:
  1.0           0.0         0.0       0.0
 -0.0616207     1.0         0.0       0.0
 -0.000885742  -0.114319    1.0       0.0
 -0.0190625    -0.0672147  -0.184693  1.0

Reduction to Reduced Row Echelon Form

Now we will reduce the augmented matrix to a reduced row echelon form, by eliminating the non-zero coefficients above the diagonal from the bottom to the top.

We concatenate the matrices


In [60]:
M = hcat(U,B)


Out[60]:
4x8 Array{Float64,2}:
 4.1951  0.908264     0.855442   0.618706   …   0.0         0.0       0.0
 0.0     4.4047       0.0212902  0.280283       1.0         0.0       0.0
 0.0     0.0          4.25307    0.0199805     -0.114319    1.0       0.0
 0.0     5.55112e-17  0.0        4.69418       -0.0672147  -0.184693  1.0

We already know that the matrix is invertible (explain why).

We normalize the last diagonal element by


In [61]:
M[n,:] = M[n,:]/M[n,n]


Out[61]:
1x8 Array{Float64,2}:
 0.0  1.18255e-17  0.0  1.0  -0.00406088  -0.0143187  -0.039345  0.21303

and we perform a for loop to eliminate the entried above the diagonal


In [62]:
for i=n-1:-1:1
    M[i,:] = M[i,:] - M[i,n]*M[n,:]
end
M


Out[62]:
4x8 Array{Float64,2}:
 4.1951   0.908264     0.855442   0.0  …   0.024343   -0.131803  
 0.0      4.4047       0.0212902  0.0      0.0110277  -0.0597085 
 0.0     -2.3628e-19   4.25307    0.0      1.00079    -0.00425643
 0.0      1.18255e-17  0.0        1.0     -0.039345    0.21303   

and for the second to last column


In [63]:
M[n-1,:] = M[n-1,:]/M[n-1,n-1]
for i=(n-1)-1:-1:1
    M[i,:] = M[i,:] - M[i,n-1]*M[n-1,:]
end
M


Out[63]:
4x8 Array{Float64,2}:
 4.1951   0.908264     0.0  0.0   1.00267      …  -0.17695     -0.130947  
 0.0      4.4047       0.0  0.0  -0.0604785        0.00601795  -0.0596872 
 0.0     -5.55551e-20  1.0  0.0  -0.000189182      0.235309    -0.00100079
 0.0      1.18255e-17  0.0  1.0  -0.00406088      -0.039345     0.21303   

and the third to last column


In [64]:
M[n-2,:] = M[n-2,:]/M[n-2,n-2]
for i=(n-2)-1:-1:1
    M[i,:] = M[i,:] - M[i,n-2]*M[n-2,:]
end
M


Out[64]:
4x8 Array{Float64,2}:
 4.1951   0.0          0.0  0.0   1.01515      …  -0.178191    -0.118639  
 0.0      1.0          0.0  0.0  -0.0137305        0.00136626  -0.0135508 
 0.0     -5.55551e-20  1.0  0.0  -0.000189182      0.235309    -0.00100079
 0.0      1.18255e-17  0.0  1.0  -0.00406088      -0.039345     0.21303   

Finally we normalize the first diagonal entry


In [65]:
M[1,:] = M[1,:]/M[1,1]
M


Out[65]:
4x8 Array{Float64,2}:
 1.0   0.0          0.0  0.0   0.241984     …  -0.0424761   -0.0282804 
 0.0   1.0          0.0  0.0  -0.0137305        0.00136626  -0.0135508 
 0.0  -5.55551e-20  1.0  0.0  -0.000189182      0.235309    -0.00100079
 0.0   1.18255e-17  0.0  1.0  -0.00406088      -0.039345     0.21303   

We can then split the augmented matrix, in the quasi identity matrix and the inverse of A,


In [66]:
I = M[1:n,1:n]
Ainv = M[1:n,n+1:end]


Out[66]:
4x4 Array{Float64,2}:
  0.241984     -0.0417997  -0.0424761   -0.0282804 
 -0.0137305     0.228071    0.00136626  -0.0135508 
 -0.000189182  -0.0268118   0.235309    -0.00100079
 -0.00406088   -0.0143187  -0.039345     0.21303   

We can check how far from the indentity I is by computing


In [67]:
println("Distance in L2 norm to the identity is ", norm(I-eye(n)))


Distance in L2 norm to the identity is 1

We can test the Ainv is the inverse of A by testing


In [68]:
@assert norm(A*Ainv- eye(n)) < n*eps(1.0)

The macro assert will test if your condition is true, if it is true it will continue otherwise it will raise an error.

We can use Ainv to solve linear systems fast (after Ainv has been computed)


In [71]:
b = randn(n,1)
x = Ainv*b;
println("Residual of the solution computed with the inverse is ", norm(A*x - b))


Residual of the solution computed with the inverse is 2.5438405243138006e-16

In particular you can solve multiple righ-hand sides simultaneously. In such cases you can stack the right-hand sides in a matrix of the form \begin{equation} B = [b_1| b_2| ...| b_R], \end{equation} and you can stack the solutions \begin{equation} X = [x_1| x_2| ...| x_R]. \end{equation} With the notation introduced above you can write the solution of the muliple linear systems as \begin{equation} AX=B. \end{equation}

If you have the inverse you can easily find the solution, by applying $A^{-1}$ to $B$.


In [73]:
B = randn(n,5)
X = Ainv*B;
println("Residual of the solution computed with the inverse is ", norm(A*X - B))


Residual of the solution computed with the inverse is 5.731849917315797e-16