Calculating large symbolic determinant is too slow
I want to calculate a series of symbolic Hankel matrices. The sizes of the matrices are from 2*2 to 50*50. I found the speed of calculating the determinants is too slow. It took 5 seconds to get the result of a 7*7 matrix, 45 seconds to get the result of an 8*8 matrix, and it was always killed when calculating the 9*9 determinants. How can I make the calculation more efficient?
Here is my full code.
import time
from functools import cache
import sage.libs.ecl
from sympy import simplify
# increase memory limit
sage.libs.ecl.ecl_eval("(ext:set-limit 'ext:heap-size 0)")
time_table = {}
l = 1
energy = var('E')
def begin(name):
time_table[name] = time.time()
def end(name):
print(f'{name}: {time.time() - time_table[name]}')
del time_table[name]
@cache
def expected_distance(n, l):
if n <= -2:
return 0
if n == -1:
return -2 * energy
if n == 0:
return 1
return simplify((
- 4 * (2 * (n + 1) - 1) * expected_distance(n - 1, l) \
- n * ((n + 1) * (n - 1) - 4 * l * (l + 1)) * expected_distance(n - 2, l)
) / (8 * (n + 1) * energy))
min_rank = 2
max_rank = 50
positive_range_graph = Graphics()
xmin = -0.7
xmax = 0.5
current_ranges = [[xmin, xmax]]
offset = 1
for rank in range(min_rank, max_rank + 1):
print(f'rank:{rank}')
begin(f'rank:{rank}')
begin('generate matrix')
hankel = matrix.hankel([expected_distance(i + offset, l) for i in range(rank)],
[expected_distance(i + offset + rank - 1, l) for i in range(1, rank)], SR)
end('generate matrix')
begin('calculate det')
# pretty slow!
# It is always killed when the matrix size is 9*9
f = hankel.det()
end('calculate det')
begin('simplify')
f = f.simplify_full()
# print(f'rank:{rank} {f}' )
end('simplify')
begin('solve inequation')
positive_solutions = solve(f >= 0, energy)
end('solve inequation')
begin('plot')
positive_ranges = []
for solution in positive_solutions:
if len(solution) > 2:
print(f'more than two expression')
if len(solution) == 1:
if solution[0].operator() in [operator.le, operator.lt]:
positive_ranges.append((xmin, solution[0].rhs()))
elif solution[0].operator() in [operator.ge, operator.gt]:
positive_ranges.append((solution[0].rhs(), xmax))
elif solution[0].operator() == operator.eq:
positive_ranges.append((solution[0].rhs(), solution[0].rhs()))
else:
print(f'Cannot recognize operator: {solution[0]}')
elif len(solution) == 2:
positive_ranges.append((solution[0].rhs(), solution[1].rhs()))
hue = (rank - min_rank) / (max_rank - min_rank)
for positive_range in positive_ranges:
dots = [(x, rank) for x in positive_range]
positive_range_graph += line(dots, hue=hue)
positive_range_graph += points(dots, color='red')
end('plot')
end(f'rank:{rank}')
print()
