Simulated  Annealing

Integer Partitions (Subset Sum Problem)

Setup

Here is the same list of 50 integers as in the last class:

integers={26,52,44,44,76,30,46,34,27,90,97,91,43,14,36,86,93,90,11,14,62,71,65,47,72,13,79,69,95,10,50,37,69,27,13,27,78,50,53,12,94,95,81,22,36,74,54,28,18,64};
n=Length[integers]

As before, a state is a list of 50 zero/one values, indicating which integers are selected for the partition.

currState=RandomInteger[{0,1},n]

We can add up the selected integers like this:

Sum[currState[[i]]*integers[[i]],{i,1,50}]

We also need an objective function to maximize:

f[state_,integers_]:=-Abs[target - Sum[state[[i]]*integers[[i]],{i,1,50}]]

A more sophisticated transition function

In the last class, our proposal transition simply involved flipping one integer. What happens if we instead flip two integers simultaneously?

getProposal[state_]:=Module[
{i,propState=state,n=Length[state]},

i=RandomInteger[{1,50}];
(*Print["flipping one number: index ",i];*)
propState[[i]]=1-propState[[i]];

propState
]

Simulated Annealing

Here is a module that runs our simulated annealing algorithm for finding a subset of the integers with the target sum.

findSubset[integers_,target_]:=Module[
(* local variables *)
{n=Length[integers],
β=1,βinc=1.001,
fVals={} ,(* we will store the values of f in this list *)
maxIter=1000,
propState,currState,
ρ,Δf
},

(* generate an initial state *)
currState=RandomInteger[{0,1},n];

(* main loop *)
Do[
(* select a proposed next state *)
propState=getProposal[currState];

(* select the next β in the annealing schedule *)
β*=βinc;

(* compute the acceptance ratio *) 
Δf=f[propState,integers]-f[currState,integers];
ρ=Exp[β*Δf];

(* decide whether to accept the proposed transition *) 
If[RandomReal[]<ρ,
currState=propState
];

(* if a partition has been found, then quit *)
AppendTo[fVals,f[currState,integers]];
If[f[currState,integers]==0,Print["Success!"];Return[] ]
,maxIter];

(* return values *) 
{β,fVals,currState}
]

{β,fVals,finalState}=findSubset[integers,1000]

Check our results

Make a list plot of the values of the function f throughout the simulated annealing process:

ListPlot[fList,Joined->True]

Here is one way to see the integers that add up to 1000:

solution={};
Do[
If[finalState[[i]]==1,AppendTo[solution,integers[[i]]]],
{i,1,50}]
solution
Total[solution]

Magic  Squares

Now let's use simulated annealing to find magic squares

Setup

We can generate a random arrangement of the integers 1, 2, ..., 9 in a 3×3 grid like this:

n=3;
square=RandomSample[Range[(n^2)]]

square=Partition[RandomSample[Range[(n^2)]],n]

Grid[square,Frame->All]

Here is the formula for the n×n magic sum:

magicSum[n_]:=((n((n^2)+1))/2)

Let's also write a function that adds up the row, column, and diagonal sums of the square:

getSums[square_]:=Module[
{sums,n=Length[square]},

(* column sums *)
sums=Total[square];

(* row sums *)
sums=Join[sums,Total[square,{2}]];

(* diagonal sums *)
AppendTo[sums,Sum[square[[i,i]],{i,1,n}]];
AppendTo[sums,Sum[square[[n+1-i,i]],{i,1,n}]];

(* return the list of sums *) 
sums
]

Objective function

What objective function should we define? We want a function that quantifies how far we are from a magic square. This function should be maximized when we find a magic square.

Transitions

What transition function should we use?

Simulated  Annealing

Now use simulated annealing to find a magic square. The following module might help.

findMagicSquare[n_,f_]:=Module[
{β=1,βinc=1.001,
propState,currState,
target=magicSum[n],
ρ,Δf,
fVals={}, (* we will store the values of f in this list *)
maxIter=10
},

(* generate a random n×n square *) 
currState=Partition[RandomSample[Range[(n^2)]],n];

(* main loop *)
Do[
(* select a proposed next state *)
propState=getProposal[currState];

(* select the next β in the annealing schedule *)
β*=βinc;

(* compute the acceptance ratio *) 
Δf=f[propState,target]-f[currState,target];
ρ=Exp[β*Δf];

(* decide whether to accept the proposed transition *) 
If[RandomReal[]<ρ,
currState=propState
];

(* if a partition has been found, then quit *)
AppendTo[fVals,f[currState,target]];
If[f[currState,target]==0,Print["Success!"];Break[] ]
,maxIter];

(* return values *) 
{β,fVals,currState}
]

Questions

Can you find a 4×4 magic squares? How about 5×5 magic squares? Larger?

On average, how many iterations are required before a magic square is found?

What annealing schedule works best?

Which function to minimize works best?

Can you think of a better proposal transition function?