Some research that I am working on requires symbolically taking the determinant of large matrices; matrices ranging from 18x18 to 318x318. The matrix entries are either numeric or polynomials of degree two in the same variable, omega.
Currently, I am trying to use the .det() method in SymPy but it is very slow; an 18x18 matrix has been running for over 45 minutes now and is still computing as I write this. I realize that determinant calculations are very intensive, but is there anything I can do to speed this up?
I've already read the post at Speeding up computation of symbolic determinant in SymPy but did not take anything away from the post as to what could actually be done to speed the process up. What can I do?
SymPy is not naive about determinants (see MatrixDeterminant class) but it appears that juggling symbolic expression throughout the computation is a slow process. When the determinant is known to be a polynomial of certain degree (because the matrix entries are), it turns out to be faster to compute its numeric values for several values of the variable, and interpolate.
My test case is a dense 15 by 15 matrix full of quadratic polynomials of variable omega, with integer coefficients. I still use SymPy's .det method for the numeric determinants, so the coefficients end up being exactly the same long integers either way.
import numpy as np
from sympy import *
import time
n = 15
omega = Symbol('omega')
A = Matrix(np.random.randint(low=0, high=20, size=(n, n)) + omega*np.random.randint(low=0, high=20, size=(n, n)) + omega**2 * np.random.randint(low=0, high=20, size=(n, n)))
start = time.time()
p1 = A.det() # direct computation
print('Time: ' + str(time.time() - start))
start = time.time()
xarr = range(-n, n+1) # 2*n+1 points to get a polynomial of degree 2*n
yarr = [A.subs(omega, x).det() for x in xarr] # numeric values
p2 = expand(interpolating_poly(len(xarr), omega, xarr, yarr)) # interpolation
print('Time: ' + str(time.time() - start))
Both p1 and p2 are the same polynomial. Running time (on a pretty slow machine, t2.nano from Amazon):
If your coefficients are floating point numbers and you don't expect exact arithmetical results when dealing with them, further speedup may be achieved by evaluating the matrix as a NumPy array, and using a NumPy method for the determinant:
Anum = lambdify(omega, A)
yarr = [np.linalg.det(Anum(x)) for x in xarr]