# How to generate all involutive permutations?

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?

• @Rolf, if I read that correctly, that skims through all possible permutationa, something that the OP already obeserved to be infeasible. (Your advantage however has the advantage that it needs less memory.) – Henrik Schumacher May 11 '18 at 18:07
• The involutive permutations are the permutations that are composed of transpositions and fixed points. So that would be one way of constructing them by hand: find permutations that fix all points (identity), find permutations that fix all but 2 points (1 transposition), etc. – anderstood May 11 '18 at 18:14

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 = DeveloperPackedArrayQ;
ToPack = DeveloperToPackedArray;

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 = CompileGetElement[idx, k, 1];
j = CompileGetElement[idx, k, 2];
x = CompileGetElement[a, i];
a[[i]] = CompileGetElement[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

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 *)


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]

• Nice, thank you! I guess the Combinatorica is not brute-forcing the permutations, right? Is the code available by any chance? – AccidentalFourierTransform May 11 '18 at 19:17
• The code can also be read directly in Combinatorica.m, and is described in the Combinatorica book. – Szabolcs May 11 '18 at 19:53
• I've been looking at Sage recently, and what's really nice is that not only is the code open (which is not big deal—every open source package has it), but it is usually well documented. The math is explained, the algorithm is described in plain English, there are references, etc. The code is not only open, but it is actually useful for learning what the function does. That is a rarity. – Szabolcs May 11 '18 at 19:55

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},
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}}