4

I want to perform operations on rational matrices. I use the modules numpy and fractions.

Here is my code:

import numpy as np
from fractions import Fraction

m=np.matrix([[Fraction(1, 6), Fraction(8, 7)], [Fraction(1, 2), Fraction(3, 2)]])
print(np.linalg.det(m))
# Gives -0.321428571429

print(m[0,0]*m[1,1] - m[0,1]*m[1,0])
# Gives -9/28

Since computing the determinant only require rational operations with the Gauss' method, the determinant of a rational matrix is rational.

So my questions are: why does numpy return a float and not a Fraction? How can I get a rational determinant?

Note that other operations on this matrix give a rational output (for instance m.trace()).

Tom Cornebize
  • 1,362
  • 15
  • 33

2 Answers2

5

NumPy computes the determinant of the matrix by a lower upper decomposition routine in LAPACK. This routine can only handle floating point numbers.

Before calculating the determinant of the matrix, linalg.det checks the types of values it has and then establishes the type of internal loop that should be run using a call to a function named _commonType(). This function will set the loop to run for either double or complex-double values.

Here is the Python part of the function linalg.det that handles the checking:

def det(a):
    a = asarray(a) # convert matrix to NumPy array
    _assertNoEmpty2d(a)
    _assertRankAtLeast2(a)
    _assertNdSquareness(a)
    t, result_t = _commonType(a) # input/output types established here
    signature = 'D->D' if isComplexType(t) else 'd->d' # signature 'float->float' chosen 
    return _umath_linalg.det(a, signature=signature).astype(result_t)

After running checks on the shape of the matrix and determining types, the return line passes the values in the array to the LAPACK implementation of the lower-upper decomposition and a float is returned.

Trying to bypass this type checking with a type signature of our own raises an error saying that no such loop is defined for object types:

>>> np.linalg._umath_linalg.det(a, signature='O->O') # 'O' is 'object'
TypeError: No loop matching the specified signature was found for ufunc det

This implies than it is not possible to keep the Fraction type as the return type when using det.

Other functions such as trace() do not do the same type checking as det and the object type may persist. trace simply sums the diagonal by calling the Fraction object's __add__ method, so a Fraction object can be kept as the return type.

If you want to calculate the determinant as a rational number, you could investigate SymPy. Matrix operations such as calculating determinants are documented here.

Alex Riley
  • 169,130
  • 45
  • 262
  • 238
  • Nice answer. I found the same thing. Is it possible to correctly cast a numpy array to an arbitrary type like `Fraction`? When I tried `m.astype(Fraction)`, the array's dtype changes to `O`, but it doesn't appear to actually change the data. – Vorticity Mar 21 '15 at 12:45
  • Thanks - I don't think it's possible to have a `dtype` of anything other than `object` for Python constructs like `Fraction`. The `O` type just means that the array's values can be references to arbitrary objects: it allows greater flexibility for some methods (e.g. `trace`, `sum`) but unexpected results for others (e.g. `det`). – Alex Riley Mar 21 '15 at 13:02
  • The suggestion to use SymPy is a good one. E.g. `M = sympy.Matrix([[3, 5], [4, 3]]); M.inv().det()` gives `-1/11`. – Warren Weckesser Mar 21 '15 at 13:08
  • ...and you can use `Fraction` objects. With `m = sympy.Matrix([[Fraction(1, 6), Fraction(8, 7)], [Fraction(1, 2), Fraction(3, 2)]])`, `m.det()` give `-9/28`. – Warren Weckesser Mar 21 '15 at 13:15
  • Thanks for your clear explanations! And sympy seems great for this, I will use it. – Tom Cornebize Mar 21 '15 at 13:50
1

It looks to me like this is not an issue that will be easily solved and may be a limitation of the fact that np.linalg relies on lapack for most of its operations. Looking at the source code for numpy.linalg it appears that a routine called _commonType is called prior to calling any lapack routine. This attempts to find the appropriate type for the data contained in the input array, but if it is unable to determine the type, it assumes the type is double. The array is the cast to the resulting type prior to being passed to the lapack routine. This was likely done since it would be next to impossible to deal with every type that could be passed.

I've never worked with the Fraction package, so I can't give you a viable solution to get back to a matrix of Fraction objects. I was going to suggest calling m.astype(Fraction), but that doesn't seem to do it either.

  • List item
Vorticity
  • 4,582
  • 4
  • 32
  • 49