How to generate all involutive permutations?

Question

Take a finite set $S$ (i.e., a list). An involutive permutation is one that squares to the identity. How can we generate all such permutations efficiently, that is, without generating all permutations first, and eliminating the non-involutive ones?

For example, the code below generates all involutive permutations on the set {1,2,3,4}, by first generating all permutations:

With[{perm = Permutations[Range[4]]}, perm[[Flatten@Position[Map[Apply[ReplaceAll], Transpose[{perm, Thread /@ Thread[ConstantArray[Range[4], 4!] -> perm]}]], Range[4]]]]]

The code returns 10 out of the 24 permutations of {1,2,3,4}. The reduction for larger lists is more extreme (e.g., 76 vs 720 for a list of size 6). I need to consider a rather large list, so generating all permutations first is unfeasible. How can we make the calculation as efficient as possible?

Answer 1

I wondered how well Involutions has been implemented, so I tried to reimplement it myself. The following implementation can be up to 15 times faster than Involutions.

Needs["Combinatorica`"];
PackedQ = Developer`PackedArrayQ;
ToPack = Developer`ToPackedArray;

myInvolutions[list_List] := Block[{data, A, n, g},
  n = Length[list];
  A = ToPack[
    UpperTriangularize[{Range[3, n]}[[ConstantArray[1, n - 1]]]]];
  A[[2 ;;]] += ToPack[LowerTriangularize[{Range[2, n - 1]}[[ConstantArray[1, n - 2]]]]];
  g[2, {{}}] = ToPack[{{{1, 2}}}];
  g[n_, data_] := With[{m = Length[data]},
    Join[
     Transpose[{{
        ConstantArray[1, {(n - 1) m}],
        Flatten[Table[ConstantArray[i, {m}], {i, 2, n}]]
        }},
      {2, 3, 1}
      ],
     ArrayReshape[
      A[[1 ;; (n - 1), Flatten[data]]], {(n - 1) m, 
      Quotient[n - 2, 2], 2}
      ], 
    2]
    ];

  data = {{}};
  Join @@ Join[
    {{list}},
    Table[
     data = g[2 i, data];
     getPermutationLists[
      list,
      ArrayReshape[
       Subsets[Range[n], {2 i}][[All, Flatten[data]]],
       {Binomial[n, 2 i] Length[data], Sequence @@ Rest[Dimensions[data]]}
       ]
      ],
     {i, 1, Quotient[n, 2]}]
    ]
  ]

getPermutationLists = Compile[{{ran, _Integer, 1}, {idx, _Integer, 2}},
   Block[{a = ran, i, j, k, x},
    Do[
     i = Compile`GetElement[idx, k, 1];
     j = Compile`GetElement[idx, k, 2];
     x = Compile`GetElement[a, i];
     a[[i]] = Compile`GetElement[a, j];
     a[[j]] = x,
     {k, 1, Length[idx]}
     ];
    a
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

Here is a test:

n = 14;
a = Range[n];
aa = myInvolutions[a]; // AbsoluteTiming // First
bb = Involutions[a]; // AbsoluteTiming // First
Sort[aa] == Sort[bb]

5.63692

83.2192

True

Answer 2

Hard to beat Henrik's answer, but I'd like to show how you could do it by hand, without using a package.

From the decomposition of permutations into cycles, you can see that involutive permutations are the compositions of fixed points and transpositions.

Accordingly, you can build all the fixed points i (for identity) and transpositions t with:

n = 10;
a = Range[n];
t = Subsets[a, {2}];
i = Subsets[a, {1}];

Then, the transpositions will be among

s = Subsets[t, Floor[n/2]];

Unfortunately, s also includes non-disjoint cycles (like {1, 2}, {1, 3}). I don't know how to assemble only disjoint cycles together, but this is an unefficient brute-force approach:

inv = Select[s, Length[Flatten@#] == Length[DeleteDuplicates@Flatten@#] &];

That's a bit stupid because it builds one million subsets for $n=10$, to keep only 1000 in the end. This can be improved. Note however that contrary to your fully brute-force approach, here it's only exhaustive for transpositions, so it's still much much better.

Now, the number of permutations is given by

Length@inv
(* 9496 *)

And you can also rebuild the permutations matrices from inv:

cycles = Cycles /@ inv;
invo = Permute[Range[n], #] & /@ cycles
And @@ Map[PermutationProduct[#, #] == Range[n] &, invo]
(* True *)

Answer 3

The partially obsolete Combinatorica package has a function called Involutions that returns all involutive permutations:

Needs["Combinatorica`"];
a = Range[10];
inv = Involutions[a];
And @@ Map[PermutationProduct[#, #] == a &, inv]

True

Note also the existence of NumberOfInvolutions. Guess what it does...

The code for Involutions is actually available and can be printed as follows:

Needs["Combinatorica`"];
Needs["GeneralUtilities`"];
PrintDefinitions[Involutions]

Answer 4

Another implementation based on the idea mentioned by @anderstood i.e. involutions are permutations that can be expressed as product of disjoint transpositions(or 2-cycles).

OrderlessPartition[set_List, part_List] := 
 Module[{f, part2, list}, Attributes[f] = Orderless;
  part2 = Sort@part;
  list = Unique["x", Temporary] & /@ part2; 
  ReplaceList[f @@ set, 
   With[{list1 = list, list2 = List /@ list}, 
    f @@ (MapThread[
         Pattern[#1, Repeated[_, {#2}]] &, {list1, part2}]) /; 
      OrderedQ[list2] -> List /@ list1]]]

 TwoSidedRule[a_] := Nothing;
 TwoSidedRule[a_, b_] := Sequence[a -> b, b -> a];

  Involutions[set_List] := (set /. # &) /@ 
   Apply[TwoSidedRule, Flatten[#, 1] &@(OrderlessPartition[set, #] & /@ 
       IntegerPartitions[Length@set, Length@set, {1, 2}]), {2}];

Now, for the sample example of $\{1,2,3,4\}$, we have:

Involutions[Range[4]]

{{2, 1, 4, 3}, {3, 4, 1, 2}, {4, 3, 2, 1}, {1, 2, 4, 3}, {1, 4, 3, 2}, {1, 3, 2, 4}, {4, 2, 3, 1}, {3, 2, 1, 4}, {2, 1, 3, 4}, {1, 2, 3, 4}}