This is an NP-complete problem - see https://mathoverflow.net/q/155468 for a reference. So, it's unlikely for a truly quick (ie. polynomial-time) way to exist. But there is also a reference to Bhat (1981) paper, which proposes some optimized algorithm for this problem.

**ADDED.**
Here is an ILP solution to this problem that for a given matrix `M`

tries to find a permutation matrix `Q`

(if one exists) such that `M*Q`

is symmetric:

```
# matrix size
n = 10
# generate a random symmetric n X n 01-matrix
M = Matrix(ZZ,n,n)
for i in range(n):
for j in range(i+1):
M[i,j] = M[j,i] = randint(0,1)
# randomly permute columns of M
p = SymmetricGroup(n).random_element()
M = M[:,list(p(i)-1 for i in (1..n))]
# problem instance is created, let's now try to solve it
# define an ILP problem
milp = MixedIntegerLinearProgram()
milp.set_objective( None )
# permutation matrix to be found
Q = milp.new_variable(binary=True)
for i in range(n):
milp.add_constraint( sum( Q[(i,j)] for j in range(n) ) == 1 )
milp.add_constraint( sum( Q[(j,i)] for j in range(n) ) == 1 )
# constraints imposing that M*Q is symmetric
for i in range(n):
for j in range(i):
milp.add_constraint( sum( M[i,k] * Q[(k,j)] for k in range(n) ) == sum( M[j,k] * Q[(k,i)] for k in range(n) ) )
milp.solve()
R = Matrix( milp.get_values(Q) )
# print the permuted matrix, which should be symmetric
print(M*R)
```