Simple code example of transforming matrix into upper triangular form

[psie]

TL;DR

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).

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))

 

[pbuk]

First, I'm a beginner at this, so sorry if my code looks stupid. I'd be grateful for any feedback concerning this.

Use spaces around binary operators and after ,s to help readability (see the PEP8 coding standard and use an IDE tool that enforces it).

[code lang=Python title=Bad]
for s in range(i,k.shape[0]):
firstentry=k[s,j]
c=-firstentry/k[j,j]
current=k[]
k[]=current+k[[j]]*c
[/code]

[code lang=Python title=Better]
for s in range(i, k.shape[0]):
firstentry = k[s, j]
c = -firstentry / k[j, j]
current = k[]
k[] = current + k[[j]] * c
[/code]

• Use descriptive names for functions.
• Use descriptive names for variables: it is conventional to use M for a matrix, and it is much easier to avoid mistake if you create a variable e.g. nRows than to refer to M.shape[0].
• Do not use redundant else conditions
• Use spaces around binary operators and after ,s to help readability (see the PEP8 coding standard and use an IDE tool that enforces it).

[code lang=Python title=Bad]
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]):
[/code]

[code lang=Python title=Better]
def diagonalize(M):

if not isinstance(M, np.ndarray):
return print("input must be a numpy array")

if len(M.shape) != 2:
return print("input must have 2 dimensions")

[nRows, nCols] = M.shape

for col in range(0, nCols):
[/code]

 

[pbuk]

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.

As far as the algorithm itself goes, working this up from first principles is a worthwhile exercice but only in the context of a structured study of numerical methods. Are you following such a course of study?

If not I suggest these lecture notes: https://people.sc.fsu.edu/~jburkardt/classes/nla_2015/numerical_linear_algebra.pdf

 

[psie]

As far as the algorithm itself goes, working this up from first principles is a worthwhile exercice but only in the context of a structured study of numerical methods. Are you following such a course of study?

If not I suggest these lecture notes: https://people.sc.fsu.edu/~jburkardt/classes/nla_2015/numerical_linear_algebra.pdf

Thanks for the feedback and help! Also the recommendation. I am currently just reading my linear algebra book (a mathematics book) where this was just an exercise in describing the procedure of how to go about it. I immediately felt for trying to do something more with it.