How to use genetic algorithm to calculate the maximum value of this matrix
Follows a classical Genetic Algorithm with the normal functionalities as
- Population initialization
- Fitness Evaluation
- Crossover
- Mutation
- 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.
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}}
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}}*)