# Optimization of function taking a permutation

I have a function $$f$$ which takes a permutation $$P$$ of the integers 1-100 in order to give a numerical value $$f(x)$$. The function is given by a black box, but is relatively "smooth", so should be amenable to optimization.

For instance, define a function $$f$$:

f[samp_?ListQ] := Total@Total@Table[Table[(-1)^(i), {i, 1, Length[samp]}]
* Reverse@Cos[Mod[samp, n]]* Mod[samp, n], {n, {3, 5, 7, 11, 13, 17, 23}}]


Now f[RandomSample[Range]] will give a numerical value, but I can't figure out how to specify this as an optimization problem only on $$P$$. I can't cast it into the form of the Travelling Salesman Problem, as the function depends on $$x$$ more generally than pairwise interactions.

Edit I mentioned in a comment that what I'm actually trying to do is find the best-scoring set of words in a line of Scrabble tiles as detailed in this puzzle. For this, this is the code for scoring a permutation (without the blanks):

nonblanks =
Sort@ToLowerCase@StringSplit[
"eeeeeeeeeeeeaaaaaaaaaiiiiiiiiioooooooonnnnnnrrrrrrttttttllllssssuuuuddddgggbbccmmppffhhvvwwyykjxqz", ""];

dictionary = Import["https://norvig.com/ngrams/enable1.txt", "List"];
dictionaryMax = Max[StringLength /@ dictionary];

pointSub = Thread[CharacterRange["a", "z"] -> {1, 3, 3, 2, 1, 4, 3, 4, 1, 8, 5,
1, 3, 1, 1, 3, 10, 1, 1, 1, 1, 4, 4, 8, 4, 10}];
score[wordlist_?ListQ] := Total[Flatten@Characters@wordlist /. pointSub];

getScore[samp_?ListQ, scoreOnly_: False] := getScore[samp, scoreOnly] =
Module[{perm, poswords, wordlist},
perm = nonblanks[[samp]];
poswords = Flatten[Table[StringJoin@perm[[i ;; j]], {i, 1, (Length@perm) - 1},
{j, i + 1, Min[(Length@perm), i + dictionaryMax]}]];
wordlist = Intersection[poswords, dictionary];
If[scoreOnly, score@wordlist, {StringJoin@perm, score@wordlist, wordlist}]
]


So given any permutation of hte integers 1-98, getScore will give a numerical value:

getScore[Range]
(* 158 *)


and you can see the words by:

getScore[Range, False]
{"rqciorwlstrndziimdfnsobtroaanikhijxieeevgesiwtpenuoustaearavhnfcdyoa\
glareiuumaploindteeaoeleetogyb", 158, {"aa", "ae", "ag", "aglare",
"an", "ani", "ar", "are", "ear", "el", "en", "es", "et", "glare",
"hi", "in", "khi", "la", "lar", "lee", "leet", "lo", "loin", "ma",
"map", "nu", "oe", "or", "oust", "pe", "pen", "re", "rei", "si",
"so", "sob", "ta", "tae", "tee", "to", "tog", "um", "us", "xi", "yo"}}

• I show some approaches for this in notes here, in particular the sections "Relative position indexing" for permutations" and "Representing and using permutations as shuffles", and "Another shuffle method". – Daniel Lichtblau Jul 1 '19 at 21:07
• Did any of the answers satisfied your need? There are things to do after your question is answered. But wait! It's a good idea to stay vigilant for some time, better approaches may come later improving over previous replies. Experienced users may point alternatives, caveats or limitations. New users should test answers before voting and wait 24 hours before accepting the best one. Participation is essential for the site, please do your part. – rhermans Jul 8 '19 at 10:58
• Should not the call to getScore above include a random permutation such as getScore[RandomSample[Range,False]? – Dominic Jul 10 '19 at 10:40

How about a Monte-Carlo-Metropolis search? I'll implement a simplistic version here. See complete universal code further down. Update: Cleaned-up code now available in the Wolfram Function Repository, so you can use ResourceFunction["MaximizeOverPermutations"] instead of a locally-defined MaximizeOverPermutations. NUG25 and NUG30 are given as applications in the documentation.

To move stochastically through permutation space, we need a random-move generator. Here I'll only use random two-permutations on M=100 list elements: given a list L of 100 elements, generate a new list that has two random elements interchanged,

M = 100;
randomperm[L_] := Permute[L, Cycles[{RandomSample[Range[M], 2]}]]


With this randomperm function we then travel stochastically through permutation-space using the Metropolis-Hastings algorithm. One step of this algorithm consists of proposing a step (with randomperm) and accepting/rejecting it depending on how much the merit function f increases/decreases:

f[samp_?ListQ] := f[samp] =     (* merit function with memoization *)
Total@Total@Table[Table[(-1)^(i), {i, 1, Length[samp]}]*
Reverse@Cos[Mod[samp, n]]*
Mod[samp, n], {n, {3, 5, 7, 11, 13, 17, 23}}]

MH[L_, β_] := Module[{L1, f0, f1, fdiff, prob},
L1 = randomperm[L];   (* proposed new position *)
f0 = f[L];            (* merit function of old position *)
f1 = f[L1];           (* merit function of proposed new position *)
fdiff = N[f1 - f0];   (* probability of accepting the move *)
prob = If[fdiff > 0, 1, E^(β*fdiff)];  (* this is Metropolis-Hastings *)
(* make the move? with calculated probability *)
If[RandomReal[] <= prob, L1, L]]


The parameter β is an effective temperature that nobody knows how to set.

Let's experiment: start with the uniform permutation Range[M] and try with β=1 to see how high we can go with f:

With[{β = 1, nstep = 30000},
Z = NestList[MH[#, β] &, Range[M], nstep];]
ZZ = {#, f[#]} & /@ Z;
ListPlot[ZZ[[All, 2]]] After only $$30\,000$$ Metropolis-Hastings steps we have already found a permutation that gives $$f=1766.64$$:

MaximalBy[ZZ, N@*Last] // DeleteDuplicates
(* {{{69, 31, 91, 2, 47, 89, 75, 37, 96, 61, 40, 22, 64, 95, 81,
10, 66, 43, 19, 82, 85, 26, 28, 62, 78, 72, 34, 54, 45, 86,
57, 60, 65, 33, 13, 74, 5, 8, 11, 68, 77, 88, 23, 15, 35,
50, 83, 3, 93, 9, 18, 53, 63, 4, 58, 56, 30, 42, 46, 55, 36,
94, 1, 87, 51, 44, 14, 21, 97, 27, 52, 49, 99, 73, 39, 71,
7, 20, 41, 48, 24, 38, 29, 84, 6, 79, 90, 16, 59, 32, 12,
70, 98, 67, 92, 100, 76, 25, 17, 80},
184 + 154 Cos - 157 Cos - 252 Cos - 194 Cos +
69 Cos + 238 Cos + 190 Cos + 8 Cos - 154 Cos -
120 Cos + 17 Cos + 94 Cos + 134 Cos + 19 Cos -
81 Cos - 76 Cos + 14 Cos + 23 Cos + 36 Cos +
4 Cos - 35 Cos - 21 Cos}} *)


We can continue along this line with (i) increasing $$\beta$$, and (ii) introducing more moves, apart from randomperm.

For example, we can raise $$\beta$$ slowly during the MH-Iteration, starting with $$\beta_{\text{min}}$$ and going up to $$\beta_{\text{max}}$$: this gives a simulated annealing advantage and tends to give higher results for f.

With[{βmin = 10^-2, βmax = 10, nstep = 10^6},
With[{γ = N[(βmax/βmin)^(1/nstep)]},
Z = NestList[{MH[#[], #[]], γ*#[]} &, {Range[M], βmin}, nstep];]]
ZZ = {#[], #[], f[#[]]} & /@ Z;
ListLogLinearPlot[ZZ[[All, {2, 3}]]] After playing around for a while, all f-values computed so far are stored as DownValues of f and we can easily determine the absolutely largest f-value seen so far: in my case, the largest value ever seen was $$f=1805.05$$,

MaximalBy[Cases[DownValues[f],
RuleDelayed[_[f[L_ /; VectorQ[L, NumericQ]]], g_] :> {L, g}],
N@*Last]
(* {{{93, 61, 1, 15, 7, 2, 51, 72, 92, 78, 59, 43, 58, 10, 63, 21, 13,
48, 76, 49, 99, 42, 35, 31, 11, 95, 69, 88, 82, 36, 57, 77, 97, 73,
47, 9, 28, 86, 24, 79, 6, 71, 39, 27, 83, 68, 40, 33, 98, 80, 75,
37, 91, 32, 19, 3, 56, 25, 84, 87, 41, 100, 52, 20, 64, 67, 34, 60,
14, 50, 70, 16, 46, 17, 90, 94, 5, 55, 23, 54, 45, 4, 85, 38, 65,
26, 18, 44, 29, 22, 81, 89, 66, 74, 96, 62, 30, 8, 12, 53},
170 + 174 Cos - 150 Cos - 282 Cos - 172 Cos +
120 Cos + 218 Cos + 191 Cos - 13 Cos - 214 Cos -
141 Cos + 22 Cos + 117 Cos + 109 Cos +
27 Cos - 60 Cos - 52 Cos + 6 Cos + 23 Cos +
43 Cos - 8 Cos - 29 Cos - 19 Cos}} *)

%[[All, 2]] // N
(* {1805.05} *)


# Complete and universal code for permutational optimization

Here is a version of the above code that is more cleaned up and emits useful error messages:

(* error messages *)
MaximizeOverPermutations::Pstart = "Starting permutation 1 is invalid.";
MaximizeOverPermutations::f = "Optimization function does not yield a real number on 1.";

(* interface for calculation at fixed β *)
MaximizeOverPermutations[f_,                        (* function to optimize *)
M_Integer /; M >= 2,       (* number of arguments of f *)
β_?NumericQ,               (* annealing parameter *)
steps_Integer?Positive,    (* number of iteration steps *)
Pstart_: Automatic] :=     (* starting permutation *)
MaximizeOverPermutations[f, M, {β, β}, steps, Pstart]

(* interface for calculation with geometrically ramping β *)
MaximizeOverPermutations[f_,                        (* function to optimize *)
M_Integer /; M >= 2,       (* number of arguments of f *)
{βstart_?NumericQ,         (* annealing parameter at start *)
βend_?NumericQ},          (* annealing parameter at end *)
steps_Integer?Positive,    (* number of iteration steps *)
Pstart_: Automatic] :=     (* starting permutation *)
Module[{P, g, Pmax, gmax, Pnew, gnew, β, γ, prob},
(* determine the starting permutation *)
P = Which[Pstart === Automatic, Range[M],
VectorQ[Pstart, IntegerQ] && Sort[Pstart] == Range[M], Pstart,
True, Message[MaximizeOverPermutations::Pstart, Pstart]; $$Failed]; If[FailureQ[P], Return[$$Failed]];
(* evaluate the function on the starting permutation *)
g = f[P] // N;
If[! Element[g, Reals], Message[MaximizeOverPermutations::f, P]; Return[$$Failed]]; (* store maximum merit function *) Pmax = P; gmax = g; (* inverse temperature: geometric progression from βstart to βend *) β = βstart // N; γ = (βend/βstart)^(1/(steps - 1)) // N; (* Metropolis-Hastings iteration *) Do[ (* propose a new permutation by applying a random 2-cycle *) Pnew = Permute[P, Cycles[{RandomSample[Range[M], 2]}]]; (* evaluate the function on the new permutation *) gnew = f[Pnew] // N; If[! Element[gnew, Reals], Message[MaximizeOverPermutations::f, Pnew]; Return[$$Failed]];
(* Metropolis-Hasting acceptance probability *)
prob = If[gnew > g, 1, Quiet[Exp[-β (g - gnew)], General::munfl]];
(* acceptance/rejection of the new permutation *)
If[RandomReal[] <= prob,
P = Pnew; g = gnew;
If[g > gmax, Pmax = P; gmax = g]];
(* update inverse temperature *)
β *= γ,
{steps}];
(* return maximum found *)
{Pmax, gmax}]


The OP's problem can be optimized with

f[samp_List] := Total[Table[(-1)^Range[Length[samp]]*Reverse@Cos[Mod[samp, n]]*
Mod[samp, n], {n, {3, 5, 7, 11, 13, 17, 23}}], 2]
MaximizeOverPermutations[f, 100, {1/100, 10}, 10^6]


A simpler problem, where we know the perfect optimum, is

SeedRandom;
MM = 100;
x = RandomVariate[NormalDistribution[], MM];
Z[L_List] := L.x


The optimum is known: put the permutation in the same order as the numbers in the list x. For this particular case of random numbers, we get

Z[Ordering[Ordering[x]]]
(*    2625.98    *)


A quick search yields something not quite as high,

MaximizeOverPermutations[Z, MM, 1, 10^4][]
(*    2597.67    *)


To track the progress of the Monte-Carlo search, use a Sow/Reap combination:

zz = Reap[MaximizeOverPermutations[Sow@*Z, MM, 1, 10^4]];
ListPlot[zz[[2, 1]], GridLines -> {None, {zz[[1, 2]]}}] zz = Reap[MaximizeOverPermutations[Sow@*Z, MM, {1/10, 10}, 10^5]];
ListPlot[zz[[2, 1]], GridLines -> {None, {zz[[1, 2]]}}] • Very nice, and it deserves more upvotes (mine was the second such). You should consider submitting this to the Wolfram Function Repository as a general optimization code for permutation problems. Inputs would be the objective function f and permutation length n, with underlying assumption that f works with permutations of Range[n]. Plus any other inputs you might want (temperatiure function, quenching rate, whatever). – Daniel Lichtblau Jul 3 '19 at 16:31
• Thanks for the kind words, @DanielLichtblau . It works, yes; and yet right now my construction feels a bit amateurish to me, because I have no intuition for the topology of the space of permutations and how one travels efficiently in this space (maybe with steps that aren't 2-cycles?). I'll be thinking and researching some more. – Roman Jul 3 '19 at 17:08
• Upon a bit more reflection it seems that it takes only $M-1$ applications of a 2-cycle (a permutation of two indices) to go from any permutation of length $M$ to any other. So the method proposed above may not be bad after all, as it can cross permutation space rather quickly despite the space's incredible vastness. – Roman Jul 3 '19 at 18:13
• Thanks Roman. You can very easily apply n-cycles with: randomperm[L_,n_] := Permute[L, Cycles[{RandomSample[Range[Length@L], n]}]], and then I'm trying drawing the $n$ from a negative binomial distribution, so it weights to small shuffles but sometimes tries large ones. – KraZug Jul 3 '19 at 20:10
• I used your annealing variant on the NUG25 problem that has been in benchmarks for a few decades (I show it in the tutorial notes I pointed to in a comment). Here's the thing. Your method found a global optimum* within 10^5 steps (less than half a minute). Asterisk: It won't find the global every time for that number of steps. But it will come close. This is quite powerful. I'd really like to see this in the WFR (and maybe some day in the Wolfram Language proper). – Daniel Lichtblau Jul 9 '19 at 17:11

Here is one approach from among the ones I allude to in a comment.

f[samp_?ListQ] :=
Total@Total@
Table[Table[(-1)^(i), {i, 1, Length[samp]}]*
Reverse@Cos[Mod[samp, n]]*
Mod[samp, n], {n, {3, 5, 7, 11, 13, 17, 23}}]


Now just define a function that takes a numeric vector, creates a permutation, and evaluates f on it.

g[ll : {_?NumberQ ..}] := N[f[Ordering[ll]]]


We can get a reasonable value with NMaximize. Restricting the range of the values seems to help here.

n = 100;
vars = Array[x, n];
AbsoluteTiming[{max, vals} =
NMaximize[{g[vars], Thread[0 <= vars <= 1]},
Map[{#, 0, 1} &, vars], MaxIterations -> 5000];]
max
best = Ordering[vars /. vals]
N[f[best]]

(* During evaluation of In:= NMaximize::cvmit: Failed to converge to the requested accuracy or precision within 5000 iterations.

Out= {62.699518, Null}

Out= 636.619153268

Out= {9, 40, 46, 2, 19, 47, 53, 77, 97, 87, 21, 33, 71, 35, 95, \
73, 39, 28, 52, 43, 6, 75, 5, 20, 27, 31, 22, 64, 49, 83, 42, 38, 92, \
58, 65, 79, 30, 11, 12, 13, 7, 66, 86, 67, 41, 4, 72, 100, 60, 10, 1, \
48, 81, 8, 84, 55, 36, 32, 25, 96, 70, 44, 80, 16, 18, 68, 29, 88, \
89, 15, 91, 69, 23, 17, 82, 90, 94, 93, 50, 99, 59, 85, 74, 62, 56, \
26, 24, 34, 78, 3, 98, 63, 14, 61, 51, 76, 45, 54, 37, 57}

Out= 636.619153268 *)


Could of course instead minimize in the same manner. Also there are numerous variations one might try, using option and method sub-option settings for NMinimize.

• Ahh, that is such a simple way to turn the permutation problem into something that the built-in functions can work on. – KraZug Jul 3 '19 at 20:28

it seems that Objective Function must return Numeric Value,not Symbolic expression.

f[samp_?ListQ] :=
Total@Total@
Table[Table[(-1)^(i), {i, 1, Length[samp]}]*
Reverse@Cos[Mod[samp, n]]*
Mod[samp, n], {n, {3, 5, 7, 11, 13, 17, 23}}]

Nf[samp_?ListQ] :=
N@Total@Total@
Table[Table[(-1)^(i), {i, 1, Length[samp]}]*
Reverse@Cos[Mod[samp, n]]*
Mod[samp, n], {n, {3, 5, 7, 11, 13, 17, 23}}]


Print[forwardDP[f, Range[1, 100]] // f // N]


-118.075

Print[forwardDP[Nf, Range[1, 100]] // Nf]


1164.08

The first thing that came to mind is the heuristic.
The other is approximated dynamic programming.

### Heuristic

Easy and Fast Heuristic Implementation.

Table[
Nest[
With[{try = RandomSample@Range},
tryvalue = f[try];
If[#2 >= tryvalue, {#1, #2},
{try, tryvalue}]] & @@ # &,
{1, -10000}, 500],
{100}
] // MaximalBy[#, #[] &] & // Flatten[#, 1] & (*no elements should be duplicate.*)
Not@*Equal @@ # & /@ Subsets[First@%, {2}] // And @@ # &


=>

True

### Dynamic Programming(forward)

forwardDP[obj_, action_?(VectorQ[#, IntegerQ] &)] :=
Block[{solution, nothing, tryaction},
solution = ConstantArray[nothing, Length@action];
Do[solution[[index]] = First[First[Table[solution[[index]] = trynum;
tryaction =
Join[DeleteCases[solution, nothing],
DeleteCases[action, x_ /; ContainsAny[solution][{x}]]];
{trynum, obj[tryaction]}, {trynum,
DeleteCases[action,
x_ /; ContainsAny[DeleteCases[solution, nothing]][{x}]]}] //
MaximalBy[#, #[] &] &]], {index, Range[1, Length@action]}];
solution];

forwardDP[f, Range[1, 100]] // AbsoluteTiming f[%]


=>

608

Not@*Equal @@ # & /@ Subsets[%%, {2}] // And @@ # &


=>

True

About feasible region of control/action,please modify the code around DeleteCases of trynum and tryaction for your problem.

• Can you explain what the dynamic programming code is doing? – KraZug Jul 2 '19 at 7:47
• for each time-step,from index=1 to index=100,We tried all possible solution from feasible region determined by the constraints and take the best one. – Xminer Jul 2 '19 at 16:15
• First re-write the solution for that time step,index with trynum, then create a candidate solutiontryaction that satisfies the constraints. – Xminer Jul 2 '19 at 16:20
• please note when we choose the best action,sometime MaximalBy returns more than two values whose scores are the same depending on the object function. in that case tweak the constraint. – Xminer Jul 2 '19 at 16:26
• Finally,this is approximated dynamic programming,there may be better optimal solution,but this is efficiency more than heuristic or NMaximize: ) – Xminer Jul 2 '19 at 16:30

Code to include the blank tiles:

Revised 7/26/19: (previous code did not include definition of cRange--added it). Also converted the addition of the two blank tiles to a function.

We have 98 lettered tiles and two blanks. We first compile a list of all possible combinations of 2-letters for the blanks. Run metropolis with M=100 for each combination. That gives 351 runs. Find the maximum from that set. Here is the code to incorporate the blanks in the list:

    cRange = CharacterRange["a", "z"]
theBlanks =
Join[Subsets[CharacterRange["a", "z"], {2}], {#, #} & /@ cRange];

scrabbleList[n_] :=
"eeeeeeeeeeeeaaaaaaaaaiiiiiiiiioooooooonnnnnnrrrrrrttttttllllssssuuuu\
ddddgggbbccmmppffhhvvwwyykjxqz" <> theBlanks[[n]]


Note: See Solution to scrabble puzzle for a 1629 score using the letter "S" twice.

• Was this intended as a response to a different question? – Daniel Lichtblau Jul 10 '19 at 14:04
• No. Just found the problem interesting. Please delete if off-topic. – Dominic Jul 10 '19 at 14:12
• I gather now that this is also intended to be a permutation optimization problem. Probably it should be a separate thread (linking to this one). But maybe the moderators will be okay with having it here. It does need a full specification however. I do not know what cRange is, or what specifically is to be optimized, so Mathematica code for the objective function would need to be added. – Daniel Lichtblau Jul 11 '19 at 16:05
• Hi Daniel, the OP stated the problem: fivethirtyeight.com/features/whats-your-best-scrabble-string. Essentially find the combination of 100 scrabble letters with the most points (2 blanks included). Roman above basically wrote the code to do the problem with the 98 pieces and I was just interested in how it would be done with the two blank pieced added. If I had time I would work it through. Maybe sometime in the future. It's a beautiful problem I think. – Dominic Jul 11 '19 at 18:16
• @DanielLichtblau, yes, I updated my question explaining that was my actual objective function. Thank you for giving the bounty btw, Roman's answer definitely deserves it. I have considered making that a separate question. – KraZug Jul 11 '19 at 18:42