How can this assert() ever trigger?

[Albert van der Horst]

I have the following code for calculating the determinant of
a matrix. It works inasfar that it gives the same result as an
octave program on a same matrix.

/ ----------------------------------------------------------------

def determinant( mat ):
''' Return the determinant of the n by n matrix mat
i row j column
Destroys mat ! '''
#print "getting determinat of", mat
n=len(mat)
nom = 1.
if n == 1: return mat[0][0]
lastr = mat.pop()
jx=-1
for j in xrange(n):
if lastr[j]:
jx=j
break
if jx==-1: return 0.
result = lastr[jx]
assert(result<>0.)
# Make column jx zero by subtracting a multiple of the last row.
for i in xrange(n-1):
pivot = mat[jx]
if 0. == pivot: continue
assert(result<>0.)
nom *= result # Compenstate for multiplying a row.
for j in xrange(n):
mat[j] *= result
for j in xrange(n):
mat[j] -= pivot*lastr[j]
# Remove colunm jx
for i in xrange(n-1):
x= mat.pop(jx)
assert( x==0 )

if (n-1+jx)%2<>0: result = -result
det = determinant( mat )
assert(nom<>0.)
return result*det/nom

/-----------------------------------------

Now on some matrices the assert triggers, meaning that nom is zero.
How can that ever happen? mon start out as 1. and gets multiplied
with a number that is asserted to be not zero.

Any hints appreciated.

Groetjes Albert

 

[Ethan Furman]

What happens if you run the same matrix through Octave? By any chance, is nom just really, really small?

 

[Gary Herron]

I have the following code for calculating the determinant of
a matrix. It works inasfar that it gives the same result as an
octave program on a same matrix.

/ ----------------------------------------------------------------

def determinant( mat ):
''' Return the determinant of the n by n matrix mat
i row j column
Destroys mat ! '''
#print "getting determinat of", mat
n=len(mat)
nom = 1.
if n == 1: return mat[0][0]
lastr = mat.pop()
jx=-1
for j in xrange(n):
if lastr[j]:
jx=j
break
if jx==-1: return 0.
result = lastr[jx]
assert(result<>0.)
# Make column jx zero by subtracting a multiple of the last row.
for i in xrange(n-1):
pivot = mat[jx]
if 0. == pivot: continue
assert(result<>0.)
nom *= result # Compenstate for multiplying a row.
for j in xrange(n):
mat[j] *= result
for j in xrange(n):
mat[j] -= pivot*lastr[j]
# Remove colunm jx
for i in xrange(n-1):
x= mat.pop(jx)
assert( x==0 )

if (n-1+jx)%2<>0: result = -result
det = determinant( mat )
assert(nom<>0.)
return result*det/nom

/-----------------------------------------

Now on some matrices the assert triggers, meaning that nom is zero.
How can that ever happen? mon start out as 1. and gets multiplied
with a number that is asserted to be not zero.

Easily due to *underflow* precision trouble. Your "result" may never be
zero, but it can be very small. Take the product of many of such tiny
values, and the result can be less then the smallest value representable
by a float, at which point it becomes zero.

To see this clearly, try this Python code:.... a = a*1.0e-50
.... print(a)
....
1e-50
1e-100
1e-150
1e-200
1e-250
1e-300
0.0

Gary Herron

 

[Peter Otten]

I have the following code for calculating the determinant of
a matrix. It works inasfar that it gives the same result as an
octave program on a same matrix.

/ ----------------------------------------------------------------

def determinant( mat ):
''' Return the determinant of the n by n matrix mat
i row j column
Destroys mat ! '''
#print "getting determinat of", mat
n=len(mat)
nom = 1.
if n == 1: return mat[0][0]
lastr = mat.pop()
jx=-1
for j in xrange(n):
if lastr[j]:
jx=j
break
if jx==-1: return 0.
result = lastr[jx]
assert(result<>0.)
# Make column jx zero by subtracting a multiple of the last row.
for i in xrange(n-1):
pivot = mat[jx]
if 0. == pivot: continue
assert(result<>0.)
nom *= result # Compenstate for multiplying a row.
for j in xrange(n):
mat[j] *= result
for j in xrange(n):
mat[j] -= pivot*lastr[j]
# Remove colunm jx
for i in xrange(n-1):
x= mat.pop(jx)
assert( x==0 )

if (n-1+jx)%2<>0: result = -result
det = determinant( mat )
assert(nom<>0.)
return result*det/nom

/-----------------------------------------

Now on some matrices the assert triggers, meaning that nom is zero.
How can that ever happen? mon start out as 1. and gets multiplied
with a number that is asserted to be not zero.

Any hints appreciated.

Floating point precision is limited:
.... x *= .1
.... assert x
....
Traceback (most recent call last):
323

 

[Alain Ketterlin]

(e-mail address removed)4all.nl (Albert van der Horst) writes:

[...]

Now on some matrices the assert triggers, meaning that nom is zero.
How can that ever happen? mon start out as 1. and gets multiplied

[several times]

with a number that is asserted to be not zero.

Finite precision. Try: 1.*1e-162*1e-162. Equals zero.

-- Alain.

 

[Albert van der Horst]

(e-mail address removed)4all.nl (Albert van der Horst) writes:

[...]

Now on some matrices the assert triggers, meaning that nom is zero.
How can that ever happen? mon start out as 1. and gets multiplied

[several times]

with a number that is asserted to be not zero.

Finite precision. Try: 1.*1e-162*1e-162. Equals zero.

-- Alain.

Thanks you Alan and all others to point this out.
That was indeed the problem. Somehow I expected that floating
underflow would raise an exception, so I had a blind spot there.

Groetjes Albert

 

[Joseph Martinot-Lagarde]

Le 10/05/2014 17:24, Albert van der Horst a écrit :

I have the following code for calculating the determinant of
a matrix. It works inasfar that it gives the same result as an
octave program on a same matrix.

/ ----------------------------------------------------------------

def determinant( mat ):
''' Return the determinant of the n by n matrix mat
i row j column
Destroys mat ! '''
#print "getting determinat of", mat
n=len(mat)
nom = 1.
if n == 1: return mat[0][0]
lastr = mat.pop()
jx=-1
for j in xrange(n):
if lastr[j]:
jx=j
break
if jx==-1: return 0.
result = lastr[jx]
assert(result<>0.)
# Make column jx zero by subtracting a multiple of the last row.
for i in xrange(n-1):
pivot = mat[jx]
if 0. == pivot: continue
assert(result<>0.)
nom *= result # Compenstate for multiplying a row.
for j in xrange(n):
mat[j] *= result
for j in xrange(n):
mat[j] -= pivot*lastr[j]
# Remove colunm jx
for i in xrange(n-1):
x= mat.pop(jx)
assert( x==0 )

if (n-1+jx)%2<>0: result = -result
det = determinant( mat )
assert(nom<>0.)
return result*det/nom

/-----------------------------------------

Now on some matrices the assert triggers, meaning that nom is zero.
How can that ever happen? mon start out as 1. and gets multiplied
with a number that is asserted to be not zero.

Any hints appreciated.

Groetjes Albert

I know it's not the question, but if you want a replacement for octave
did you try numpy (and scipy) ? The determinant would be computer faster
and with less memory than with your function.

 

[Albert van der Horst]

Le 10/05/2014 17:24, Albert van der Horst a écrit :

I have the following code for calculating the determinant of
a matrix. It works inasfar that it gives the same result as an
octave program on a same matrix.

/ ----------------------------------------------------------------

def determinant( mat ): ...
result = lastr[jx]
assert(result<>0.) ....
assert(result<>0.)
nom *= result # Compenstate for multiplying a row. ....
assert(nom<>0.) ...

/-----------------------------------------

Now on some matrices the assert triggers, meaning that nom is zero.
How can that ever happen? mon start out as 1. and gets multiplied
with a number that is asserted to be not zero.

Any hints appreciated.

Groetjes Albert

I know it's not the question, but if you want a replacement for octave
did you try numpy (and scipy) ? The determinant would be computer faster
and with less memory than with your function.

I'm using several programming languages in a mix to solve Euler problems.
This is about learning how octave compares to python for a certain kind of
problem as anything.
The determinant program I had lying around, but it was infinite precision
with integer only arithmetic. Then I made a simple modification and got
mad because I didn't understand why it didn't work.

I have used numpy and its det before, but I find it difficult to
remember how to create a matrix in numpy. This is the kind of thing
that is hard to find in the docs. Now I looked it up in my old
programs: you start a matrix with the zeroes() function.

I expect the built in determinant of octave to be on a par with corresponding
python libraries.

---

Groetjes Albert

 

[Joseph Martinot-Lagarde]

Le 13/05/2014 11:56, Albert van der Horst a écrit :

Le 10/05/2014 17:24, Albert van der Horst a écrit :

I have the following code for calculating the determinant of
a matrix. It works inasfar that it gives the same result as an
octave program on a same matrix.

/ ----------------------------------------------------------------

def determinant( mat ): ..
result = lastr[jx]
assert(result<>0.) ...
assert(result<>0.)
nom *= result # Compenstate for multiplying a row. ...
assert(nom<>0.) ..

/-----------------------------------------

Now on some matrices the assert triggers, meaning that nom is zero.
How can that ever happen? mon start out as 1. and gets multiplied
with a number that is asserted to be not zero.

Any hints appreciated.

Groetjes Albert

I know it's not the question, but if you want a replacement for octave
did you try numpy (and scipy) ? The determinant would be computer faster
and with less memory than with your function.

I'm using several programming languages in a mix to solve Euler problems.
This is about learning how octave compares to python for a certain kind of
problem as anything.
The determinant program I had lying around, but it was infinite precision
with integer only arithmetic. Then I made a simple modification and got
mad because I didn't understand why it didn't work.

I have used numpy and its det before, but I find it difficult to
remember how to create a matrix in numpy. This is the kind of thing
that is hard to find in the docs. Now I looked it up in my old
programs: you start a matrix with the zeroes() function.

I expect the built in determinant of octave to be on a par with corresponding
python libraries.

---

Groetjes Albert

You can use numpy.zeros(), but you can also use the same list of lists
that you use for your problem.

Transform a list of lists into a numpy array:

np.asarray([[1, 2],[3, 4]])

array([[1, 2],
[3, 4]])

Use a numpy function directly on a list of lists (works for must numpy
functions):

np.linalg.det([[1, 2],[3, 4]])

-2.0000000000000004

More info on array creation:
http://wiki.scipy.org/Tentative_NumPy_Tutorial#head-d3f8e5fe9b903f3c3b2a5c0dfceb60d71602cf93