ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Sun, 01 Jan 2023 02:35:32 +0100Calculating large symbolic determinants is too slowhttps://ask.sagemath.org/question/65646/calculating-large-symbolic-determinants-is-too-slow/# 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. According to [Complexity of calculating the determinant - Mathematics Stack Exchange](https://math.stackexchange.com/questions/3561347/complexity-of-calculating-the-determinant), the time complexities for calculating determinants can lower to O(n^2.38), however, Sage seems to run in O(n!) from my measurement. 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()IvanaGyroSun, 01 Jan 2023 02:35:32 +0100https://ask.sagemath.org/question/65646/How to take the index of a variable in sagemath?https://ask.sagemath.org/question/64755/how-to-take-the-index-of-a-variable-in-sagemath/I used the method in the [post](https://ask.sagemath.org/question/7715/symbolic-linear-algebra/) to generate a matrix with entries which have indices.
def symbolic_matrix(root,m,n):
mlist=[]
for i in range(m):
for j in range(n):
mlist.append(root+'_'+str(i)+'_'+str(j))
return Matrix(SR,m,n,mlist)
k,n=4,8
r1=symbolic_matrix('y',k,n)
r1
Now I would like to get the index of the entries of the matrix r1. For example, the index of y_1_2 is [1,2]. How to do this in sagemath? I tried something like y.index(y_1_2) but it doesn't work. Thank you very much.lijr07Thu, 03 Nov 2022 18:39:29 +0100https://ask.sagemath.org/question/64755/LU decomposition for symbolic matrices with many different variables.https://ask.sagemath.org/question/62885/lu-decomposition-for-symbolic-matrices-with-many-different-variables/I asked a question about symbolic matrices and LU decompositions in https://ask.sagemath.org/question/62878/how-to-define-a-matrix-with-variables-in-sagemath/
I have another question about LU decomposition. I would like to define a matrix $x_i(t)$ as a matrix with $1$'s on the diagonal and with $t$ at $(i,i+1)$ position, and all other entries are $0$. Let $g = x_{i_1}(t_1) \cdots x_{i_m}(t_m) T$, where $i_1, \ldots, i_m$ are some given integers, $T$ is a diagonal matrix with entries $\lambda_1, \ldots, \lambda_n$. I would like to define this matrix in Sage and apply LU decomposition to it. Now $t_1, \ldots, t_m, \lambda_1, \ldots, \lambda_n$ are symbols. How to do this in Sage? Thank you very much.lijr07Sun, 19 Jun 2022 09:35:28 +0200https://ask.sagemath.org/question/62885/Matrix constructor runs out of memoryhttps://ask.sagemath.org/question/60673/matrix-constructor-runs-out-of-memory/I am trying to construct a matrix with entries in $\mathbb{Q}(q)$, from a list of Sage vectors. I do this using
`A = matrix(QQ['q'].fraction_field(), relations)
`
where the vectors in the list `relations` are Sage vectors over `QQ['q'].fraction_field()`.
The list `relations` is generated by a subroutine involving some randomness, and in many cases the above works fine but sometimes my script uses too much memory and is killed by the operating system. Tracing the memory allocation it is clear that the problem occurs at the given line where I construct the matrix from the list of vectors.
The sorts of matrices for which the script fails are tend to have about 100 columns and 50-100 rows. They are sparse (each row has <= 4 nonzero entries) and the nonzero entries are of the form `-2/q^x + 2` where `x` is in the range `10000 - 100000`.
It seems clear that when some combination of the order of magnitude of the exponents `x` and the matrix size is reached, then the matrix constructor uses too much memory. I wondered if there is a solution anyone can suggest, or even just explain a bit more about why this constructor seems to max out on memory?
(P.S. I plan to re-implement things using sparse matrices, and maybe this will help, though I'd still like to understand better why the dense version fails.)PatrickKinnearThu, 13 Jan 2022 19:01:24 +0100https://ask.sagemath.org/question/60673/Are there fast(er) ways to compute inverses of Hermitian matrices?https://ask.sagemath.org/question/51638/are-there-faster-ways-to-compute-inverses-of-hermitian-matrices/I'm dealing with some Hermitian matrices valued in symbolic expressions. I need to be able to compute inverses of these matrices, but they seem to get big pretty fast. That is it would nice to be able to do this in a reasonable amount of time with 10x10 matrices at least, and hopefully larger than that.
I'm currently running a cell where the inverses of 2 10x10 matrices are being computed, and it's taken well over 10 minutes.
Are there ways to exploit the fact that these matrices are Hermitian to compute the inverses faster, and/or are there better ways to compute inverses of symbolic matrices than the standard .inverse()??
EDIT: I originally avoided an example, because their construction is somewhat convoluted as you'll see below.
First we have two variables which for which I want to run the larger program for various choices thereof.
p = 1, q= 5
X = {(i): var("X_{}".format(i), latex_name="X_{}") for i in range(2*p*q)}
e = {(i,j,k): var("e_{}{}{}".format(i,j,k), latex_name="e_{{{}{}{}}}".format(i,j,k)) for i in range(2) for j in range(q) for k in range((p+q))}
The following creates a list L such that L[i] is the collection of lists of length i whose elements are integers between 0 and q-1 in strictly increasing order. This is used throughout the larger program.
L = [[[]]]
LL = []
for i in range(q):
LL.append([i])
L.append(LL)
if q>=2:
for k in range(2,q+1):
LLL = []
for i in range(binomial(q, k-1)):
for j in range(q):
if L[k-1][i][len(L[k-1][i])-1]<j:
mm = []
for t in L[k-1][i]:
mm.append(t)
mm.append(j)
LLL.append(mm)
L.append(LLL)
Here, a matrix is defined which we be used to construct the matrices I'm interested in.
HE = matrix(SR, q, q)
for i in range(q):
for j in range(q):
prod = 0
for k in range(p):
prod -= (e[(0, i, k)]*e[(1, j, k)])
for l in range(p+q):
if l>=p:
prod += (e[(0, i, l)]*e[(0, j, l)])
HE[i,j] = prod
h = {(i): var("h_{}".format(i)) for i in range(q+1)}
for i in range(q+1):
h[i] = matrix(SR, binomial(q, i), binomial(q, i))
h[0][0,0] = 1
for i in range(1, q+1):
for z in L[i]:
for w in L[i]:
h[i][L[i].index(z), L[i].index(w)] = HE[z, w].det()
Finally, we come to the inverse matrices I would like to compute. It is absolutely necessary that I have these inverse matrices.
hinv = {(i): var("hinv_{}".format(i)) for i in range(q+1)}
for i in range(q+1):
hinv[i] = h[i].inverse()sum8tionSat, 30 May 2020 07:51:15 +0200https://ask.sagemath.org/question/51638/<symbolic matrix>.simplify_rational() accepts no keywords?https://ask.sagemath.org/question/49934/symbolic-matrixsimplify_rational-accepts-no-keywords/I have a symbolic matrix whose entries I would like to simplify for display. Minimal working example (in the `sage` REPL):
```
a,b=var('a,b'); m=matrix(SR,2,[[(a*b-a)/a,0],[0,1]]).simplify_rational()
```
works fine and upon calling `m` produces
```[b - 1 0]
[ 0 1]
```
as expected. On the other hand, running
```
a,b=var('a,b'); m=matrix(SR,2,[[(a*b-a)/a,0],[0,1]]).simplify_rational(algorithm='simple')
```
as per the documentation produces
TypeError Traceback (most recent call last)
<ipython-input-30-3edf7642200d> in <module>()
----> 1 a,b=var('a,b'); m=matrix(SR,Integer(2),[[(a*b-a)/a,Integer(0)],[Integer(0),Integer(1)]]).simplify_rational(algorithm='simple')
TypeError: simplify_rational() takes no keyword argument
If I apply the method to a *single* expression it works fine:
```
((a*b-a)/a).simplify_rational(algorithm='simple')
```
gives back `b-1`.
grobberSat, 15 Feb 2020 18:03:18 +0100https://ask.sagemath.org/question/49934/Optimal method to invert symbolic matrixhttps://ask.sagemath.org/question/47401/optimal-method-to-invert-symbolic-matrix/ If I read the source code correctly, a matrix defined over SR is inverted by the Moore-Penrose pseudo-inversion. For a real, symmetric, symbolic, full-rank matrix, is this algorithm optimal? Or is there another algorithm that would be a) less memory intensive, and/or b) less time-consuming. [In my particular case, the attempted inversion of a 4x4 matrix runs out of memory (8 GB) after approximately 15 hours.]Richard_LTue, 06 Aug 2019 01:04:16 +0200https://ask.sagemath.org/question/47401/Formal determinant of symbolic matrixhttps://ask.sagemath.org/question/46624/formal-determinant-of-symbolic-matrix/ I have some sparse symbolic matrices, and want to compute their formal determinant (without cancellation of terms). In other words, if I have the matrix
x,y = var('x,y')
M = Matrix(SR, [[x,y],[x,y]])
I would like the result of
M.determinant()
to be x*y - x*y, rather than just 0. The variables in each monomial are allowed to commute with each other, but on the other hand I would like all monomials containing a 0 to vanish (i.e if in the example above M = Matrix(SR, [[x,0],[x,y]])), then the determinant should be just x*y, rather than x*y - 0*x).
Is there a way to achieve this (without using the expansion of the determinant as permutations, since the dimension of the matrices gets quite big!)? Thanks in advance!
danieleCWed, 22 May 2019 16:32:20 +0200https://ask.sagemath.org/question/46624/