Optimization of function taking a permutation - Mathematica Stack Exchange most recent 30 from mathematica.stackexchange.com 2019-11-11T23:24:12Z https://mathematica.stackexchange.com/feeds/question/201333 https://creativecommons.org/licenses/by-sa/4.0/rdf https://mathematica.stackexchange.com/q/201333 13 Optimization of function taking a permutation KraZug https://mathematica.stackexchange.com/users/23105 2019-07-01T20:11:37Z 2019-07-26T20:06:06Z <p>I have a function <span class="math-container">$f$</span> which takes a permutation <span class="math-container">$P$</span> of the integers 1-100 in order to give a numerical value <span class="math-container">$f(x)$</span>. The function is given by a black box, but is relatively "smooth", so should be amenable to optimization.</p> <p>For instance, define a function <span class="math-container">$f$</span>:</p> <pre><code>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}}] </code></pre> <p>Now <code>f[RandomSample[Range]]</code> will give a numerical value, but I can't figure out how to specify this as an optimization problem only on <span class="math-container">$P$</span>. I can't cast it into the form of the Travelling Salesman Problem, as the function depends on <span class="math-container">$x$</span> more generally than pairwise interactions.</p> <p>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 <a href="https://fivethirtyeight.com/features/whats-your-best-scrabble-string/" rel="noreferrer">as detailed in this puzzle</a>. For this, this is the code for scoring a permutation (without the blanks):</p> <pre><code>nonblanks = Sort@ToLowerCase@StringSplit[ "eeeeeeeeeeeeaaaaaaaaaiiiiiiiiioooooooonnnnnnrrrrrrttttttllllssssuuuuddddgggbbccmmppffhhvvwwyykjxqz", ""]; dictionary = Import["https://norvig.com/ngrams/enable1.txt", "List"]; dictionaryMax = Max[StringLength /@ dictionary]; pointSub = Thread[CharacterRange["a", "z"] -&gt; {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}] ] </code></pre> <p>So given any permutation of hte integers 1-98, getScore will give a numerical value:</p> <pre><code>getScore[Range] (* 158 *) </code></pre> <p>and you can see the words by:</p> <pre><code>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"}} </code></pre> https://mathematica.stackexchange.com/questions/201333/-/201341#201341 3 Answer by Xminer for Optimization of function taking a permutation Xminer https://mathematica.stackexchange.com/users/61541 2019-07-02T00:57:25Z 2019-07-02T20:48:11Z <p>it seems that Objective Function must return Numeric Value,not Symbolic expression.</p> <pre><code>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}}] </code></pre> <p><br></p> <pre><code>Print[forwardDP[f, Range[1, 100]] // f // N] </code></pre> <blockquote> <p>-118.075</p> </blockquote> <pre><code>Print[forwardDP[Nf, Range[1, 100]] // Nf] </code></pre> <blockquote> <p>1164.08</p> </blockquote> <hr> <p>The first thing that came to mind is the heuristic.<br> The other is approximated dynamic programming.</p> <h3><em>Heuristic</em></h3> <p>Easy and Fast Heuristic Implementation.</p> <pre><code>Table[ Nest[ With[{try = RandomSample@Range}, tryvalue = f[try]; If[#2 &gt;= tryvalue, {#1, #2}, {try, tryvalue}]] &amp; @@ # &amp;, {1, -10000}, 500], {100} ] // MaximalBy[#, #[] &amp;] &amp; // Flatten[#, 1] &amp; </code></pre> <p><img src="https://i.stack.imgur.com/2N67r.png" alt="Mathematica graphics"></p> <pre><code>(*no elements should be duplicate.*) Not@*Equal @@ # &amp; /@ Subsets[First@%, {2}] // And @@ # &amp; </code></pre> <p>=> </p> <blockquote> <p>True</p> </blockquote> <h3><em>Dynamic Programming(forward)</em></h3> <pre><code>forwardDP[obj_, action_?(VectorQ[#, IntegerQ] &amp;)] := 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[#, #[] &amp;] &amp;]], {index, Range[1, Length@action]}]; solution]; forwardDP[f, Range[1, 100]] // AbsoluteTiming </code></pre> <p><img src="https://i.stack.imgur.com/mmo59.png" alt="Mathematica graphics"></p> <pre><code>f[%] </code></pre> <p>=></p> <blockquote> <pre><code>608 </code></pre> </blockquote> <pre><code>Not@*Equal @@ # &amp; /@ Subsets[%%, {2}] // And @@ # &amp; </code></pre> <p>=></p> <blockquote> <p>True</p> </blockquote> <p>About feasible region of control/action,please modify the code around <code>DeleteCases</code> of <code>trynum</code> and <code>tryaction</code> for your problem.</p> https://mathematica.stackexchange.com/questions/201333/-/201423#201423 22 Answer by Roman for Optimization of function taking a permutation Roman https://mathematica.stackexchange.com/users/26598 2019-07-02T21:57:43Z 2019-07-13T05:20:04Z <p>How about a <a href="https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm" rel="noreferrer">Monte-Carlo-Metropolis</a> search? I'll implement a simplistic version here. See complete universal code further down. <strong>Update:</strong> Cleaned-up code now <a href="https://resources.wolframcloud.com/FunctionRepository/resources/MaximizeOverPermutations" rel="noreferrer">available in the Wolfram Function Repository</a>, so you can use <code>ResourceFunction["MaximizeOverPermutations"]</code> instead of a locally-defined <code>MaximizeOverPermutations</code>. NUG25 and NUG30 are given as applications in the documentation.</p> <p>To move stochastically through permutation space, we need a random-move generator. Here I'll only use random two-permutations on <code>M=100</code> list elements: given a list <code>L</code> of 100 elements, generate a new list that has two random elements interchanged,</p> <pre><code>M = 100; randomperm[L_] := Permute[L, Cycles[{RandomSample[Range[M], 2]}]] </code></pre> <p>With this <code>randomperm</code> function we then travel stochastically through permutation-space using the <a href="https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm" rel="noreferrer">Metropolis-Hastings algorithm</a>. One step of this algorithm consists of proposing a step (with <code>randomperm</code>) and accepting/rejecting it depending on how much the merit function <code>f</code> increases/decreases:</p> <pre><code>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 &gt; 0, 1, E^(β*fdiff)]; (* this is Metropolis-Hastings *) (* make the move? with calculated probability *) If[RandomReal[] &lt;= prob, L1, L]] </code></pre> <p>The parameter <code>β</code> is an effective temperature that nobody knows how to set.</p> <p>Let's experiment: start with the uniform permutation <code>Range[M]</code> and try with <code>β=1</code> to see how high we can go with <code>f</code>:</p> <pre><code>With[{β = 1, nstep = 30000}, Z = NestList[MH[#, β] &amp;, Range[M], nstep];] ZZ = {#, f[#]} &amp; /@ Z; ListPlot[ZZ[[All, 2]]] </code></pre> <p><a href="https://i.stack.imgur.com/VNdNk.png" rel="noreferrer"><img src="https://i.stack.imgur.com/VNdNk.png" alt="enter image description here"></a></p> <p>After only <span class="math-container">$30\,000$</span> Metropolis-Hastings steps we have already found a permutation that gives <span class="math-container">$f=1766.64$</span>:</p> <pre><code>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}} *) </code></pre> <p>We can continue along this line with (i) increasing <span class="math-container">$\beta$</span>, and (ii) introducing more moves, apart from <code>randomperm</code>.</p> <p>For example, we can raise <span class="math-container">$\beta$</span> slowly during the MH-Iteration, starting with <span class="math-container">$\beta_{\text{min}}$</span> and going up to <span class="math-container">$\beta_{\text{max}}$</span>: this gives a <a href="https://en.wikipedia.org/wiki/Simulated_annealing" rel="noreferrer">simulated annealing</a> advantage and tends to give higher results for <code>f</code>.</p> <pre><code>With[{βmin = 10^-2, βmax = 10, nstep = 10^6}, With[{γ = N[(βmax/βmin)^(1/nstep)]}, Z = NestList[{MH[#[], #[]], γ*#[]} &amp;, {Range[M], βmin}, nstep];]] ZZ = {#[], #[], f[#[]]} &amp; /@ Z; ListLogLinearPlot[ZZ[[All, {2, 3}]]] </code></pre> <p><a href="https://i.stack.imgur.com/tAGl4.png" rel="noreferrer"><img src="https://i.stack.imgur.com/tAGl4.png" alt="enter image description here"></a></p> <p>After playing around for a while, all <code>f</code>-values computed so far are stored as <code>DownValues</code> of <code>f</code> and we can easily determine the absolutely largest <code>f</code>-value seen so far: in my case, the largest value ever seen was <span class="math-container">$f=1805.05$</span>,</p> <pre><code>MaximalBy[Cases[DownValues[f], RuleDelayed[_[f[L_ /; VectorQ[L, NumericQ]]], g_] :&gt; {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} *) </code></pre> <h1>Complete and universal code for permutational optimization</h1> <p>Here is a version of the above code that is more cleaned up and emits useful error messages:</p> <pre><code>(* 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 &gt;= 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 &gt;= 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] &amp;&amp; Sort[Pstart] == Range[M], Pstart, True, Message[MaximizeOverPermutations::Pstart, Pstart]; <span class="math-container">$Failed]; If[FailureQ[P], Return[$</span>Failed]]; (* evaluate the function on the starting permutation *) g = f[P] // N; If[! Element[g, Reals], Message[MaximizeOverPermutations::f, P]; Return[<span class="math-container">$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[$</span>Failed]]; (* Metropolis-Hasting acceptance probability *) prob = If[gnew &gt; g, 1, Quiet[Exp[-β (g - gnew)], General::munfl]]; (* acceptance/rejection of the new permutation *) If[RandomReal[] &lt;= prob, P = Pnew; g = gnew; If[g &gt; gmax, Pmax = P; gmax = g]]; (* update inverse temperature *) β *= γ, {steps}]; (* return maximum found *) {Pmax, gmax}] </code></pre> <p>The OP's problem can be optimized with</p> <pre><code>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] </code></pre> <p>A simpler problem, where we know the perfect optimum, is</p> <pre><code>SeedRandom; MM = 100; x = RandomVariate[NormalDistribution[], MM]; Z[L_List] := L.x </code></pre> <p>The optimum is known: put the permutation in the <a href="https://mathematica.stackexchange.com/q/194094">same order as the numbers in the list <code>x</code></a>. For this particular case of random numbers, we get</p> <pre><code>Z[Ordering[Ordering[x]]] (* 2625.98 *) </code></pre> <p>A quick search yields something not quite as high,</p> <pre><code>MaximizeOverPermutations[Z, MM, 1, 10^4][] (* 2597.67 *) </code></pre> <p>To track the progress of the Monte-Carlo search, use a <a href="https://reference.wolfram.com/language/tutorial/CollectingExpressionsDuringEvaluation.html" rel="noreferrer"><code>Sow</code>/<code>Reap</code> combination</a>:</p> <pre><code>zz = Reap[MaximizeOverPermutations[Sow@*Z, MM, 1, 10^4]]; ListPlot[zz[[2, 1]], GridLines -&gt; {None, {zz[[1, 2]]}}] </code></pre> <p><a href="https://i.stack.imgur.com/JFqGe.png" rel="noreferrer"><img src="https://i.stack.imgur.com/JFqGe.png" alt="enter image description here"></a></p> <pre><code>zz = Reap[MaximizeOverPermutations[Sow@*Z, MM, {1/10, 10}, 10^5]]; ListPlot[zz[[2, 1]], GridLines -&gt; {None, {zz[[1, 2]]}}] </code></pre> <p><a href="https://i.stack.imgur.com/aPTi6.png" rel="noreferrer"><img src="https://i.stack.imgur.com/aPTi6.png" alt="enter image description here"></a></p> https://mathematica.stackexchange.com/questions/201333/-/201424#201424 5 Answer by Daniel Lichtblau for Optimization of function taking a permutation Daniel Lichtblau https://mathematica.stackexchange.com/users/51 2019-07-02T22:13:22Z 2019-07-02T22:13:22Z <p>Here is one approach from among the ones I allude to in a comment.</p> <pre><code>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}}] </code></pre> <p>Now just define a function that takes a numeric vector, creates a permutation, and evaluates <code>f</code> on it.</p> <pre><code>g[ll : {_?NumberQ ..}] := N[f[Ordering[ll]]] </code></pre> <p>We can get a reasonable value with <code>NMaximize</code>. Restricting the range of the values seems to help here.</p> <pre><code>n = 100; vars = Array[x, n]; AbsoluteTiming[{max, vals} = NMaximize[{g[vars], Thread[0 &lt;= vars &lt;= 1]}, Map[{#, 0, 1} &amp;, vars], MaxIterations -&gt; 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 *) </code></pre> <p>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 <code>NMinimize</code>.</p> https://mathematica.stackexchange.com/questions/201333/-/201853#201853 0 Answer by Dominic for Optimization of function taking a permutation Dominic https://mathematica.stackexchange.com/users/47466 2019-07-10T12:00:55Z 2019-07-26T20:06:06Z <p>Code to include the blank tiles: </p> <p>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.</p> <p>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:</p> <pre><code> cRange = CharacterRange["a", "z"] theBlanks = Join[Subsets[CharacterRange["a", "z"], {2}], {#, #} &amp; /@ cRange]; scrabbleList[n_] := "eeeeeeeeeeeeaaaaaaaaaiiiiiiiiioooooooonnnnnnrrrrrrttttttllllssssuuuu\ ddddgggbbccmmppffhhvvwwyykjxqz" &lt;&gt; theBlanks[[n]] </code></pre> <p>Note: See <a href="https://fivethirtyeight.com/features/whats-the-optimal-way-to-play-horse/" rel="nofollow noreferrer">Solution to scrabble puzzle</a> for a 1629 score using the letter "S" twice. </p>