Class, standalone functions and miscellaneous methods (1/3)

The Mat(rix) Class

First off, we need a class 'Mat' (as in 'matrix'). Matrix like data will be passed to the Mat class in the form of a list of lists, with each sub-list acting as a row. We will gradually build up a suit of methods to do all the needed operations.

class Mat:
    def __init__(self, data):
        self.data = data

# define a 2x3 matrix
Mat([[1, 2, 3], [4, 5, 6]])

Standalone functions

Printing matrices

Here is a function to print a matrix, one row beneath the other. This is helpful for debugging. We also include an option to round the values to a specified precision, which is useful to ignore very tiny rounding errors on what would otherwise be integer numbers:

Code implementation

def print_mat(A, round_dp=99):
    for row in A.data:
        rounded = [round(j,round_dp) for j in row]
        print(rounded)
    print()

Generate matrices

We make a method to generate a matrix of a given size, populated with a particular set of values (default all zeros) and of a particular type (full by default, or optionally diagonal, upper, lower):

First, we ensure that we won't run out of values when populating the matrix. If a single value is passed, then we repeat that value. If more than one value is passed, but is less than the largest row/column dimension of the matrix, then we pad the value list with zeros. This scheme will allow us to make some quite useful types of matrices such as tri-diagonal etc.

def gen_mat(size, values=[0], family='full'):
    if type(size) is int:
        size = [size, size]
    if len(values) == 1:
        values = [values[0] for val in range(max(size))]
    elif len(values) < max(size):
        values += [0 for val in range(max(size)-len(values))]

Now that we have our list of values, we loop over each row/column of the matrix and decide whether to put a value from the list or a zero (in the case of diagonal/upper/lower matrices, we know some parts of the matrix will be zero).

If the matrix is to have a single diagonal, we place a value from the list once per row, in the jth column.

An interesting case is when we provide more than one value with family='full'. In this case we check which diagonal we're on (i.e. main, 1st off-diagonal, 2nd off-diagonal etc.) and pull the corresponding item (0th, 1st, 2nd) from the value list. If the original value list contained more than one value it will have been padded with zeros, and the result will be a matrix with some recurring number along each diagonal (see demos below). If only a single value was supplied then this value will have been repeated and the result will be a uniformly populated matrix

generated_mat = []
for i in range(size[0]):
    row = []
    for j in range(size[1]):
        if (family == 'diag' and j!=i) or (family == 'upper' and j<=i) or (family == 'lower' and j>=i):
            row.append(0)
        elif family == 'diag':
            row.append(values[j])
        elif j>=i:
            row.append(values[j-i])
        elif j<i:
            row.append(values[i-j])
    generated_mat.append(row)
return Mat(generated_mat)

A useful matrix to generate immediately is the identity matrix. We use the previous function and populate the diagonal with 1's.

def eye(size):
    return gen_mat(size, values=[1], family='diag')

Demo

import linalg as la

la.print_mat(la.gen_mat([7,7], values=[1], family='full'))
la.print_mat(la.gen_mat([7,7], values=[1], family='upper'))
la.print_mat(la.gen_mat([7,7], values=[1,2,0,4,5], family='diag'))
la.print_mat(la.gen_mat([7,7], values=[2,-1], family='full'))

Outputs:

>>> la.print_mat(la.gen_mat([7,7], values=[1], family='full'))
[1, 1, 1, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1, 1]

>>> la.print_mat(la.gen_mat([7,7], values=[1], family='upper'))

[0, 1, 1, 1, 1, 1, 1]
[0, 0, 1, 1, 1, 1, 1]
[0, 0, 0, 1, 1, 1, 1]
[0, 0, 0, 0, 1, 1, 1]
[0, 0, 0, 0, 0, 1, 1]
[0, 0, 0, 0, 0, 0, 1]
[0, 0, 0, 0, 0, 0, 0]

>>> la.print_mat(la.gen_mat([7,7], values=[1,2,0,4,5], family='diag'))

[1, 0, 0, 0, 0, 0, 0]
[0, 2, 0, 0, 0, 0, 0]
[0, 0, 0, 0, 0, 0, 0]
[0, 0, 0, 4, 0, 0, 0]
[0, 0, 0, 0, 5, 0, 0]
[0, 0, 0, 0, 0, 0, 0]
[0, 0, 0, 0, 0, 0, 0]

>>> la.print_mat(la.gen_mat([7,7], values=[2,-1], family='full'))

[2, -1, 0, 0, 0, 0, 0]
[-1, 2, -1, 0, 0, 0, 0]
[0, -1, 2, -1, 0, 0, 0]
[0, 0, -1, 2, -1, 0, 0]
[0, 0, 0, -1, 2, -1, 0]
[0, 0, 0, 0, -1, 2, -1]
[0, 0, 0, 0, 0, -1, 2]

The vandamonde matrix

The Vandamonde matrix stores the components of an nth degree polynomial curve. It has as many rows as points on the curve and one column per term of the polynomial - the 1st column is all ones, and is the basis for the intercept ($x^0 = 1$); The second incrementally increases and forms the basis for the slope ($x^1 = x$); The third column is for the squared terms ($x^2 = x^2$) and so on. When a vector $b$ is [projected onto](../projection_and_regression.md) the column space of this matrix, the resulting projection is the 'best fitting' polynomial curve through the points $b$:

def vandermonde(n_rows, order=1):
    A = gen_mat([n_rows, 1])
    for i in range(n_rows):
        orders = []
        for exponent in range(order+1):
            orders.append(i**exponent)
        A.data[i] = orders
    return A

Demo

We generate a vandermonde matrix with terms up to cubic and print the result:

V = la.vandermonde(5, order=3)

la.print_mat(V)

Outputs:

>>> la.print_mat(V)
[1, 0, 0, 0]
[1, 1, 1, 1]
[1, 2, 4, 8]
[1, 3, 9, 27]
[1, 4, 16, 64]

Miscellaneous methods

Size of a matrix

The size of a matrix is simply its number of rows and columns.

Code implementation

Here we return the number of rows, or columns, or both depending on the axis argument:

def size(self, axis=2):
    if axis == 0:
        return len(self.data)
    elif axis == 1:
        return len(self.data[0])
    elif axis == 2:
        return [len(self.data), len(self.data[0])]

Transpose of a matrix

To transpose a matrix, we need a new matrix in which the rows are the former columns and the columns are the former rows:

\[\begin{equation} A =% \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{bmatrix}; A^T =% \begin{bmatrix} 1 & 4 \\ 2 & 5 \\ 3 & 6 \end{bmatrix} \end{equation}\]

Code implementation

We loop over the rows and columns of the input matrix, keeping a record of the row and column subscript indices. The transposed matrix is built up gradually with a set of new rows (each row being the values found in the corresponding column):

def transpose(self):
    transposed = []
    for i, row in enumerate(self.data):
        for j, col in enumerate(row):
            # first time through, make new row for each old column
            if i == 0:
                transposed.append([col])
            else:
                # append to newly created rows
                transposed[j].append(col)
    return Mat(transposed)

Demo

We create a matrix, call the transpose method and print the result:

import linalg as la

A = la.Mat([[1, 2, 3],
            [4, 5, 6]])

transposed = A.tr()

la.print_mat(transposed)

Outputs:

>>> la.print_mat(transposed)
[1, 4]
[2, 5]
[3, 6]

Indexing a matrix

While coding other functions in this library, I noticed I was often extracting a single row, column, or particular value from the matrix by doing:

my_value = A.data[0][0]

This seemed a bit clunky so I defined a method with a more readable interface to do it:

Code implementation

We check to see if a row index 'i' and if a column index 'j' were passed. If only a row index was passed, we return that entire row. If only a column index was passed, we return the entire column (by transposing the matrix to make indexing a column easy) and if both were passed we return the value in the 'ith' and 'jth' column (as an integer):

def ind(self, i=None, j=None):
    if isinstance(i, int) and not isinstance(j, int):
        return Mat([self.data[i]])
    elif isinstance(j, int) and not isinstance(i, int):
        return Mat([self.tr().data[j]]).tr()
    elif isinstance(i, int) and isinstance(j, int):
        return self.data[i][j]

Make 1x1 matrix a scalar

In case we end up with a 1x1 matrix, we probably want to use it as a scalar (that's all a 1x1 matrix is...) and so we make a method that indexes the element for us (since the intent is more understandable than indexing directly):

def make_scalar(self):
    if max(self.size()) == 1:
        return self.ind(0,0)

Get the diagonal of a matrix

This function works along the diagonal of the matrix, starting in the top left and stopping after running out of rows or columns.

def diag(self):
    diag_vals = []
    for idx in range(min(self.size())):
        diag_vals.append(self.ind(idx,idx))
    return diag_vals

Demo

We index the entire 2nd row, the entire 3rd column, or the value in row 2 column 3:

import linalg as la

A = la.Mat([[1, 2, 3],
            [-2, 1, 4],
            [0, 1, 2],
            [3, 6, 1]])

la.print_mat(A.ind(1,''))

la.print_mat(A.ind('', 2))

print(A.ind(1, 2))

Outputs:

 

>>> la.print_mat(A.ind(1,''))
[-2, 1, 4]

>>> la.print_mat(A.ind('', 2))
[3]
[4]
[2]
[1]

>>> print(A.ind(1, 2))
4

Is matrix triangular, diagonal, symmetric? >

back to project main page
back to home