1

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?

Community
  • 1
  • 1
Nick-H
  • 342
  • 1
  • 3
  • 15
  • 5
    What makes you think this is possible? The determinant of a dense symbolic NxN matrix has `math.factorial(N)` terms. So, e.g., for an 18x18 matrix, that's 6,402,373,705,728,000 terms. I'd say we'll be dead before that finishes, except you'll run out of RAM long before then ;-) – Tim Peters Apr 20 '17 at 04:30
  • @TimPeters: I think that should be an answer. – Wrzlprmft Apr 20 '17 at 05:30
  • 3
    @TimPeters, I'm not sure where you are getting N factorial from, but there are far more efficient algorithms to calculate a determinant than standard co-factor expansion. For example, using LU decomposition can DRASTICALLY reduce the computational expense required in calculating a determinant of a large matrix. For a 15x15 matrix, calculating a determinant using co-factor expansion will require ~2.3 trillion multiplications. Using LU factorization will require only around 1100 multiplications to calculate the determinant. There must be something I can change to make this run faster. – Nick-H Apr 20 '17 at 06:28
  • @metalheadzone Your matrix is __symbolic__. The determinant of a numeric answer is a single number. The determinant of a generic symbolic matrix of size N by N is an expression with N! terms, as Tim Peters pointed out. The __answer__ is not even close to fitting into RAM, let alone whatever computation it takes to obtain it. –  Apr 20 '17 at 13:49
  • @Gerry, as far as I understand, the code simplifies the algebra as it goes. for a 318x318 matrix i will end up with a 318th order polynomial so there will be 319 terms, not N!. I know this computation is possible because Mathematica can do it for a symbolic 318 x 318 matrix in relatively short time. I'm just not sure which algorithm they implement or what they may be doing to speed up the process. – Nick-H Apr 20 '17 at 14:11
  • 2
    So it seems your matrix has special structure: it's not like [[a,b],[c,d]] with symbols a,b,c,d but more like [[x+1, x], [x-1, x-3]], with one symbol in the matrix. This is important because then the answer is indeed much smaller, and one can hope to optimize the computations. Please edit your question to describe the structure of the symbolic matrix. –  Apr 20 '17 at 14:41
  • @Gerry I'm sorry, I was not clear enough on that. Yes the matrix is symbolic in only one symbol. There is one symbol, omega, and a bunch of entries are purely numeric. – Nick-H Apr 20 '17 at 17:38
  • 1
    Okay, but if those are totally generic functions of omega, one cannot expect any simplifications in the determinant. Are they linear functions of omega, or polynomials of some other degree? Are the coefficients integers or floats? –  Apr 20 '17 at 19:59
  • @Gerry, all of the omega terms are second order polynomials. – Nick-H Apr 20 '17 at 22:48

2 Answers2

2

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):

  • 74.6 seconds for direct computation,
  • 5.4 seconds for the interpolation.

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]
1

As a follow up to anyone else looking at this thread: Since trying to solve this problem a few years ago, I've learned a lot more about numerical methods and general computation and realized just how infeasible taking a symbolic determinant of a matrix that large is. I ended up solving this problem numerically by converting it to an eigenvalue problem. Moral of the story... there's usually multiple ways of solving a problem and some may be more feasible than others.

Nick-H
  • 342
  • 1
  • 3
  • 15