Getting the Determinant of a Matrix in PyRTL

11 Jun 2021

In my recent kick to understand matrices a little better, I decided to implement a function to get the determinant of an n by n matrix.

Here's the outline of what we'll be covering:

Pure Python
PyRTL Version
LU Decomposition

Pure Python

Like usual, I started by making a pure Python implementation:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34 def determinant(matrix):
    ''' Get the determinant of the matrix.

    :param list[list[int]] matrix: Matrix to operate on.
    :return: An integer determinant.

    Note: This only works for square matrices.
    '''
    from functools import reduce
    from itertools import permutations

    if len(matrix) != len(matrix[0]):
        raise Exception("Must be square")

    def perms(n):
        return list(permutations(range(n)))

    def sgn(l):
        count = 0
        for i in range(len(l)):
            for j in range(i + 1, len(l)):
                if l[i] > l[j]:
                    count += 1
        return 1 if (count % 2) == 0 else -1

    res = 0
    n = len(matrix)
    for perm in perms(n):
        s = sgn(perm)
        p = reduce(lambda a1, a2: a1 * a2,
                   [matrix[i][perm[i]] for i in range(n)])
        res += s * p

    return res

The above approach uses the Liebniz formula, which has a complexity of O(n!) (not great…). But as long as we’re testing with relatively modest-sized matrices, we’ll be okay.

Let’s test it out:

>>> determinant([[6, 1, 1], [4, -2, 5], [2, 8, 7]])
-306

(That is correct, by the way).

PyRTL Version

So let’s make it in PyRTL!

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36 from functools import reduce
from itertools import permutations
import pyrtl

def det(matrix):
    ''' Get the determinant of the matrix.

    :param Matrix matrix: Matrix to operate on.
    :return: A WireVector containing the determinant.

    Note: This only works for square matrices.
    '''

    if matrix.rows != matrix.columns:
        raise pyrtl.PyrtlError('Matrix must be square.')

    def perms(n):
        return list(permutations(range(n)))

    def sgn(ls):
        count = 0
        for i in range(len(ls)):
            for j in range(i + 1, len(ls)):
                if ls[i] > ls[j]:
                    count += 1
        return 1 if (count % 2) == 0 else -1

    res = 0
    n = matrix.rows
    for perm in perms(n):
        s = sgn(perm)
        p = reduce(lambda a1, a2: pyrtl.signed_mult(a1, a2),
                   [matrix[i, perm[i]] for i in range(n)])
        res = pyrtl.signed_add(res, pyrtl.signed_mult(s, p))

    return res

As you can probably tell, very little needed to change. The main changes are in lines 32 and 33. Essentially, we’ve replaced the primitive addition and multiplication calls with PyRTL’s signed_add and signed_mult. In particular, note that the return value from sgn is an integer (possibly negative). By passing it to pyrtl.signed_mult on line 32, PyRTL automatically converts the integer into a PyRTL Const with signed=True — very convenient!

Now let’s test this:

>>> import pyrtl.rtllib.matrix as Matrix
>>>
>>> m = Matrix.Matrix(3, 3, bits=5, signed=True, value=[[6, 1, 1], [4, -2, 5], [2, 8, 7]])
>>> result = Matrix.det(m)
>>> pyrtl.probe(result, 'result')
>>>
>>> sim = pyrtl.Simulation()
>>> sim.step({})
>>>
>>> expected = pyrtl.formatted_str_to_val('-306', 's' + str(result.bitwidth))
>>> assert(sim.inspect('result') == expected)

Since the value returned from inspecting a wire in simulation is not interpreted as signed, we first called formatted_str_to_val to convert our -306 expected value into the two’s complement equivalent for comparison (we could also have converted the inspected value instead). And it works!

LU Decomposition

Like I mentioned previously, the method used was Leibniz’s formula, which isn’t terribly efficient. It would be better to use LU decomposition to compute the determinant, but we’ll leave that for another day.