I have encountered what looks like a bug relating to tensor symmetry computations (or maybe I am missing something?)

Take the following example:

Table[E^(-(((i - j)/2 - p + q)^2/2)), {i, 5}, {j, 5}, {p, 3}, {q, 3}]

This array is symmetric if you swap slots 1 with 2 and 3 with 4. TensorSymmetry should give {{Cycles[{{1, 2}, {3, 4}}], 1}}. However, it doesn't detect any symmetries at all. Checking

\[Element] Arrays[{5, 5, 3, 3}, Reals, {{Cycles[{{1, 2}, {3, 4}}], 1}}]

Also returns false. Passing to SymmetrizedArray also returns an array with no symmetries. Note that if you explicitly specify the symmetry it will work.

SymmetrizedArray[{i_, j_, p_, q_} :> E^(-(((i - j)/2 - p + q)^2/2)), {5, 5, 3, 3}, {{{2, 1, 4, 3}, 1}}]

The array you get will be symmetric and works with \[Element] and TensorSymmetry. In fact, it will pass an equality check with the non-symmetrized version, proving the array is indeed symmetric.

SymmetrizedArray[{i_, j_, p_, q_} :> E^(-(((i - j)/2 - p + q)^2/2)), 
 {5, 5, 3, 3}, {{{2, 1, 4, 3}, 1}}] == 
 Table[E^(-(((i - j)/2 - p + q)^2/2)), {i, 5}, {j, 5}, {p, 3}, {q, 3}]

Gives True, and

Symmetrize[
  Table[E^(-(((i - j)/2 - p + q)^2/2)), {i, 5}, {j, 5}, {p, 3}, {q, 3}], {{{2, 1, 4, 3}, 1}}] == 
 Table[E^(-(((i - j)/2 - p + q)^2/2)), {i, 5}, {j, 5}, {p, 3}, {q, 3}]

also gives True.

I suspect there is a bug in the symmetry finding or checking code. Here I am using exact numbers, so there shouldn't be any tolerance issues. But if you use approximate numbers you get the same issue. This is on Wolfram Desktop 13.2.1.0 on MacOS ARM.