How to use genetic algorithm to calculate the maximum value of this matrix

1-9 the 9 numbers are arranged into 3 * 3 matrix without repetition.How to select the matrix with the largest determinant value:

list = Permutations[Range[9], {9}];
matrix = Partition[#, 3] & /@ list;
matrix[[#]] & /@ pos
Det[%[[1]]]


I've used permutation to get results, but how can I use genetic algorithms to achieve similar results.

• Why not try NMaximum with Method->”DifferentialEvolution” – M.R. Jan 31 '20 at 3:54

Follows a classical Genetic Algorithm with the normal functionalities as

1. Population initialization
2. Fitness Evaluation
3. Crossover
4. Mutation
5. Offspring Selection

The main difficulty found to implement this procedure was the crossover implementation, because the offspring should preserve always the same elements (1,2,3,4,5,6,7,8,9) without absences or repetitions. By this reason, the corrections needed after the one-point crossover, are considered as mutations. The program bulk was extracted from one repository but were included the essential modifications. The script can be utilized for a generic matrix dimension dim.

Clear[doMutation];
doMutation[string_] := Module[{tempstring, i, ind1, ind2, atom}, tempstring = string;
If[Random[] < mutationRate,
ind1 = RandomInteger[{1, length}];
ind2 = RandomInteger[{1, length}];
atom = tempstring[[ind1]];
tempstring[[ind1]] = tempstring[[ind2]];
tempstring[[ind2]] = atom];
Return[tempstring]]

Clear[correct]
correct[lista_] := Module[{out, ok = Table[0, length],
list = lista, i, k, ind},
out = Complement[Range[length], list];
If[Length[out] == 0, Return[list],
For[i = 1; k = 1, i <= length, i++, ind = list[[i]];
If[ok[[ind]] == 0, ok[[ind]] = 1, list[[i]] = out[[k]]; k = k + 1]]];
Return[list]]

Clear[fitnessFunction];
fitnessFunction[list_] := Max[0, Det[ArrayReshape[list, {dim, dim}]]]

Clear[doSingleCrossover];
doSingleCrossover[ string1_, string2_] := Module[{cut, temp1, temp2},
cut = RandomInteger[{1, length}];
temp1 = Join[ Take[string1, cut], Drop[string2, cut] ];
temp2 = Join[ Take[string2, cut], Drop[string1, cut] ];
{correct[temp1], correct[temp2]} ]

Clear[doCumSumOfFitness];
doCumSumOfFitness := Module[{temp},
temp = 0.0;
Table[ temp += popFitness[[i]], {i, popSize} ]]

Clear[doSingleSelection];
doSingleSelection := Module[{rfitness, ind},
rfitness = RandomReal[{0, cumFitness[[popSize]]}];
ind = 1;
While[ rfitness > cumFitness[[ind]], ind++ ];
Return[ind]]

Clear[selectPair];
selectPair := Module[{ind1, ind2},
ind1 = doSingleSelection;
While[ (ind2 = doSingleSelection) == ind1 ];
{ind1, ind2}]

Clear[pickRandomPair];
pickRandomPair := Module[{ind1, ind2},
ind1 = RandomInteger[{1, popSize}];
While[ (ind2 = RandomInteger[ {1, popSize}]) == ind1 ];
{ind1, ind2}]

Clear[exchangeString];
exchangeString[ind_, newstring_, newF_] := Module[{},
popStrings[[ind]] = newstring;
popFitness[[ind]] = newF]

Clear[renormalizeFitness];
renormalizeFitness[fitness0_List] := Module[{minF, maxF, a, b, fitness = fitness0, i},
minF = Min[fitness];
maxF = Max[fitness];
a = 0.5*maxF/(maxF + minF);
b = (1 - a)*maxF;
Map[  a # + b &, fitness]]

Clear[bestDet]
bestDet := Module[{bestFitness = -1, i, ibest = 1},
For[i = 1, i <= popSize, i++,
If[popFitness[[i]] > bestFitness, bestFitness = popFitness[[i]]; ibest = i]];
If[bestFitness > bestOfAll, bestOfAll = bestFitness;
bestIndividual = popStrings[[ibest]]];
Return[popStrings[[ibest]]]]

Clear[doInitialize];
doInitialize := Module[{i},
popFitness = Table[fitnessFunction[popStrings[[i]]], {i, popSize} ];
popFitness = renormalizeFitness[popFitness];
cumFitness = doCumSumOfFitness;
listOfCumFitness = {cumFitness[[popSize]]};
historyOfPop = {bestDet}]

Clear[updateGenerationSync];
updateGenerationSync := Module[{parentsid, children, ip},
parentsid = {};
Do[AppendTo[parentsid, selectPair], {popSize/2}];
children = {};
Do[AppendTo[children,
doSingleCrossover[popStrings[[parentsid[[ip, 1]]]],
popStrings[[parentsid[[ip, 2]]]]]], {ip, popSize/2}];
popStrings = Flatten[ children, 1];
popStrings = Map[doMutation, popStrings];
popFitness = Map[fitnessFunction, popStrings];
popFitness = renormalizeFitness[popFitness];
cumFitness = doCumSumOfFitness]


And now the main program

SeedRandom[4];
bestOfAll = -1;
dim = 6;
length = dim^2;
popSize = 100; (* should be even *)
numberOfEpochs = 500;
mutationRate = 0.005;
popStrings = Table[RandomSample[Range[length], length], {popSize} ];
doInitialize;

Do[updateGenerationSync;
AppendTo[historyOfPop,bestDet];
AppendTo[ listOfCumFitness, cumFitness[[popSize]] ],
{numberOfEpochs} ];

ListLinePlot[listOfCumFitness, PlotRange -> All ]


ListLinePlot[Map[fitnessFunction, historyOfPop], PlotRange -> All]


bestIndividual
fitnessFunction[bestIndividual]

(*{27, 11, 36, 29, 6, 14, 23, 16, 22, 4, 34, 10, 18, 33, 1, 32, 13, 8, 31, 3, 2, 9, 12, 25, 5, 7, 20, 26, 19, 28, 21, 35, 24, 17, 15, 30}*)

(*1181916347*)


NOTE

This script can be enhanced including elitism.

• Thank you very much. How can you modify this algorithm to quickly find the maximal determinant in the case of 6×6 matrix using the integers 1 to 36 quickly ? oeis.org/… – A little mouse on the pampas Feb 1 '20 at 7:05
• The script now is general for any dimension dim – Cesareo Feb 1 '20 at 12:14

Using the function MaximizeOverPermutations from this great answer by Roman:

ClearAll[f]
f[samp_List] := Det[Partition[samp, 3]]
MaximizeOverPermutations[f, 9, {1/100, 10}, 10^4] // AbsoluteTiming


{0.664876, {{5,7,1,3,6,8,9,2,4}, 412.}}

ResourceFunction["MaximizeOverPermutations"][f, 9]


{{{1, 4, 8, 7, 2, 6, 5, 9, 3}, {1, 5, 7, 8, 3, 6, 4, 9, 2}, {1, 7, 5, 4, 2, 9, 8, 6, 3}, {1, 8, 4, 5, 3, 9, 7, 6, 2}, {2, 4, 9, 7, 1, 5, 6, 8, 3}, {2, 6, 7, 9, 3, 5, 4, 8, 1}, {2, 7, 6, 4, 1, 8, 9, 5, 3}, {2, 9, 4, 6, 3, 8, 7, 5, 1}, {3, 5, 9, 8, 1, 4, 6, 7, 2}, {3, 6, 8, 9, 2, 4, 5, 7, 1}, {3, 8, 6, 5, 1, 7, 9, 4, 2}, {3, 9, 5, 6, 2, 7, 8, 4, 1}, {4, 1, 8, 9, 5, 3, 2, 7, 6}, {4, 2, 9, 8, 6, 3, 1, 7, 5}, {4, 8, 1, 2, 6, 7, 9, 3, 5}, {4, 9, 2, 1, 5, 7, 8, 3, 6}, {5, 1, 7, 9, 4, 2, 3, 8, 6}, {5, 3, 9, 7, 6, 2, 1, 8, 4}, {5, 7, 1, 3, 6, 8, 9, 2, 4}, {5, 9, 3, 1, 4, 8, 7, 2, 6}, {6, 2, 7, 8, 4, 1, 3, 9, 5}, {6, 3, 8, 7, 5, 1, 2, 9, 4}, {6, 7, 2, 3, 5, 9, 8, 1, 4}, {6, 8, 3, 2, 4, 9, 7, 1, 5}, {7, 1, 5, 6, 8, 3, 2, 4, 9}, {7, 2, 6, 5, 9, 3, 1, 4, 8}, {7, 5, 1, 2, 9, 4, 6, 3, 8}, {7, 6, 2, 1, 8, 4, 5, 3, 9}, {8, 1, 4, 6, 7, 2, 3, 5, 9}, {8, 3, 6, 4, 9, 2, 1, 5, 7}, {8, 4, 1, 3, 9, 5, 6, 2, 7}, {8, 6, 3, 1, 7, 5, 4, 2, 9}, {9, 2, 4, 5, 7, 1, 3, 6, 8}, {9, 3, 5, 4, 8, 1, 2, 6, 7}, {9, 4, 2, 3, 8, 6, 5, 1, 7}, {9, 5, 3, 2, 7, 6, 4, 1, 8}}, 412}

For 4X4 matrices:

ClearAll[f2]
f2[samp_List] := Det[Partition[samp, 4]]

MaximizeOverPermutations[f2, 16, {1/100, 10}, 10^4] // AbsoluteTiming


{1.0341155, {{11, 4, 5, 15, 1, 9, 14, 10, 8, 16, 3, 7, 13, 6, 12, 2}, 40800.}}

and 5X5:

ClearAll[f3]
f3[samp_List] := Det[Partition[samp, 5]]

MaximizeOverPermutations[f3, 25, {1/100, 10}, 10^4] // AbsoluteTiming


{1.3268, {{1, 8, 17, 20, 19, 15, 13, 24, 11, 2, 22, 4, 12, 6, 21, 9, 25, 10, 5, 16, 18, 14, 3, 23, 7}, 6.83813*^6}}

• MOP is a wonderful function. Note that at its core it is simulated annealing and not a genetic method. (Presumably the original poster really wants for any working method, not necessarily one from the evolutionary algorithms family.) – Daniel Lichtblau Jan 31 '20 at 16:31
• This yields 6.838.13*  for n==5. But http://oeis.org/search?q=10%2C412%2C40800%2C&sort=&language=english&go=Search gives 6.839.492 as maximum for n==5. What is true? – Akku14 Jan 31 '20 at 16:43
• @Akku14, I think changing 10^4 in the 4th argument to a larger number we can get closer to the OEIS number. – kglr Jan 31 '20 at 16:49
• Very nice! Glad to see my code is finding usage. – Roman Feb 1 '20 at 13:06
• You can save some time by defining f[samp_List] := Det[Partition[samp, n] // N]: for some reason the floating-point algorithm is faster than the integer algorithm. – Roman Feb 1 '20 at 18:41

Maximize works for 2*2 matrix, but gets out of memory for 3*3 on my 16 GB computer.

max4 = Maximize[{Det@Partition[Array[a, 4], 2],
And @@ Thread[0 < Array[a, 4] < 5] &&
Unequal @@ Array[a, 4]},
Array[a, 4], Integers]

(* {10, {a[1] -> 3, a[2] -> 1, a[3] -> 2, a[4] -> 4}} *)

nmax9 = Maximize[{Det@Partition[Array[a, 9], 3],
And @@ Thread[0 < Array[a, 9] < 10] &&
Unequal @@ Array[a, 9]},
Array[a, 9], Integers]

(*   No more memory available.
Mathematica kernel has shut down.   *)


Edit

NMaximize with method "SimulatedAnnealing" gives you at least one of the 36 solutions within one minute.

nmax9 = NMaximize[{Det@Partition[Array[a, 9], 3],
And @@ Thread[0 < Array[a, 9] < 10] && Unequal @@ Array[a, 9] &&
Element[Array[a, 9], Integers]}, Array[a, 9],
MaxIterations -> 10000, Method -> "SimulatedAnnealing"]

(*   {412., {a[1] -> 1, a[2] -> 8, a[3] -> 4,
a[4] -> 5, a[5] -> 3, a[6] -> 9, a[7] -> 7, a[8] -> 6, a[9] -> 2}}   *)


Method "DifferentialEvolution" finds an other of the 36 solutions.

(*   {412., {a[1] -> 3, a[2] -> 8, a[3] -> 6, a[4] -> 5, a[5] -> 1,
a[6] -> 7, a[7] -> 9, a[8] -> 4, a[9] -> 2}}*)

• i don't get it: how can the kernel run out of memory if the total population is fixed in a genetic algorithm? from what i know the new generation take the place of the lower ranked members of the previous generation – Alucard Jan 31 '20 at 8:42
• I must confess, I don't know, what exactly a "genetic" algorithm is. I copied it from the OP. – Akku14 Jan 31 '20 at 8:46
• towardsdatascience.com/… i am a newbie too, i googled it because it sounds interesting – Alucard Jan 31 '20 at 8:50
• @Alucard There is a simple genetic algorithm code in the following answer. – A little mouse on the pampas Jan 31 '20 at 10:09
• @alucarD It was Maximize that ran out of memory, not NMaximize. The former uses exact methods, in this case either involving Lagrange multipliers and exact polynomial system solving, or else some form of cylindric decomposition. Both can be explosively consumptive of RAM. – Daniel Lichtblau Jan 31 '20 at 16:28

A "genetic algorithm" of this question:

(*QQ:2636051698*)
maxDet[n_, size_, iterations_] :=
Module[{fitness, choose, mutation, mutationInGroup, result},
fitness =
Compile[{{list, _Real, 1}},
Evaluate@Det@Quiet@Array[list[[n # + #2 - n]] &, {n, n}],
RuntimeAttributes -> Listable];
choose[group_] := group[[Ordering[Abs@fitness@group, -size]]];
mutation =
Compile[{{gene, _Real, 1}},
Module[{list = gene, changePos = RandomSample[Range[n^2], 3]},
list[[changePos]] = RandomSample[list[[changePos]]];
list], RuntimeAttributes -> Listable];
mutationInGroup[group_] := group~Join~(mutation /@ group);
result =
Nest[choose@mutationInGroup@# &,
Table[RandomSample[Range[n^2]], {i, 1, 10^5}], iterations][[1]];
{MatrixForm[Partition[result, n]], Round@fitness@result}];
maxDet[4, 100, 1000] // AbsoluteTiming

• Sorry, but this code failes for n >= 5 on my MMA Version 8.0. Maximia are { 1, 10, 412, 40800, 6839492, 1865999570, 762150368499} according to http://oeis.org/search?q=10%2C412%2C40800%2C&sort=&language=english&go=Search . – Akku14 Jan 31 '20 at 10:17
• @Akku14 This needs to wait for the others of the SE community to revise and speed up or open a new post tomorrow. – A little mouse on the pampas Jan 31 '20 at 10:36
• @Akku14 Someone has used genetic algorithm to answer this question. – A little mouse on the pampas Feb 1 '20 at 2:24