diff --git a/matrix/matrix_operation.py b/matrix/matrix_operation.py
index dd7c01582681..b32a4dcf7af3 100644
--- a/matrix/matrix_operation.py
+++ b/matrix/matrix_operation.py
@@ -1,64 +1,131 @@
-from __future__ import print_function
+"""
+function based version of matrix operations, which are just 2D arrays
+"""
+
 
 def add(matrix_a, matrix_b):
-    rows = len(matrix_a)
-    columns = len(matrix_a[0])
-    matrix_c = []
-    for i in range(rows):
-        list_1 = []
-        for j in range(columns):
-            val = matrix_a[i][j] + matrix_b[i][j]
-            list_1.append(val)
-        matrix_c.append(list_1)
-    return matrix_c
-
-def scalarMultiply(matrix , n):
+    if _check_not_integer(matrix_a) and _check_not_integer(matrix_b):
+        rows, cols = _verify_matrix_sizes(matrix_a, matrix_b)
+        matrix_c = []
+        for i in range(rows[0]):
+            list_1 = []
+            for j in range(cols[0]):
+                val = matrix_a[i][j] + matrix_b[i][j]
+                list_1.append(val)
+            matrix_c.append(list_1)
+        return matrix_c
+
+
+def subtract(matrix_a, matrix_b):
+    if _check_not_integer(matrix_a) and _check_not_integer(matrix_b):
+        rows, cols = _verify_matrix_sizes(matrix_a, matrix_b)
+        matrix_c = []
+        for i in range(rows[0]):
+            list_1 = []
+            for j in range(cols[0]):
+                val = matrix_a[i][j] - matrix_b[i][j]
+                list_1.append(val)
+            matrix_c.append(list_1)
+        return matrix_c
+
+
+def scalar_multiply(matrix, n):
     return [[x * n for x in row] for row in matrix]
 
+
 def multiply(matrix_a, matrix_b):
-    matrix_c = []
-    n = len(matrix_a)
-    for i in range(n):
-        list_1 = []
-        for j in range(n):
-            val = 0
-            for k in range(n):
-                val = val + matrix_a[i][k] * matrix_b[k][j]
-            list_1.append(val)
-        matrix_c.append(list_1)
-    return matrix_c
+    if _check_not_integer(matrix_a) and _check_not_integer(matrix_b):
+        matrix_c = []
+        rows, cols = _verify_matrix_sizes(matrix_a, matrix_b)
+
+        if cols[0] != rows[1]:
+            raise ValueError(f'Cannot multiply matrix of dimensions ({rows[0]},{cols[0]}) '
+                             f'and ({rows[1]},{cols[1]})')
+        for i in range(rows[0]):
+            list_1 = []
+            for j in range(cols[1]):
+                val = 0
+                for k in range(cols[1]):
+                    val = val + matrix_a[i][k] * matrix_b[k][j]
+                list_1.append(val)
+            matrix_c.append(list_1)
+        return matrix_c
+
 
 def identity(n):
+    """
+    :param n: dimension for nxn matrix
+    :type n: int
+    :return: Identity matrix of shape [n, n]
+    """
+    n = int(n)
     return [[int(row == column) for column in range(n)] for row in range(n)] 
 
-def transpose(matrix):
-    return map(list , zip(*matrix))
+
+def transpose(matrix, return_map=True):
+    if _check_not_integer(matrix):
+        if return_map:
+            return map(list, zip(*matrix))
+        else:
+            # mt = []
+            # for i in range(len(matrix[0])):
+            #     mt.append([row[i] for row in matrix])
+            # return mt
+            return [[row[i] for row in matrix] for i in range(len(matrix[0]))]
+
 
 def minor(matrix, row, column):
     minor = matrix[:row] + matrix[row + 1:]
     minor = [row[:column] + row[column + 1:] for row in minor]
     return minor
 
+
 def determinant(matrix):
-    if len(matrix) == 1: return matrix[0][0]
+    if len(matrix) == 1:
+        return matrix[0][0]
     
     res = 0
     for x in range(len(matrix)):
-        res += matrix[0][x] * determinant(minor(matrix , 0 , x)) * (-1) ** x
+        res += matrix[0][x] * determinant(minor(matrix, 0, x)) * (-1) ** x
     return res
 
+
 def inverse(matrix):
     det = determinant(matrix)
-    if det == 0: return None
+    if det == 0:
+        return None
 
-    matrixMinor = [[] for _ in range(len(matrix))]
+    matrix_minor = [[] for _ in range(len(matrix))]
     for i in range(len(matrix)):
         for j in range(len(matrix)):
-            matrixMinor[i].append(determinant(minor(matrix , i , j)))
+            matrix_minor[i].append(determinant(minor(matrix, i, j)))
     
-    cofactors = [[x * (-1) ** (row + col) for col, x in enumerate(matrixMinor[row])] for row in range(len(matrix))]
+    cofactors = [[x * (-1) ** (row + col) for col, x in enumerate(matrix_minor[row])] for row in range(len(matrix))]
     adjugate = transpose(cofactors)
-    return scalarMultiply(adjugate , 1/det)
+    return scalar_multiply(adjugate, 1/det)
+
+
+def _check_not_integer(matrix):
+    try:
+        rows = len(matrix)
+        cols = len(matrix[0])
+        return True
+    except TypeError:
+        raise TypeError("Cannot input an integer value, it must be a matrix")
+
+
+def _shape(matrix):
+    return list((len(matrix), len(matrix[0])))
+
+
+def _verify_matrix_sizes(matrix_a, matrix_b):
+    shape = _shape(matrix_a)
+    shape += _shape(matrix_b)
+    if shape[0] != shape[2] or shape[1] != shape[3]:
+        raise ValueError(f"operands could not be broadcast together with shape "
+                         f"({shape[0], shape[1]}), ({shape[2], shape[3]})")
+    return [shape[0], shape[2]], [shape[1], shape[3]]
+
 
 def main():
     matrix_a = [[12, 10], [3, 9]]
@@ -68,9 +135,10 @@ def main():
     print('Add Operation, %s + %s = %s \n' %(matrix_a, matrix_b, (add(matrix_a, matrix_b))))
     print('Multiply Operation, %s * %s = %s \n' %(matrix_a, matrix_b, multiply(matrix_a, matrix_b)))
     print('Identity:  %s \n' %identity(5))
-    print('Minor of %s = %s \n' %(matrix_c, minor(matrix_c , 1 , 2)))
+    print('Minor of %s = %s \n' %(matrix_c, minor(matrix_c, 1, 2)))
     print('Determinant of %s = %s \n' %(matrix_b, determinant(matrix_b)))
     print('Inverse of %s = %s\n'%(matrix_d, inverse(matrix_d)))
 
+
 if __name__ == '__main__':
     main()
diff --git a/matrix/tests/pytest.ini b/matrix/tests/pytest.ini
new file mode 100644
index 000000000000..8a978b56ef8b
--- /dev/null
+++ b/matrix/tests/pytest.ini
@@ -0,0 +1,3 @@
+[pytest]
+markers =
+    mat_ops: tests for matrix operations
diff --git a/matrix/tests/test_matrix_operation.py b/matrix/tests/test_matrix_operation.py
new file mode 100644
index 000000000000..8b81b65d0fc8
--- /dev/null
+++ b/matrix/tests/test_matrix_operation.py
@@ -0,0 +1,112 @@
+"""
+Testing here assumes that numpy and linalg is ALWAYS correct!!!!
+
+If running from PyCharm you can place the following line in "Additional Arguments" for the pytest run configuration
+-vv -m mat_ops -p no:cacheprovider
+"""
+
+# standard libraries
+import sys
+import numpy as np
+import pytest
+import logging
+
+# Custom/local libraries
+from matrix import matrix_operation as matop
+
+mat_a = [[12, 10], [3, 9]]
+mat_b = [[3, 4], [7, 4]]
+mat_c = [[3, 0, 2], [2, 0, -2], [0, 1, 1]]
+mat_d = [[3, 0, -2], [2, 0, 2], [0, 1, 1]]
+mat_e = [[3, 0, 2], [2, 0, -2], [0, 1, 1], [2, 0, -2]]
+mat_f = [1]
+mat_h = [2]
+
+logger = logging.getLogger()
+logger.level = logging.DEBUG
+stream_handler = logging.StreamHandler(sys.stdout)
+logger.addHandler(stream_handler)
+
+
+@pytest.mark.mat_ops
+@pytest.mark.parametrize(('mat1', 'mat2'), [(mat_a, mat_b), (mat_c, mat_d), (mat_d, mat_e),
+                                            (mat_f, mat_h)])
+def test_addition(mat1, mat2):
+    if (np.array(mat1)).shape < (2, 2) or (np.array(mat2)).shape < (2, 2):
+        with pytest.raises(TypeError):
+            logger.info(f"\n\t{test_addition.__name__} returned integer")
+            matop.add(mat1, mat2)
+    elif (np.array(mat1)).shape == (np.array(mat2)).shape:
+        logger.info(f"\n\t{test_addition.__name__} with same matrix dims")
+        act = (np.array(mat1) + np.array(mat2)).tolist()
+        theo = matop.add(mat1, mat2)
+        assert theo == act
+    else:
+        with pytest.raises(ValueError):
+            logger.info(f"\n\t{test_addition.__name__} with different matrix dims")
+            matop.add(mat1, mat2)
+
+
+@pytest.mark.mat_ops
+@pytest.mark.parametrize(('mat1', 'mat2'), [(mat_a, mat_b), (mat_c, mat_d), (mat_d, mat_e),
+                                            (mat_f, mat_h)])
+def test_subtraction(mat1, mat2):
+    if (np.array(mat1)).shape < (2, 2) or (np.array(mat2)).shape < (2, 2):
+        with pytest.raises(TypeError):
+            logger.info(f"\n\t{test_subtraction.__name__} returned integer")
+            matop.subtract(mat1, mat2)
+    elif (np.array(mat1)).shape == (np.array(mat2)).shape:
+        logger.info(f"\n\t{test_subtraction.__name__} with same matrix dims")
+        act = (np.array(mat1) - np.array(mat2)).tolist()
+        theo = matop.subtract(mat1, mat2)
+        assert theo == act
+    else:
+        with pytest.raises(ValueError):
+            logger.info(f"\n\t{test_subtraction.__name__} with different matrix dims")
+            assert matop.subtract(mat1, mat2)
+
+
+@pytest.mark.mat_ops
+@pytest.mark.parametrize(('mat1', 'mat2'), [(mat_a, mat_b), (mat_c, mat_d), (mat_d, mat_e),
+                                            (mat_f, mat_h)])
+def test_multiplication(mat1, mat2):
+    if (np.array(mat1)).shape < (2, 2) or (np.array(mat2)).shape < (2, 2):
+        logger.info(f"\n\t{test_multiplication.__name__} returned integer")
+        with pytest.raises(TypeError):
+            matop.add(mat1, mat2)
+    elif (np.array(mat1)).shape == (np.array(mat2)).shape:
+        logger.info(f"\n\t{test_multiplication.__name__} meets dim requirements")
+        act = (np.matmul(mat1, mat2)).tolist()
+        theo = matop.multiply(mat1, mat2)
+        assert theo == act
+    else:
+        with pytest.raises(ValueError):
+            logger.info(f"\n\t{test_multiplication.__name__} does not meet dim requirements")
+            assert matop.subtract(mat1, mat2)
+
+
+@pytest.mark.mat_ops
+def test_scalar_multiply():
+    act = (3.5 * np.array(mat_a)).tolist()
+    theo = matop.scalar_multiply(mat_a, 3.5)
+    assert theo == act
+
+
+@pytest.mark.mat_ops
+def test_identity():
+    act = (np.identity(5)).tolist()
+    theo = matop.identity(5)
+    assert theo == act
+
+
+@pytest.mark.mat_ops
+@pytest.mark.parametrize('mat', [mat_a, mat_b, mat_c, mat_d, mat_e, mat_f])
+def test_transpose(mat):
+    if (np.array(mat)).shape < (2, 2):
+        with pytest.raises(TypeError):
+            logger.info(f"\n\t{test_transpose.__name__} returned integer")
+            matop.transpose(mat)
+    else:
+        act = (np.transpose(mat)).tolist()
+        theo = matop.transpose(mat, return_map=False)
+        assert theo == act