How to speed up the code?

edit

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.

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