# How to speed up the code?

Given a matrix M, to find all 1,-1 matrices which satisfy the condition given in the code. For n=3, I am getting the outputs easily. For n=6, it ran for around 12 hours and got "flint exception error". I understand that len(T) is huge so there might be some memory issues. Is there anyway to get over it?

 n=3
T=Tuples((1,-1),n^2)
A=[matrix(n,n,t) for t in T]
for a in A:
if 3*a^2==M:
show(a)

edit retag close merge delete

Can you give an example of matrix M resulting in an error?

( 2022-10-10 20:27:40 +0200 )edit
1

As for speeding up the code it may be worth to compute the square root of matrix M/3and check if it composed of $\pm1$.

( 2022-10-10 20:34:06 +0200 )edit

An example of M is:

  M = -n *ones_matrix(n) + (4*n)*(identity_matrix(n))

( 2022-10-11 04:47:15 +0200 )edit

Sort by » oldest newest most voted

A simple modification would be to not create the list of matrices:

 n = M.nrows()
for t in Tuples((1,-1), n^2):
a = matrix(n, n, t)
if 3*a^2 == M:
show(a)


With n=5, this at least enters the loop immediately, whereas the original code spends a while creating the list of matrices. I gave up waiting for the original code after about 5 minutes. (You can add some print statements to see its progress.)

By the way, you can see that Tuples is fast, but creating a list from its entries can be slow:

sage: %time Tuples((1, -1), 16)
CPU times: user 24 µs, sys: 8 µs, total: 32 µs
Wall time: 46 µs
Tuples of (1, -1) of length 16

sage: %time L = [t for t in Tuples((1, -1), 16)] # slow
CPU times: user 453 ms, sys: 23.6 ms, total: 477 ms
Wall time: 479 ms

sage: %time L = [matrix(4, 4, t) for t in Tuples((1, -1), 16)] # slower
CPU times: user 988 ms, sys: 15.6 ms, total: 1 s
Wall time: 1.01 s


Going through all of these matrices is going to take a long time, in any case:

sage: len(Tuples((1, -1), 25))
33554432
sage: len(Tuples((1, -1), 36))
68719476736


Edit: here is a modified version of your code:

def f(M):
results = []
n = M.nrows()
j = 0
for t in Tuples((1,-1), n^2):
a = matrix(n, n, t)
if j % 10000 == 0:
print(j)
j += 1
if 3*a^2 == M:
results.append(a)
# if you want to stop after finding one solution, you could add "break" here to stop the loop.
return(results)


You can then call f(M) for any given matrix M. It prints something after every 10,000 matrices, and it saves the results and returns them at the end.

more

Ok, thanks for the information

( 2022-10-11 04:51:42 +0200 )edit
1

In fact I would guess that your original error occurred in the creation of the list A: I don't know if Python will successfully create a list of 2**36 6-by-6 matrices.

( 2022-10-11 20:09:37 +0200 )edit

I started my computer on the 6x6 case with the matrix you mentioned above (-n *ones_matrix(n) + (4*n)*(identity_matrix(n)), and at the rate it's going, it looks like it would take a few months.

( 2022-10-11 23:34:23 +0200 )edit

Your idea behind "modified version" is really very good . Thanks . I also ran it on my laptop a few hours back. Yes, it looks it is going to take a lot of time.

( 2022-10-12 06:12:18 +0200 )edit

One can employ matrix square root to significantly speed up things. The following implementation illustrates this idea and supports only matrices with distinct eigenvalues:

def sqrt_pm1(M):
if not M.det().is_square():    # quick test for insolubility
return []

n = M.nrows()

d, L = M.eigenmatrix_right()
d = d.diagonal()
assert len(set(d)) == n          # only this case is supported
L = L.change_ring(QQbar)

pm1 = {QQbar(-1),QQbar(1)}
return [R for signs in Tuples((-1,1),n) if set((R := L * diagonal_matrix([sqrt(e)*s for e,s in zip(d,signs)]) * L.inverse()).list()).issubset(pm1)]


For example,

sqrt_pm1(Matrix([[0,-6],[6,0]]) / 3)


gives the following list of two $\pm 1$-matrices that are square roots of Matrix([[0,-6],[6,0]]) / 3:

[
[-1  1]  [ 1 -1]
[-1 -1], [ 1  1]
]

more

Okay...Thanks a lot for sparing your time. For M = -n ones_matrix(n) + (4n)*(identity_matrix(n)), mentioned above , eigen values are not distinct. Can anything be done in that case?

( 2022-10-12 06:14:57 +0200 )edit
1

At very least you can quickly filter out unsolvable M by checking (M/3).det().is_square()

( 2022-10-12 12:08:53 +0200 )edit

Okay. Thanks

( 2022-10-12 19:00:01 +0200 )edit