Test if matrix is invertible over finite field

Question

I would like to test if a particular type of random matrix is invertible over a finite field, in particular F_2. I can test if a matrix is invertible over the reals using the following simple code.

import random
from scipy.linalg import toeplitz
import numpy as np
n=10
column = [random.choice([0,1]) for x in xrange(n)]
row = [column[0]]+[random.choice([0,1]) for x in xrange(n-1)]
matrix = toeplitz(column, row)
if (np.linalg.matrix_rank(matrix) < n):
    print "Not invertible!"

Is there some way to achieve the same thing but over F_2?

You could do it with Sage easily enough ([example](http://aleph.sagemath.org/?z=eJzzDVawVfBNLCnKrAguSExO1XB30zDS1FEwBiJNXq7yjMycVIWQotJUK14uBSDwBSkP1itKzEvJz41PzUnNTc0r0dCESGamKfjqZRbHZ-aVpRaVZCblpGpoQvWBQFJRamI2gsvLVVCUmVeioO5rpQ5j-yIEgYYgieuBzSxOBVkFU6GFpkZBC1UdABH6PRM=&lang=sage)). I'll be interested to see if there's a slick solution on the science stack (numpy/scipy/sympy/mpmath/pandas etc.), though. — DSM, Apr 27 '13 at 17:09
I think that if you regard the matrix over F_2 as a matrix over Z using only 0 and 1, then the determinant over F_2 should be the determinant over Z modulo 2 (i.e. the check becomes if the determinant over Z is even or odd). This might not be algorithmically optimal. — Armin Rigo, Apr 27 '13 at 19:01
@ArminRigo Unfortunately I can't get this idea to work. Set n =100 in the code above and print linalg.det(matrix), linalg.det(matrix)%2. I always get 0 for linalg.det(matrix)%2 which is presumably because of floating point problems. Is there an exact integer determinant function? — marshall, Apr 27 '13 at 22:36
@marshall: I noticed you added the `sympy` tag. It'll do the full integer determinant, but it'll be awfully slow. — DSM, Apr 27 '13 at 22:51
@DSM Yes, I don't really want to use sympy for that. Ultimately I wonder which part of the science stack *should* have this in it? — marshall, Apr 28 '13 at 06:13

pv. · Accepted Answer · 2013-04-27T23:35:30.443

It would be better to use Sage or some other proper tool for this.

The following is just unsophisticated non-expert attempt at doing something, but pivoted Gaussian elimination should give the exact result for invertibility:

import random
from scipy.linalg import toeplitz
import numpy as np

def is_invertible_F2(a):
    """
    Determine invertibility by Gaussian elimination
    """
    a = np.array(a, dtype=np.bool_)
    n = a.shape[0]
    for i in range(n):
        pivots = np.where(a[i:,i])[0]
        if len(pivots) == 0:
            return False

        # swap pivot
        piv = i + pivots[0]
        row = a[piv,i:].copy()
        a[piv,i:] = a[i,i:]
        a[i,i:] = row

        # eliminate
        a[i+1:,i:] -= a[i+1:,i,None]*row[None,:]

    return True

n = 10
column = [random.choice([0,1]) for x in xrange(n)]
row = [column[0]]+[random.choice([0,1]) for x in xrange(n-1)]
matrix = toeplitz(column, row)

print(is_invertible_F2(matrix))
print(int(np.round(np.linalg.det(matrix))) % 2)

Note that np.bool_ is analogous to F_2 only in a restricted sense --- the binary operation + in F_2 is - for bool, and the unary op - is +. Multiplication is the same, though.

>>> x = np.array([0, 1], dtype=np.bool_)
>>> x[:,None] - x[None,:]
array([[False,  True],
       [ True, False]], dtype=bool)
>>> x[:,None] * x[None,:]
array([[False, False],
       [False,  True]], dtype=bool)

The gaussian elimination above uses only these operations, so it works.

Thank you. I don't mind importing external libraries for this specific task if that is the right thing to do. I have never used sage and have no idea how well it interacts with scipy matrices for example. — marshall, Apr 28 '13 at 09:34

score 1 · Answer 2 · answered Apr 28 '13 at 20:01

Unfortunately, SymPy can't yet handle finite fields in matrices, though support is planned.

As some commenters noted, though, you can just check the determinant over the integers. If it's 1 (mod 2), the matrix is invertible. To actually find the inverse, you can just take the normal inverse over the integers, multiply by the determinant (so that you don't have fractions), and mod each element by 2. I can't imagine it would be too efficient, and you could probably use any matrix library, even a numerical one (rounding to the nearest integer). SymPy can also do each of these steps.

I should point out that in general cyclic finite fields, the "multiply by the determinant" part will need to be undone by multiplying by the inverse mod p (it is unnecessary mod 2 because the only possibility is 1).

Thanks that's interesting. I suppose a nice addition would be scipy could compute the determinant over the integers. — marshall, Apr 29 '13 at 07:46

Test if matrix is invertible over finite field

2 Answers2

Linked