# 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.
Commented Jan 31, 2020 at 3:54
– kglr
Commented Jan 31, 2020 at 8:33

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/… Commented Feb 1, 2020 at 7:05
• The script now is general for any dimension dim Commented Feb 1, 2020 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.) Commented Jan 31, 2020 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? Commented Jan 31, 2020 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
Commented Jan 31, 2020 at 16:49
• Very nice! Glad to see my code is finding usage. Commented Feb 1, 2020 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. Commented Feb 1, 2020 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 Commented Jan 31, 2020 at 8:42
• I must confess, I don't know, what exactly a "genetic" algorithm is. I copied it from the OP. Commented Jan 31, 2020 at 8:46
• towardsdatascience.com/… i am a newbie too, i googled it because it sounds interesting Commented Jan 31, 2020 at 8:50
• @Alucard There is a simple genetic algorithm code in the following answer. Commented Jan 31, 2020 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. Commented Jan 31, 2020 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 . Commented Jan 31, 2020 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. Commented Jan 31, 2020 at 10:36
• @Akku14 Someone has used genetic algorithm to answer this question. Commented Feb 1, 2020 at 2:24