Stuck at matrix inverse algorithm

Question

I want to get inverse of the matrix in pure Python and im stuck at inverse for inverse I'm using linear row reduction method. I want to know how I can perform operations on my matrix so that it becomes an identity matrix and the same operations are performed on an identity matrix so that it becomes inverse. Here is the code, note that other functions are for support

matrix1 = array([
    [2,3,4],
    [8,9,4],
    [3,4,8]
    ])


def inverse():
    # 0,0  1,1   2,2
    if len(matrix1) != len(matrix1[0]):
        return ZeroDivisionError
    #list1 = []
    flag = 0
    for i in range(len(matrix1)):
        if matrix1[i][i]  != 1:
            matrix1[i][i]  //= matrix1[i][i]
            #list1.append((i,i)) 

        # 1st col to make zero values :  1,0 , 2,0 
        # 2nd col to make zero values :  0,1 , 2,1
        # 3rd col to make zero values :  0,2 , 1,2

        for j in range(len(matrix1)):

            if i != j:

                if matrix1[i][j] != 0:
                    rowsub(matrix1[i],matrix[i][j])  

                #if flag == len(matrix1):
                #break

    return matrix1
def give_identity_mat(imatrix):
    identitymat = [[0 for i in range(len(imatrix))] for i in range(len(imatrix))]

    for i in range( 0, len(imatrix)):
        identitymat[i][i] = 1

def rowadd(row1,row2,const=1):
    #  row1 + const(row2)
    # add
    for i in range(len(matrix1[row1-1])):
        matrix1[row1-1][i] += const*(matrix1[row2-1][i])

def rowsub(row1,row2,const=1):
    for i in range(len(matrix1[row1-1])):
        matrix1[row1-1][i] -= const*(matrix1[row2-1][i])

def issquare(M):
    if len(M[0]) == len(M):
        return True
    else:
        return False

def multiply(row,num):
    for i in range(len(matrix1[row-1])):
        matrix1[row-1][i] *= num

def divide(row,num):
    for i in range(len(matrix1[row-1])):
        matrix1[row-1][i] = matrix1[row-1][i] // num

JohanC · Accepted Answer

Here is your code, with errors eliminated. It's hard to describe all the problems with the code. You seem to try to learn Python by implementing an involved algorithm. Probably it is best to start with simpler exercises.

Some advice:

Try not to directly change global variables (such as matrix1) inside other functions. Try to always give them as parameters.
In plain Python matrices are just represented by list of lists.
In Python3, "//" is used for integer division. For integer division 5 // 2 equals 2, being just the integer part of the division. For elements of a matrix you need float division, which uses one "/". 5 / 2 equals 2.5, the float version of the result.
Indices start at 0, not at 1. It is very confusing when parameters are numbered starting from 1. Best just stick to starting with 0.
I didn't implement yet the step where you have to search for a non-zero element in the column.
For the algorithm you need two matrices: one given input matrix and an output matrix that is built up step by step.

def test_calculating_the_inverse_matrix():
    matrix1 = [
        [2, 3, 4],
        [8, 9, 4],
        [3, 4, 8]
    ]
    inv = inverse(matrix1)
    print("The inverse matrix is:")
    for row in inv:
        print ("  ", row)

def inverse(matrix1):
    # calculate the inverse matrix of matrix1
    # algorithm:
    # - initialize invMatrix with a identity matrix of the same size as matrix1
    # - do gaussian elimination on both matrices simultaneously, converting matrix1 to an identity matrix while
    #     buidling up the inverse matrix

    if not issquare(matrix1):
        return ZeroDivisionError
    numRows = len(matrix1)
    invMatrix = give_identity_mat(numRows)

    for i in range(numRows):
        if matrix1[i][i] == 0:
            print("TODO: when a zero on the diagonal is encountered, search for a non-zero below in the same column")
            print("TODO: if a non-zero is found, exchange both rows")
            print("TODO: if no non-zero is found, the matrix has no inverse")
            return ZeroDivisionError

        if matrix1[i][i] != 1:   # we have to divide the complete row by matrix1[i][i], for both matrices
            factor = matrix1[i][i]
            divide(matrix1, i, factor)
            divide(invMatrix, i, factor)
        print(f"after dividing row {i}:")
        print("  ", matrix1)
        print("  ", invMatrix)

        for j in range(numRows):
            if i != j:
                if matrix1[j][i] != 0:
                    factor = matrix1[j][i]
                    rowsub(matrix1, i, j, factor)
                    rowsub(invMatrix, i, j, factor)
                    print(f"after subtracting row {i} from row {j}:")
                    print("  ", matrix1)
                    print("  ", invMatrix)
    return invMatrix

def give_identity_mat(num_rows):
    identitymat = [[0 for i in range(num_rows)] for i in range(num_rows)]
    for i in range( 0, num_rows):
        identitymat[i][i] = 1
    return identitymat

def rowsub(matrix, row1, row2, factor):
    # subtract row1 from row2
    for i in range(len(matrix[row2])):
        matrix[row2][i] -= factor * matrix[row1][i]

def issquare(M):
    return len(M[0]) == len(M)

def multiply(matrix, row, factor):
    for i in range(len(matrix[row])):
        matrix[row][i] *= factor

def divide(matrix, row, factor):
    for i in range(len(matrix[row])):
        matrix[row][i] /= factor


test_calculating_the_inverse_matrix()

Output:

The inverse matrix is:
   [-2.3333333333333335, 0.33333333333333337, 1.0]
   [2.1666666666666665, -0.16666666666666666, -1.0]
   [-0.20833333333333334, -0.041666666666666664, 0.25]

Note that due to inevitable rounding errors, the product of the original matrix with its inverse is almost never the exact identity matrix. It will be an approximation.

Stuck at matrix inverse algorithm

Answers (1)

Related Questions