Wolfram Function Repository

Instant-use add-on functions for the Wolfram Language

Function Repository Resource:

 MaximizeOverPermutations

Source Notebook

Find the permutation that maximizes a function

Contributed by: Roman Schmied, University of Basel

ResourceFunction["MaximizeOverPermutations"][f,n]

maximizes f numerically with respect to permutations of {1,…,n}.

ResourceFunction["MaximizeOverPermutations"] takes a Method option that accepts either "Enumerate" or "MonteCarlo".
The number of permutations of n objects scales very quickly with n. The default "Enumerate" method should not be used with n≳12 in order to prevent memory overflow, as the time and space required for this enumeration both scale with n!. A warning is issued, but the user needs to decide whether or not to call for full enumeration of larger permutation spaces.
With the default option setting Method→"Enumerate", the maximum is determined by applying the function f to all permutations of the numbers 1… n. The result is guaranteed to be optimal.
With the option setting Method→"MonteCarlo", the maximum is determined by a Monte Carlo Metropolis–Hastings stochastic search. There is no guarantee of the result’s optimality. This method makes it possible to search very large permutation spaces efficiently.
With the suboption settings Method→{"MonteCarlo",options…}, the Monte Carlo method can be fine-tuned as follows:
The suboption "Iterations"→m specifies the number of iterations that the Monte Carlo algorithm performs in the stochastic search. The algorithm's runtime scales linearly with m. The default value is "Iterations"→104.
The suboption "AnnealingParameter"→β specifies the annealing parameter of the Metropolis–Hastings acceptance/rejection step. Smaller values of β give a wider search; larger values of β tend to restrict the search more toward a monotonically increasing strategy. The default value is "AnnealingParameter"→1.
With "AnnealingParameter"→{βstart,βend}, the annealing parameter progresses exponentially from βstart at the simulation’s beginning to βend at the end. An increasing annealing parameter is useful to narrow the search progressively during the course of the simulation.
The suboption "StartingPermutation"→σ specifies the initial point of the stochastic search in permutation space and must be given as a list σ={s1,s2,…,sn} that is a permutation of the numbers 1…n. The default value is "StartingPermutation"→Automatic, which uses the trivial permutation {1,2,…,n} as the starting point.
The Monte Carlo Metropolis–Hastings stochastic search relies on the annealing parameter β, which can be kept uniform or varied exponentially during the simulation. A proposed random step from a permutation σi to a new permutation σi+1 is accepted with probability Max[1,Exp[β(f[σi+1]-f[σi])]]. The annealing parameter β must therefore be chosen in relation to the typically occurring differences in the function f to be maximized. The value of this parameter is crucial for the efficiency of the simulation and must be determined by experimentation.

Examples

Basic Examples (3) 

Maximize the second value:

In[1]:=
ResourceFunction["MaximizeOverPermutations"][#[[2]] &, 3]
Out[1]=

Given a function that takes a permutation of the numbers {1,2,3}, determine the maximum function value:

In[2]:=
f[{1, 2, 3}] = 2; f[{1, 3, 2}] = 9; f[{2, 1, 3}] = 1; f[{2, 3, 1}] = -3; f[{3, 1, 2}] = 9; f[{3, 2, 1}] = 5;

The result shows that there are two permutations, {1,3,2} and {3,1,2}, that both yield the maximal value f=9:

In[3]:=
ResourceFunction["MaximizeOverPermutations"][f, 3]
Out[3]=

Construct a function that takes a permutation of the numbers {1,2,…,100} and whose global maximum is known:

In[4]:=
m = 100; x = RandomReal[{0, 1}, m]; g = # . x &;

The global maximum of g is given by the permutation that is in the same numerical order as the list x:

In[5]:=
g[Ordering[Ordering[x]]]
Out[5]=

Find the maximum of g by a Monte Carlo search though the space of m!≈9×10157 possible permutations:

In[6]:=
ResourceFunction["MaximizeOverPermutations"][g, m, Method -> {"MonteCarlo", "Iterations" -> 10^5, "AnnealingParameter" -> 10}]
Out[6]=

This result is not quite optimal, but close. It is much larger than the typical values of g found by applying random permutations directly:

In[7]:=
Histogram[Table[g[RandomSample[Range[m]]], 10^4]]
Out[7]=

Minimizing a function f can be achieved by maximizing -f:

In[8]:=
f[{1, 2, 3}] = 2; f[{1, 3, 2}] = 9; f[{2, 1, 3}] = 1; f[{2, 3, 1}] = -3; f[{3, 1, 2}] = 9; f[{3, 2, 1}] = 5;
In[9]:=
{#[[1]], -#[[2]]} &@ ResourceFunction["MaximizeOverPermutations"][-f[#] &, 3]
Out[9]=

Find the extremum of the NUG25 standard Quadratic Assignment Problem:

In[10]:=
S = 25; M1 = \!\(\* TagBox[ RowBox[{"(", "", GridBox[{ {"0", "1", "2", "3", "4", "1", "2", "3", "4", "5", "2", "3", "4", "5", "6", "3", "4", "5", "6", "7", "4", "5", "6", "7", "8"}, {"1", "0", "1", "2", "3", "2", "1", "2", "3", "4", "3", "2", "3", "4", "5", "4", "3", "4", "5", "6", "5", "4", "5", "6", "7"}, {"2", "1", "0", "1", "2", "3", "2", "1", "2", "3", "4", "3", "2", "3", "4", "5", "4", "3", "4", "5", "6", "5", "4", "5", "6"}, {"3", "2", "1", "0", "1", "4", "3", "2", "1", "2", "5", "4", "3", "2", "3", "6", "5", "4", "3", "4", "7", "6", "5", "4", "5"}, {"4", "3", "2", "1", "0", "5", "4", "3", "2", "1", "6", "5", "4", "3", "2", "7", "6", "5", "4", "3", "8", "7", "6", "5", "4"}, {"1", "2", "3", "4", "5", "0", "1", "2", "3", "4", "1", "2", "3", "4", "5", "2", "3", "4", "5", "6", "3", "4", "5", "6", "7"}, {"2", "1", "2", "3", "4", "1", "0", "1", "2", "3", "2", "1", "2", "3", "4", "3", "2", "3", "4", "5", "4", "3", "4", "5", "6"}, {"3", "2", "1", "2", "3", "2", "1", "0", "1", "2", "3", "2", "1", "2", "3", "4", "3", "2", "3", "4", "5", "4", "3", "4", "5"}, {"4", "3", "2", "1", "2", "3", "2", "1", "0", "1", "4", "3", "2", "1", "2", "5", "4", "3", "2", "3", "6", "5", "4", "3", "4"}, {"5", "4", "3", "2", "1", "4", "3", "2", "1", "0", "5", "4", "3", "2", "1", "6", "5", "4", "3", "2", "7", "6", "5", "4", "3"}, {"2", "3", "4", "5", "6", "1", "2", "3", "4", "5", "0", "1", "2", "3", "4", "1", "2", "3", "4", "5", "2", "3", "4", "5", "6"}, {"3", "2", "3", "4", "5", "2", "1", "2", "3", "4", "1", "0", "1", "2", "3", "2", "1", "2", "3", "4", "3", "2", "3", "4", "5"}, {"4", "3", "2", "3", "4", "3", "2", "1", "2", "3", "2", "1", "0", "1", "2", "3", "2", "1", "2", "3", "4", "3", "2", "3", "4"}, {"5", "4", "3", "2", "3", "4", "3", "2", "1", "2", "3", "2", "1", "0", "1", "4", "3", "2", "1", "2", "5", "4", "3", "2", "3"}, {"6", "5", "4", "3", "2", "5", "4", "3", "2", "1", "4", "3", "2", "1", "0", "5", "4", "3", "2", "1", "6", "5", "4", "3", "2"}, {"3", "4", "5", "6", "7", "2", "3", "4", "5", "6", "1", "2", "3", "4", "5", "0", "1", "2", "3", "4", "1", "2", "3", "4", "5"}, {"4", "3", "4", "5", "6", "3", "2", "3", "4", "5", "2", "1", "2", "3", "4", "1", "0", "1", "2", "3", "2", "1", "2", "3", "4"}, {"5", "4", "3", "4", "5", "4", "3", "2", "3", "4", "3", "2", "1", "2", "3", "2", "1", "0", "1", "2", "3", "2", "1", "2", "3"}, {"6", "5", "4", "3", "4", "5", "4", "3", "2", "3", "4", "3", "2", "1", "2", "3", "2", "1", "0", "1", "4", "3", "2", "1", "2"}, {"7", "6", "5", "4", "3", "6", "5", "4", "3", "2", "5", "4", "3", "2", "1", "4", "3", "2", "1", "0", "5", "4", "3", "2", "1"}, {"4", "5", "6", "7", "8", "3", "4", "5", "6", "7", "2", "3", "4", "5", "6", "1", "2", "3", "4", "5", "0", "1", "2", "3", "4"}, {"5", "4", "5", "6", "7", "4", "3", "4", "5", "6", "3", "2", "3", "4", "5", "2", "1", "2", "3", "4", "1", "0", "1", "2", "3"}, {"6", "5", "4", "5", "6", "5", "4", "3", "4", "5", "4", "3", "2", "3", "4", "3", "2", "1", "2", "3", "2", "1", "0", "1", "2"}, {"7", "6", "5", "4", "5", "6", "5", "4", "3", "4", "5", "4", "3", "2", "3", "4", "3", "2", "1", "2", "3", "2", "1", "0", "1"}, {"8", "7", "6", "5", "4", "7", "6", "5", "4", "3", "6", "5", "4", "3", "2", "5", "4", "3", "2", "1", "4", "3", "2", "1", "0"} }, GridBoxAlignment->{"Columns" -> {{Center}}, "Rows" -> {{Baseline}}}, GridBoxSpacings->{"Columns" -> { Offset[0.27999999999999997`], { Offset[0.7]}, Offset[0.27999999999999997`]}, "Rows" -> { Offset[0.2], { Offset[0.4]}, Offset[0.2]}}], "", ")"}], Function[BoxForm`e$, MatrixForm[BoxForm`e$]]]\); M2 = \!\(\* TagBox[ RowBox[{"(", "", GridBox[{ {"0", "3", "2", "0", "0", "10", "5", "0", "5", "2", "0", "0", "2", "0", "5", "3", "0", "1", "10", "0", "2", "1", "1", "1", "0"}, {"3", "0", "4", "0", "10", "0", "0", "2", "2", "1", "5", "0", "0", "0", "0", "0", "1", "6", "1", "0", "2", "2", "5", "1", "10"}, {"2", "4", "0", "3", "4", "5", "5", "5", "1", "4", "0", "4", "0", "4", "0", "3", "2", "5", "5", "2", "0", "0", "3", "1", "0"}, {"0", "0", "3", "0", "0", "0", "2", "2", "0", "6", "2", "5", "2", "5", "1", "1", "1", "2", "2", "4", "2", "0", "2", "2", "5"}, {"0", "10", "4", "0", "0", "2", "0", "0", "0", "0", "0", "0", "0", "0", "2", "0", "0", "2", "0", "5", "0", "2", "1", "0", "2"}, {"10", "0", "5", "0", "2", "0", "10", "10", "5", "10", "6", "0", "0", "10", "2", "10", "1", "5", "5", "2", "5", "0", "2", "0", "1"}, {"5", "0", "5", "2", "0", "10", "0", "1", "3", "5", "0", "0", "2", "4", "5", "10", "6", "0", "5", "5", "5", "0", "5", "5", "0"}, {"0", "2", "5", "2", "0", "10", "1", "0", "10", "2", "5", "2", "0", "3", "0", "0", "0", "4", "0", "5", "0", "5", "2", "2", "5"}, {"5", "2", "1", "0", "0", "5", "3", "10", "0", "5", "6", "0", "1", "5", "5", "5", "2", "3", "5", "0", "2", "10", "10", "1", "5"}, {"2", "1", "4", "6", "0", "10", "5", "2", "5", "0", "0", "1", "2", "1", "0", "0", "0", "0", "6", "6", "4", "5", "3", "2", "2"}, {"0", "5", "0", "2", "0", "6", "0", "5", "6", "0", "0", "2", "0", "4", "2", "1", "0", "6", "2", "1", "5", "0", "0", "1", "5"}, {"0", "0", "4", "5", "0", "0", "0", "2", "0", "1", "2", "0", "2", "1", "0", "3", "10", "0", "0", "4", "0", "0", "4", "2", "5"}, {"2", "0", "0", "2", "0", "0", "2", "0", "1", "2", "0", "2", "0", "4", "5", "0", "1", "0", "5", "0", "0", "0", "5", "1", "1"}, {"0", "0", "4", "5", "0", "10", "4", "3", "5", "1", "4", "1", "4", "0", "0", "0", "2", "2", "0", "2", "5", "0", "5", "2", "5"}, {"5", "0", "0", "1", "2", "2", "5", "0", "5", "0", "2", "0", "5", "0", "0", "2", "0", "0", "0", "6", "3", "5", "0", "0", "5"}, {"3", "0", "3", "1", "0", "10", "10", "0", "5", "0", "1", "3", "0", "0", "2", "0", "0", "5", "5", "1", "5", "2", "1", "2", "10"}, {"0", "1", "2", "1", "0", "1", "6", "0", "2", "0", "0", "10", "1", "2", "0", "0", "0", "5", "2", "1", "1", "5", "6", "5", "5"}, {"1", "6", "5", "2", "2", "5", "0", "4", "3", "0", "6", "0", "0", "2", "0", "5", "5", "0", "4", "0", "0", "0", "0", "5", "0"}, {"10", "1", "5", "2", "0", "5", "5", "0", "5", "6", "2", "0", "5", "0", "0", "5", "2", "4", "0", "5", "4", "4", "5", "0", "2"}, {"0", "0", "2", "4", "5", "2", "5", "5", "0", "6", "1", "4", "0", "2", "6", "1", "1", "0", "5", "0", "4", "4", "1", "0", "2"}, {"2", "2", "0", "2", "0", "5", "5", "0", "2", "4", "5", "0", "0", "5", "3", "5", "1", "0", "4", "4", "0", "1", "0", "10", "1"}, {"1", "2", "0", "0", "2", "0", "0", "5", "10", "5", "0", "0", "0", "0", "5", "2", "5", "0", "4", "4", "1", "0", "0", "0", "0"}, {"1", "5", "3", "2", "1", "2", "5", "2", "10", "3", "0", "4", "5", "5", "0", "1", "6", "0", "5", "1", "0", "0", "0", "0", "0"}, {"1", "1", "1", "2", "0", "0", "5", "2", "1", "2", "1", "2", "1", "2", "0", "2", "5", "5", "0", "0", "10", "0", "0", "0", "2"}, {"0", "10", "0", "5", "2", "1", "0", "5", "5", "2", "5", "5", "1", "5", "5", "10", "5", "0", "2", "2", "1", "0", "0", "2", "0"} }, GridBoxAlignment->{"Columns" -> {{Center}}, "Rows" -> {{Baseline}}}, GridBoxSpacings->{"Columns" -> { Offset[0.27999999999999997`], { Offset[0.7]}, Offset[0.27999999999999997`]}, "Rows" -> { Offset[0.2], { Offset[0.4]}, Offset[0.2]}}], "", ")"}], Function[BoxForm`e$, MatrixForm[BoxForm`e$]]]\);

Here is the merit function:

In[11]:=
F = -Total[M1*M2[[#, #]], 2] &;

The absolute maximum of f is known:

In[12]:=
R = {5, 11, 20, 15, 22, 2, 25, 8, 9, 1, 18, 16, 3, 6, 19, 24, 21, 14, 7, 10, 17, 12, 4, 23, 13}; F[R]
Out[12]=

Monte Carlo search for the maximum (with no guarantee of maximality):

In[13]:=
SeedRandom[1234]; ResourceFunction["MaximizeOverPermutations"][F, S, Method -> {"MonteCarlo", "Iterations" -> 10^6, "AnnealingParameter" -> {0.0001, 1}}]
Out[13]=

Find the extremum of the NUG30 standard Quadratic Assignment Problem:

In[14]:=
(* Evaluate this cell to get the example input *) CloudGet["https://www.wolframcloud.com/obj/29a19086-4dc6-4776-84e0-9c20a84f5c37"]

Here is the merit function:

In[15]:=
F = -Total[M1*M2[[#, #]], 2] &;

The absolute maximum of f is known:

In[16]:=
R = {5, 12, 6, 13, 2, 21, 26, 24, 10, 9, 29, 28, 17, 1, 8, 7, 19, 25, 23, 22, 11, 16, 30, 4, 15, 18, 27, 3, 14, 20}; F[R]
Out[16]=

Monte Carlo search for the maximum (with no guarantee of maximality):

In[17]:=
SeedRandom[1234]; ResourceFunction["MaximizeOverPermutations"][F, S, Method -> {"MonteCarlo", "Iterations" -> 10^7, "AnnealingParameter" -> {0.001, 1}}]
Out[17]=

Publisher

Roman Schmied

Version History

  • 1.0.0 – 12 July 2019

Related Resources

As implemented, the Monte-Carlo algorithm traverses permutation space solely by applying random 2-cycles (that is, by interchanging two random elements at a time). Even though there are n! possible permutations of n objects, any pair of permutations is connected by only n-1 such 2-cycles: permutation space is very well connected and surprisingly small in diameter. For this reason the present implementation does not use larger random cycles to traverse permutation space, and we believe that even highly inhomogeneous maximization problems can be addressed with the present implementation by starting with a sufficiently low annealing parameter β.

License Information