0

I am trying to implement rabin's information dispersal algorithm which involves multiplication of a NxM matrix into a file. So I need to work with matrices and their inverses throughout my implementation. Along the way though I started testing generation of matrices and their inverses.

As far as I am concerned I know the result of multiplication of a matrix with its inverse must be an identiy matrix namely, its main diagonal being ones and having zeros elsewhere. So I tried generating a Vandermonde matrix as follows:

import numpy as np
from numpy.linalg import inv

x = np.array([1, 2, 3, 5])
x1 = np.vander(x)
x1
array([[  1,   1,   1,   1],
       [  8,   4,   2,   1],
       [ 27,   9,   3,   1],
       [125,  25,   5,   1]])

and calculated the inverse as follows:

x1_inv = inv(x1_inv)
inv(x1_inv)
array([[-0.125     ,  0.33333333, -0.25      ,  0.04166667],
       [ 1.25      , -3.        ,  2.        , -0.25      ],
       [-3.875     ,  7.66666667, -4.25      ,  0.45833333],
       [ 3.75      , -5.        ,  2.5       , -0.25      ]])

So far so good, but by by multiplying x1 and x1_inv theoretically I would be expecting an identity matrix like follows:

array([[ 1.,  0.,  0.,  0.],
   [ 0.,  1.,  0.,  0.],
   [ 0.,  0.,  1.,  0.],
   [ 0.,  0.,  0.,  1.]])

but here's what I get:

x1_inv.dot(x1)
array([[  1.00000000e+00,   0.00000000e+00,  -1.66533454e-16,
         -1.38777878e-17],
       [  3.55271368e-15,   1.00000000e+00,   4.44089210e-16,
          2.77555756e-17],
       [ -1.42108547e-14,  -1.77635684e-15,   1.00000000e+00,
         -8.32667268e-16],
       [ -7.10542736e-15,  -8.88178420e-16,   2.22044605e-16,
          1.00000000e+00]])

I'm getting ones on the main diagonal, fine, but not zeros elsewhere. What am I missing here? To my surprise though, when I multiply the supposedly yielded identity matrix with the main matrix I get the main matrix while I shouldn't since the yielded identity matrix actually isn't one.

and when I try np.allclose(x1.dot(inv(x1), np.eye(4)) I'm returned true.

My question being: What am I missing here? Why the result of x1.dot(x1_inv) doesn't seem what I would expect it to be (an identity matrix with ones on main diagonal and zeros elsewhere)?

Amir Afianian
  • 2,679
  • 4
  • 22
  • 46
  • 1
    these are the limitations of floating point calculations – Daniel Dec 27 '15 at 08:47
  • @Daniel is it confident that I pursue such way considering the fact that I am actually implementing some erasure coding? – Amir Afianian Dec 27 '15 at 08:49
  • 1
    erasure coding with floating points? – Daniel Dec 27 '15 at 08:51
  • @Daniel That's my point, I didn't expect x1.dot(x1_inv) not to be an identity matrix. By the way, Is it possible to move numpy matrix multiplications under a galois field or I have to implement the multiplications myself? – Amir Afianian Dec 27 '15 at 08:57
  • There is already a similar question. Maybe the answers can help: http://stackoverflow.com/questions/15707215/how-to-arrive-at-the-unit-matrix-from-numpy-dota-a-inv – LaPriWa Dec 27 '15 at 09:32
  • They might be not zeroes but they're quite damn close to it as it can be. – Loïc Faure-Lacroix Dec 27 '15 at 13:49

0 Answers0