Here is my solution for Lab 7:
from numpy import *
def print_matrix(M):
print("The matrix is currently:")
print(array(M))
print("="*80)
def get_lead_ind(arr):
for i in range(len(arr)):
if arr[i] != 0:
return i
return len(arr)
def get_row_to_swap(M, start_i):
best_lead_ind = len(M[0])
best_i = len(M[0])
for i in range(start_i, len(M)):
lead_ind = get_lead_ind(M[i])
if lead_ind < best_lead_ind:
best_i = i
best_lead_ind = lead_ind
return best_i, best_lead_ind
def subtract_rows_coefs(r1, c1, r2, c2):
r = [0]*len(r1)
for i in range(len(r1)):
r[i] = c1*r1[i]-c2*r2[i]
return r
def subtract_rows(M, row_to_sub, best_lead_ind):
for row1 in range(row_to_sub+1, len(M)):
if M[row1][best_lead_ind] != 0:
M[row1] = subtract_rows_coefs(M[row1], 1,
M[row_to_sub], M[row1][best_lead_ind]/M[row_to_sub][best_lead_ind] )
def div_lead_coef(r):
i = get_lead_ind(r)
if i < len(r):
ri = r[i]
for j in range(i, len(r)):
r[j] /= ri
def forward_step(M):
for row in range(len(M)):
print("Now looking at row %d" % (row))
best_i, best_lead_ind = get_row_to_swap(M, row)
if best_i == len(M[0]):
break
print("Swapping rows %d and %d so that entry %d in the current row is non-zero" % (row, best_i, best_lead_ind))
M[row], M[best_i] = M[best_i], M[row]
print_matrix(M)
print("Adding row %d to rows below it to eliminate coefficients in column %d" % (row, best_lead_ind))
subtract_rows(M, row, best_lead_ind)
print_matrix(M)
def subtract_rows_back(M, row_to_sub, best_lead_ind):
for row1 in range(row_to_sub):
if M[row1][best_lead_ind] != 0:
M[row1] = subtract_rows_coefs(M[row1], 1,
M[row_to_sub], M[row1][best_lead_ind]/M[row_to_sub][best_lead_ind] )
def backward_step(M):
for row in range(len(M)-1, -1, -1):
lead_ind = get_lead_ind(M[row])
print("Adding row %d to rows above it to eliminate coefficients in column %d" % (row, lead_ind))
subtract_rows_back(M, row, lead_ind)
print_matrix(M)
##set b to the basis vectors
def solve(M, b):
aug = []
for i in range(len(M)):
aug.append(M[i]+[b[i]])
print_matrix(aug)
print("Now performing the forward step")
forward_step(aug)
print_matrix(aug)
print("Now performing the backward step")
backward_step(aug)
print("Now dividing each row by the leading coefficient")
for i in range(len(aug)):
div_lead_coef(aug[i])
print_matrix(aug)
x = []
for i in range(len(aug)):
x.append(aug[i][-1])
return x
Let's try it with a regular matrix:
M = array([[1,-2,3],
[3,10,1],
[1,5,3]])
x = array([75,10,-11])
b = dot(M,x)
print("x =", solve(M.tolist(), b.tolist()))
Everything seems to work fine!
Now, let's try this with a different matrix:
n = 3
M = zeros((n, n)) #[[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
for i in range(n):
for j in range(n):
M[i][j] = 1/(i+j+1)
M is a special kind of matrix -- called a Hilbert Matrix. Let's try to solve a system of equations where the coefficients form a Hilbert Matrix
x = array([75,10,-11])
b = dot(M,x)
print("x =", solve(M.tolist(), b.tolist()))
Oops! What happenned? When we were "Adding row 1 to rows below it to eliminate coefficients in column 1", there was a glitch: because floats are stored imprecisely, we didn't eliminate the coefficient in the last row.
That we can fix:
def subtract_rows(M, row_to_sub, best_lead_ind):
for row1 in range(row_to_sub+1, len(M)):
if M[row1][best_lead_ind] != 0:
M[row1] = subtract_rows_coefs(M[row1], 1,
M[row_to_sub], M[row1][best_lead_ind]/M[row_to_sub][best_lead_ind] )
M[row1][best_lead_ind] = 0 #It's already supposed to be 0, but we're just making extra-sure.
#This makes it so that imprecision that makes M[row1][best_lead_ind]
#not quite 0 doesn't affect further calculations
Let's try again:
x = array([75,10,-11])
b = dot(M,x)
print("x =", solve(M.tolist(), b.tolist()))
Better, but there is still a problem! Imprecision that was introduced because floats cannot be stored precisely at some stage of the computation made it so that the final answer isn't quite precise.
Why did it happen particularly with the Hilbert Matrix? Because the computation of the inverse of the Hilbert Matrix is ill-conditioned -- small changes in the input make for large changes in the output.
Let's take a step back, and thing more about solving $Mx = b$. Let's solve three equations:
$M x_1 = \begin{pmatrix}1\\0\\0\end{pmatrix}$, $M x_2 = \begin{pmatrix}0\\1\\0\end{pmatrix}$, $M x_3 = \begin{pmatrix}0\\0\\1\end{pmatrix}$
You can see (after a bit of thought) that $M (x1 | x2 | x3) = ( x1 | x2 | x3) M = \begin{pmatrix}1 & 0 & 0\\0 & 1& 0\\0 & 0 & 1\end{pmatrix}$
Here, $(x1 | x2 | x3)$ is a matrix whose columns are $x_1$, $x_2$, and $x_3$. That means we can write
$$M^{-1} = (x1 | x2 | x3)$$And further, if $Mx = b$, then $x = M^{-1}b$.
So one way to solve the equation $Mx = b$ is to first solve for $x_1, x_2, x_3$, and then use $M^{-1}$. What we do with Gaussian elimination is actually similar in spirite.
Let's see what $x_1, x_2, x_3$ are.
print("x =", solve(M.tolist(), [1, 0, 0]))
print("x =", solve(M.tolist(), [0, 1, 0]))
print("x =", solve(M.tolist(), [0, 0, 1]))
As you can see, the numbers are quite lairge. Small changes in $M$ will produce large changes in $M^{-1}$ because of that -- that is a property of Hilbert Matrices. A rough analogy is thinking about computing $y^{-1} = 1/y$ for $y$ that is close to 0. Small changes in $y$ will produces large changes in $1/y$.
Some matrices are ill-conditioned -- that means small changes in them will produce large changes in their inverse. This makes it hard (and sometimes very hard) to apply standard algorithms such as Gaussian Elimination to them, since small errors that are due to the imprecision in floats in the intermediate steps will tend to accumulate.