Computing Determinant of Square Matrix using Condensation Method

determinant.py

#!/usr/bin/python3

import numpy as np
from copy import deepcopy
from functools import reduce
from time import time

N = 1 << 8

'''
Matrix of following form

[
    [1/1, 1/2, 1/3, 1/4], 
    [1/2, 1/3, 1/4, 1/5], 
    [1/3, 1/4, 1/5, 1/6], 
    [1/4, 1/5, 1/6, 1/7]
]
'''


def hilbert():
    mat = np.zeros((N, N))

    for i in range(N):
        for j in range(N):
            mat[i][j] = 1 / (i+j+1)

    return mat


# Multiplies successive pivots, obtained during condensation method


def mult_pivots(pivots):
    s1 = list(filter(lambda e: e >= -1. and e <= 1., pivots))
    s2 = list(filter(lambda e: not (e >= -1. and e <= 1.), pivots))

    while s1 and s2:
        q1 = s1.pop()
        q2 = s2.pop()
        q = q1 * q2
        if q >= -1. and q <= 1.:
            s1.append(q)
        else:
            s2.append(q)

    if s1:
        r = s1.copy()

        while len(r) > 1:
            l1 = reduce(lambda acc, cur: (cur[0], abs(cur[1] - 0.)) if acc[1] > abs(
                cur[1] - 0.) else acc, enumerate(r), (-1, float(1 << 32)))[0]
            l2 = reduce(lambda acc, cur: (cur[0], abs(cur[1] - 1.)) if cur[1] != 0 and acc[1] > abs(
                cur[1] - 1.) else acc, enumerate(r), (-1, float(1 << 32)))[0]

            r1 = r.pop(l1)
            r2 = 0
            if l1 == l2:
                r2 = r.pop()
            else:
                r2 = r.pop(l2-1 if l2 > l1 else l2)

            _r = r1 * r2
            r.append(_r)

        return r
    else:
        r = s2.copy()

        while len(r) > 1:
            r1 = r.pop()
            r2 = r.pop()

            _r = r1 * r2
            r.append(_r)

        return r

# Read: https://doi.org/10.1088/1742-6596%2F341%2F1%2F012031


def condense(mat):
    b = deepcopy(mat)
    pivots = []

    n = N
    for _ in range(N-2):
        l = reduce(
            lambda acc, cur: cur[0] if acc == -1 and cur[1] != 0 else acc, enumerate(b[0]), -1)
        if l == -1:
            break

        p = b[0][l]
        for i in range(n):
            b[i][l] /= p

        b_ = np.zeros((n-1, n-1))
        for i in range(n-1):
            for j in range(n-1):
                if j >= l:
                    b_[i][j] = b[0][l] * b[i+1][j+1] - b[0][j+1] * b[i+1][l]
                else:
                    b_[i][j] = -1 * b[i+1][j] * b[0][l]

        pivots.append(p)
        n -= 1
        b = b_

    if b.shape != (2, 2):
        return 0

    lst = b[0][0] * b[1][1] - b[0][1] * b[1][0]
    return lst * mult_pivots(pivots)[0]


if __name__ == '__main__':
    # a = hilbert()
    a = np.random.random((N, N))
    start = time() * (10 ** 3)
    print(f"determinant: {condense(a)}\t[computed]\tin {time() * (10 ** 3) - start} ms")
    start = time() * (10 ** 3)
    print(f"determinant: {np.linalg.det(a)}\t[numpy]\tin {time() * (10 ** 3) - start} ms")