psie
- 315
- 40
- TL;DR Summary
- I've been working on an elementary linear algebra exercise, namely that every matrix can be transformed into an upper triangular matrix by using elementary row operations of type 1 and 3 only, that is, by swapping rows and adding a multiple of a row to another. Now, I've tried to implement this in python, and after some struggling, I think I've succeeded, though I don't really know how to be sure 100%.
First, I'm a beginner at this, so sorry if my code looks stupid. I'd be grateful for any feedback concerning this. I wonder if anything's redundant or can be improved?
The idea of the "algorithm" is to start with the first column of the matrix, in particular the first entry below the diagonal entry. Then check the entries in this column to see if there's a nonzero entry. If there is, we swap it with the row containing the diagonal entry (by using type 1 operation) and then use type 3 operation to subtract this nonzero entry (multiplied by the appropriate constant) from all the other nonzero entries in this column. Then we proceed to the next column and repeat this process.
I've checked by hand that the matrices p and q below give the correct outputs (also trying other, bigger matrices).
The idea of the "algorithm" is to start with the first column of the matrix, in particular the first entry below the diagonal entry. Then check the entries in this column to see if there's a nonzero entry. If there is, we swap it with the row containing the diagonal entry (by using type 1 operation) and then use type 3 operation to subtract this nonzero entry (multiplied by the appropriate constant) from all the other nonzero entries in this column. Then we proceed to the next column and repeat this process.
I've checked by hand that the matrices p and q below give the correct outputs (also trying other, bigger matrices).
Python:
import numpy as np
p=np.array([[1,0,1],
[1,1,1],
[0,0,1]])
q=np.array([[1,2],
[2,3],
[2,2]])
def f(k):
if not isinstance(k,np.ndarray):
return print("input is not a numpy array")
elif len(k.shape)!=2:
return print("numpy array is not correct shape")
else:
for j in range(0,k.shape[1]):
#if matrix is square, we don't need j=k.shape[1]-1
if k.shape[0]==k.shape[1] and j==k.shape[1]-1:
break
else:
for i in range(j+1,k.shape[0]):
if k[i,j]!=0:
#swap rows
k[[j,i]]=k[[i,j]]
#this for-loop subtracts the nonzero entry from the rest in a given column
for s in range(i,k.shape[0]):
firstentry=k[s,j]
c=-firstentry/k[j,j]
current=k[[s]]
k[[s]]=current+k[[j]]*c
#maybe this is redudant, but if k[i,j]=0 we should continue unless i=k.shape[0]-1
elif k[i,j]==0 and i==k.shape[0]-1:
break
else:
continue
return k
print(f(p),'\n')
print(f(q))