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}
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,
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]:
and we suppose a $n \times n$ Identity matrix
In [48]:
I_n = eye(n)
Out[48]:
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]:
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]:
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]:
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]:
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]:
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]:
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]:
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]:
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]:
We concatenate the matrices
In [60]:
M = hcat(U,B)
Out[60]:
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]:
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]:
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]:
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]:
Finally we normalize the first diagonal entry
In [65]:
M[1,:] = M[1,:]/M[1,1]
M
Out[65]:
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]:
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)))
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))
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))