Archive for simulated annealing

stochastic magnetic bits, simulated annealing and Gibbs sampling

Posted in Statistics with tags , , , , , , , , , on October 17, 2019 by xi'an

A paper by Borders et al. in the 19 September issue of Nature offers an interesting mix of computing and electronics and optimisation. With two preparatory tribunes! One [rather overdone] on Feynman’s quest. As a possible alternative to quantum computers for creating probabilistic bits. And making machine learning (as an optimisation program) more efficient. And another one explaining more clearly what is in the paper. As well as the practical advantage of the approach over quantum computing. As for the paper itself, the part I understood about factorising an integer F via minimising the squared difference between a product of two integers and F and using simulated annealing sounded rather easy, while the part I did not about constructing a semi-conductor implementing this stochastic search sounded too technical (especially in the métro during rush hour). Even after checking the on-line supplementary material. Interestingly, the paper claims for higher efficiency thanks to asynchronicity than a regular Gibbs simulation of Boltzman machines, quoting Roberts and Sahu (1997) without further explanation and possibly out of context (as the latter is not concerned with optimisation).

ABC-SAEM

Posted in Books, Statistics, University life with tags , , , , , , , , , , , , , , , , , , on October 8, 2019 by xi'an

In connection with the recent PhD thesis defence of Juliette Chevallier, in which I took a somewhat virtual part for being physically in Warwick, I read a paper she wrote with Stéphanie Allassonnière on stochastic approximation versions of the EM algorithm. Computing the MAP estimator can be done via some adapted for simulated annealing versions of EM, possibly using MCMC as for instance in the Monolix software and its MCMC-SAEM algorithm. Where SA stands sometimes for stochastic approximation and sometimes for simulated annealing, originally developed by Gilles Celeux and Jean Diebolt, then reframed by Marc Lavielle and Eric Moulines [friends and coauthors]. With an MCMC step because the simulation of the latent variables involves an untractable normalising constant. (Contrary to this paper, Umberto Picchini and Adeline Samson proposed in 2015 a genuine ABC version of this approach, paper that I thought I missed—although I now remember discussing it with Adeline at JSM in Seattle—, ABC is used as a substitute for the conditional distribution of the latent variables given data and parameter. To be used as a substitute for the Q step of the (SA)EM algorithm. One more approximation step and one more simulation step and we would reach a form of ABC-Gibbs!) In this version, there are very few assumptions made on the approximation sequence, except that it converges with the iteration index to the true distribution (for a fixed observed sample) if convergence of ABC-SAEM is to happen. The paper takes as an illustrative sequence a collection of tempered versions of the true conditionals, but this is quite formal as I cannot fathom a feasible simulation from the tempered version and not from the untempered one. It is thus much more a version of tempered SAEM than truly connected with ABC (although a genuine ABC-EM version could be envisioned).

Le Monde puzzle [#1092]

Posted in Statistics with tags , , , , , , , on April 18, 2019 by xi'an

A Latin square Le Monde mathematical puzzle that I found rather dreary:

A hidden 3×3 board contains all numbers from 1 to 9. Anselm wants to guess the board and makes two proposals. Berenicke tells him how many entries are in the right rows and colums for each proposal, along with the information that no entry is at the right location. Anselm deduces the right board.

Which I solved by brute force and not even simulated annealing, first defining a target

ordoku1=ordoku2=matrix(1,9,2)
ordoku1[,1]=c(1,1,1,2,2,2,3,3,3)
ordoku1[,2]=rep(1:3,3)
ordoku2[,1]=c(3,2,3,1,2,3,2,1,1)
ordoku2[,2]=c(2,2,3,2,3,1,1,3,1)
fitz=function(ordo){
  (sum(ordo[c(1,4,7),2]==1)==1)+(sum(ordo[c(2,5,8),2]==2)==1)+
  (sum(ordo[c(3,6,9),2]==3)==0)+(sum(ordo[c(1,2,3),1]==1)==1)+
  (sum(ordo[c(4,5,6),1]==2)==1)+(sum(ordo[c(7,8,9),1]==3)==2)+
  (sum(ordo[c(6,7,9),2]==1)==2)+(sum(ordo[c(1,2,4),2]==2)==1)+
  (sum(ordo[c(3,5,8),2]==3)==2)+(sum(ordo[c(4,8,9),1]==1)==1)+
  (sum(ordo[c(7,2,5),1]==2)==1)+(sum(ordo[c(1,3,6),1]==3)==0)+
  (!(0%in%apply((ordo-ordoku1)^2,1,sum)))+(!(0%in%apply((ordo-ordoku2)^2,1,sum)))
}

on a 9×9 board entry reproducing all items of information given by Berenicke. If all constraints are met, the function returns 14. And then searched for a solution at random:

temp=1
randw=function(ordo){
  for (t in 1:1e6){
    chlg=sample(1:9,2)
    temp=ordo[chlg[1],]
    ordo[chlg[1],]=ordo[chlg[2],]
    ordo[chlg[2],]=temp
    if (fitz(ordo)==14){
      print(ordo);break()}}}

which produces the correct board

4 3 5
6 7 1
9 2 8

Le Monde puzzle [#1088]

Posted in Books, Kids, R with tags , , , , , , , , on March 29, 2019 by xi'an

A board (Ising!) Le Monde mathematical puzzle in the optimisation mode, again:

On a 7×7 board, what is the maximal number of locations that one can occupy when imposing at least two empty neighbours ?

Which I tried to solve by brute force and simulated annealing (what else?!), first defining a target

targ=function(tabz){
  sum(tabz[-c(1,9),-c(1,9)]-1.2*(tabz[-c(1,9),-c(1,9)]*tabz[-c(8,9),-c(1,9)]
      +tabz[-c(1,9),-c(1,9)]*tabz[-c(1,2),-c(1,9)]
      +tabz[-c(1,9),-c(1,9)]*tabz[-c(1,9),-c(8,9)]
      +tabz[-c(1,9),-c(1,9)]*tabz[-c(1,9),-c(1,2)]>2))}

on a 9×9 board where I penalise prohibited configuration by a factor 1.2 (a wee bit more than empty nodes). The perimeter of the 9×9 board is filled with ones and never actualised. (In the above convoluted products, the goal is to count how many neighbours of the entries equal to one are also equal to one. More than 2 is penalised.) The simulated annealing move is then updating the 9×9 grid gridz:

temp=1
maxarg=curarg=targ(gridz)
for (t in 1:1e3){
  for (v in 1:1e4){
    i=sample(2:8,1);j=sample(2:8,1)
    newgrid=gridz;newgrid[i,j]=1-gridz[i,j]
    newarg=targ(newgrid)
    if (log(runif(1))<temp*(newarg-curarg)){
      gridz=newgrid;curarg=newarg}}
temp=temp+.01}

and calls to the procedure always return 28 entries as the optimum, as in

     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    1    0    1    0    1    0    1
[2,]    0    1    1    0    1    1    0
[3,]    1    1    0    1    0    1    1
[4,]    0    0    1    0    1    0    0
[5,]    1    1    0    1    0    1    1
[6,]    0    1    1    0    1    1    0
[7,]    1    0    1    0    1    0    1

As it happens, I had misread the wording of the original puzzle, which considered a dynamic placement of the units on the board, one at a time with two free neighbours imposed.

Le Monde puzzle [#1087]

Posted in Books, Kids, R, Statistics with tags , , , , , on February 25, 2019 by xi'an

A board-like Le Monde mathematical puzzle in the digit category:

Given a (k,m) binary matrix, what is the maximum number S of entries with only one neighbour equal to one? Solve for k=m=2,…,13, and k=6,m=8.

For instance, for k=m=2, the matrix

\begin{matrix} 0 &0\\ 1 &1\\ \end{matrix}

is producing the maximal number 4. I first attempted a brute force random filling of these matrices with only a few steps of explorations and got the numbers 4,8,16,34,44,57,… for the first cases. Since I was convinced that the square k² of a number k previously exhibited to reach its maximum as S=k² was again perfect in this regard, I then tried another approach based on Gibbs sampling and annealing (what else?):

gibzit=function(k,m,A=1e2,N=1e2){
  temp=1 #temperature
  board=sole=matrix(sample(c(0,1),(k+2)*(m+2),rep=TRUE),k+2,m+2)
  board[1,]=board[k+2,]=board[,1]=board[,m+2]=0 #boundaries
  maxol=counter(board,k,m) #how many one-neighbours?
  for (t in 1:A){#annealing
    for (r in 1:N){#basic gibbs steps
      for (i in 2:(k+1))
        for (j in 2:(m+1)){
          prop=board
          prop[i,j]=1-board[i,j]
          u=runif(1)
          if (log(u/(1-u))<temp*(counter(prop,k,m)-val)){ 
             board=prop;val=counter(prop,k,m) 
             if (val>maxol){
               maxol=val;sole=board}}
      }}
    temp=temp*2}
  return(sole[-c(1,k+2),-c(1,m+2)])}

which leads systematically to the optimal solution, namely a perfect square k² when k is even and a perfect but one k²-1 when k is odd. When k=6, m=8, all entries can afford one neighbour exactly,

> gibzbbgiz(6,8)
[1] 48
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    0    0    1    1    0    0    1
[2,]    1    0    0    0    0    0    0    1
[3,]    0    0    1    0    0    1    0    0
[4,]    0    0    1    0    0    1    0    0
[5,]    1    0    0    0    0    0    0    1
[6,]    1    0    0    1    1    0    0    1

but this does not seem feasible when k=6, m=7, which only achieves 40 entries with one single neighbour.

Le Monde puzzle [#1075]

Posted in Books, Kids, R with tags , , , , , , , on November 22, 2018 by xi'an

A Le Monde mathematical puzzle from after the competition:

A sequence of five integers can only be modified by subtracting an integer N from two neighbours of an entry and adding 2N to the entry.  Given the configuration below, what is the minimal number of steps to reach non-negative entries everywhere? Is this feasible for any configuration?

As I quickly found a solution by hand in four steps, but missed the mathematical principle behind!, I was not very enthusiastic in trying a simulated annealing version by selecting the place to change inversely proportional to its value, but I eventually tried and also obtained the same solution:

      [,1] [,2] [,3] [,4] [,5]
   -3    1    1    1    1
    1   -1    1    1   -1
    0    1    0    1   -1
   -1    1    0    0    1
    1    0    0    0    0

But (update!) Jean-Louis Fouley came up with one step less!

      [,1] [,2] [,3] [,4] [,5]
   -3    1    1    1    1
    3   -2    1    1   -2
    2    0    0    1   -2
    1    0    0    0    0

The second part of the question is more interesting, but again without a clear mathematical lead, I could only attempt a large number of configurations and check whether all admitted “solutions”. So far none failed.

Le Monde puzzle [#1071]

Posted in Books, Kids with tags , , , , , , , , , , , on October 18, 2018 by xi'an

A “he said she said” Le Monde mathematical puzzle sixth competition problem that reminded me of the “Singapore birthday problem” (nothing to do with the original birthday problem!):

Arwen and Brandwein are privately and respectively given the day and month of Caradoc’s birthday [in the Gregorian calendar] with the information that the month number is at least the day number. Arwen starts by stating she knows Brandwein cannot deduce the birthday, followed by Brandwein who says the same about Arwen. If this “she says he says” goes on for the largest possible number of steps before Arwen says she knows, when is Caradoc’s birthday? Arwen and Brandwein are later given two numbers corresponding to Deirdre’s birthday with no indication of which one is the day and which one is the month. They know both numbers end up with the same digit and that the month number is strictly less than the day number. Arwen states she does not know the date and she knows Brandwein cannot know either. Then Brandwein says he indeed does not the date but he knows whether he got the day or the month. This prompts Arwen to conclude she knows, then Brandwein to do the same. When is Deirdre’s birthday?

Since this was a fairly easy puzzle (and since I had spent too much time debugging the previous R code!), I did not try coding this one but instead drew the possibilities and remove the impossibilities on a blackboard. The first question is quite simple actually since the day numbers stand between 1 and 12 and that each “I cannot know” excludes one of the remaining endpoints, removing first excludes 1 from both lists, then 12, then 2, then …. 8, ending up with 7. And 07/07 as Caradoc’s birthday. The second case sees 13,…,20,23,…,30 eliminated from Arwen’s numbers, then 3,…,10 as well, which eliminates the same numbers from Brandwein’s possibilities. That he knows whether it is a month or a day leaves only 1,2,21,22,31 as possible numbers. Then Arwen’s certainty reduces her numbers to be 2, 21, 22, or 31, and since Brandwein is also sure, the only possible cases are (2,22) and (22,2). Meaning Deirdre’s birthday is on 22/02. I dunno if this symmetry was to be expected! (And I cannot fathom why this puzzle is awarded so many points, when compared with the others.)