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