Search code examples
pythonnumpymatrixrational-number

Why doesn't numpy determinant return a Fraction when given a Fraction matrix?


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()).


Solution

  • 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.