Matrix inversion tool

[Matan Ziv-Av]

On Mon, 10 Aug 2015, Oleg Goldshmidt wrote:

> I assure you I read and understood Omer's and your posts. If you go back
> to my reply you will surely realize that at no point I contradicted
> you. I just pointed out that you didn't need to divide (by det(A)==1),
> which would lead to the problem you correctly point out.

Here's what Omer actually asked:

> What happens if you use the regular matrix inversion tool which works 
> on real numbers?

And I explained that a regular matrix inversion tool will return a 
result that cannot be converted back to Z_2. If you are using another 
tool, you do not have the option to "not divide by the determinant", 
since that too returns the actual inverse matrix after the "division by 
the determinant" (though, it never actually divides by the determinant).

> Friends again? ;-)

Obviously, I have no personal problems with you. I just don't like 
incorrect information (about a subject I really care about) available, 
since it has a way of returning very unexpectedly.

To show that there is no animosity, I'll answer your other questions:

> 273x273 seems quite
> large for a regular computer. It is not clear to me (I'll appreciate a
> comment just for curiousity's sake as well as for education) that a
> straightforward algorithm will be practical.

At the end of this message is a python program with simple 
implementations of both algorithms (Gauss eliminiation and recursive). 
Run them for sizes 10 to 16 consecuitively to see how the difference 
between exponential and polynomial is very significant in this case.

Then run for 273 (or 2730, for that matter), to see that with the right 
algorithms, 273x273 is very small for a computer.

If you save it as x.py, run as
./x.py n k s
Where n is the number of rows, k is the expected number of non-zeros in 
a row, and the optional s is a seed for the pseudo-random number 
generator.

The output is the determinant by GE and the time it takes, and then the 
determinant by the recursive method and then the time it takes.

The optimizations you suggested are possible (for both algorithms), but 
they are all just scaling the time needed. Even if you speed up by a 
factor of 10^9, the recursive algorithm will just not work.

The sparsity of the matrix does not help much. If there are (on average) 
k non-zero entries per line, you get O(k^n) instead of O(n!).

> [*] I strongly suspect you read too much into Omer's use of the word
>    "real" (also used in the OP).

Actually he said real numbers. How can you read it in any other way? 
especially in this context.

>
>

-- 
Matan.

#!/usr/bin/python

import random
import sys
import time
import copy

def mat(n,p):
     return [[int(random.random()<p) for e in range(n)] for e in range(n)]

def id(n):
     return [[int(i==j) for i in range(n)] for j in range(n)]

def ExchangeLines(m,i,j):
     m[i],m[j] = m[j],m[i]

def AddLines(m,i,j):
     n=len(m)
     m[i] = [ (m[i][k]+m[j][k])%2 for k in range(n)]

def ge(a,b):
     n=len(a)
     for i in range(n):
         for j in range(i,n):
             if a[j][i]!=0:
                 ExchangeLines(a,i,j)
                 ExchangeLines(b,i,j)
                 break
         if a[i][i]==0:
             return 0
         for j in range(n):
             if j!=i:
                 if a[j][i]!=0:
                     AddLines(a,j,i)
                     AddLines(b,j,i)
     return 1

def minor(a,i,j):
     b=copy.deepcopy(a)
     b.pop(i)
     for x in b:
         x.pop(j)
     return b

def recdet(a):
     n=len(a)
     if n==1:
         return a[0][0]
     s=0
     for i in range(n):
         if a[0][i]!=0:
             s=(s+recdet(minor(a,0,i)))%2
     return s


n=int(sys.argv[1])
p=float(sys.argv[2])
if len(sys.argv)>=4:
     random.seed(int(sys.argv[3]))

a1=mat(n,p/n)
a2=copy.deepcopy(a1)
b=id(n)

s=time.clock()
print "d1=",ge(a1,b)
e=time.clock()
print e-s

s=time.clock()
print "d2=",recdet(a2)
e=time.clock()
print e-s