Home › Blogs › Conor's blog

Testing optimisation algorithms: Part III - Global optimisation in R

Sun, 09/18/2011 - 20:12 — Conor

In part I of this series of posts, I presented the maths behind posing circle packing as a high-dimensional, non-linear, constrained optimisation problem. I am using this problem as an example to try out various optimisation algorithms. In part II I demonstrated an implementation of this problem in R, together with an R function for visualising candidate solutions.

In this post, I will describe some findings from attempts to optimise the circle packing objective function with R's default, general-purpose optimisation package: optim and two other packages for performing global optimisation using genetic algorithms with derivative-based intermediate searches.

I have updated the code from the last post slightly (added two more small functions). You can download the updated version here: CirclePacking.R. I've also written a second script demonstrating the use of the functions in CirclePacking.R in combination with the 3 optimisation packages mentioned above: TestingOptimisation.R (duplicated at the bottom of this post). Running this code generates plots like the one below, visualising solutions, together with some statistics about running time distributions and packing density distributions for ensembles of solutions.

L-BFGS-BN4_200
64 solutions for packing 200 circles into a square with the optim package.

The optim package includes an implementation of the L-BFGS-B algorithm, which is an efficient, derivative-based local optimiser allowing lower and upper bounds to be placed on parameters.

The rgenoud package uses genetic algorithms and derivatives to solve global optimisation problems. I used it successfully to fit a mixture model to pixel intensity distributions as part of the segmentation step in the image analysis tool: Colonyzer. It also has some interesting functionality for distributing computing load across multiple CPUs, by using a snow cluster which I wanted to investigate.

The DEoptim package also uses genetic algorithms and derivatives to solve global optimisation problems. I have previously found that this algorithm was faster than rgenoud when automatically fitting a logistic model to yeast growth curves in Quantitative Fitness Analysis.

Results

My findings are problem-specific, and depend on the number of circles allowed (which you can vary by modifying the code below) but so far, they are as follows:

1) For all N I've tried (up to 200), the L-BFGS-B algorithm was fastest. It achieved the same packing densities, but in shorter times.

2) Using the rgenoud cluster functionality did not speed up the genoud algorithm, or improve the packing densities of its solutions.

3) DEoptim crashes with mysterious errors for N>200 (more than 600 parameters) on a 64-bit machine with 12 Gb RAM. However, this is probably many more parameters than are required for most practical purposes.

4) rgenoud was much faster than DEoptim in this case (achieved the same or higher packing densities more quickly).

It's difficult to draw solid conclusions from this analysis, as I now suspect that for the circle packing problem, there are many local minima with largely similar packing densities and so it may be difficult or impossible to distinguish between local optima and a global optimum (if it exists). This would explain result 1) and possibly result 2). Result 1) was particularly surprising since I expected this high-dimensional problem to be very vulnerable to getting stuck in local minima, and so I expected global optimisers to be more efficient. As mentioned above, result 4) is reversed for at least one other optimisation problem.

I'm happy with the solutions (see image above), they have very pleasing organic patterns.


  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
 
library(DEoptim)
library(rgenoud)
library(snow)
source("CirclePacking.R")

# Number of circles and dimensions of bounding rectangle 
N=100; W=10; H=10
root="N64_100"
arrnum=64
arrdim=ceiling(sqrt(arrnum))

# Create objective function
obj=createObj(N,W,H)

# Lowest acceptable values for dimensions (x,y) is zero
# Minimum radius (to prevent circles being optimised out of existance)
low=rep(c(0,0,0.01*min(W,H)),N)
# x can go as high as W, y as high as H
# Radius of biggest possible circle contained in box is the smaller of W/2 and H/2
up=rep(c(W,H,min(W,H)/2),N)

# Generate initial population of guesses
NumPart=10*3*N
pop=matrix(0,nrow=NumPart,ncol=3*N)
for(i in 1:NumPart) pop[i,]=genGuess(N,W,H)

# L-BFGS-B test
results=c()
times=c()
pdf(paste("L-BFGS-B",root,".pdf",sep=""))
op<-par(mfrow=c(arrdim,arrdim))
for(i in 1:arrnum){
	print(i)
	z=genGuess(N,W,H); #plotCircles(z,W,H,numbers=TRUE) # Starting guess
	tim=system.time({
		out=optim(par=z,fn=obj,method="L-BFGS-B",lower=low,upper=up,control=list(maxit=300)) # Optimise
		z=as.numeric(out$par);  results=c(results,plotCircles(z,W,H,numbers=FALSE))# Result
	})
	times=c(times,as.numeric(tim[1]))
}
par(op)
op<-par(plt=c(0.125,0.95,0.13,0.9))
plot(times,results,pch=16)
hist(times)
hist(results)
par(op)
dev.off()

# parallel genoud Differential optimisation test
doms=matrix(c(low,up),nrow=3*N,ncol=2,byrow=FALSE) # Format bounds
results=c()
times=c()
cl <- makeCluster(12, type = "SOCK") # Initiate cluster
# Unfortunately, snow clusters don't seem to respect local R variables, so we pass them as global variables to each node, like this...
clusterExport(cl,list("N","W","H"))
pdf(paste("Parallelrgenoud",root,".pdf",sep=""))
op<-par(mfrow=c(arrdim,arrdim))
for(i in 1:arrnum){
	print(i)
	tim=system.time({
		z=genGuess(N,W,H); #plotCircles(z,W,H,numbers=TRUE) # Starting guess
		out=genoud(obj,3*N,starting.values=z,boundary.enforcement=2,Domains=doms, cluster=cl, max.generations=10,print.level=0,control=list(maxit=300))# Optimise
		z=out$par; results=c(results,plotCircles(z,W,H,numbers=FALSE))# Result
	})
	times=c(times,as.numeric(tim[1]))
}
stopCluster(cl)
par(op)
op<-par(plt=c(0.125,0.95,0.13,0.9))
plot(times,results,pch=16)
hist(times)
hist(results)
par(op)
dev.off()

# genoud Differential optimisation test
doms=matrix(c(low,up),nrow=3*N,ncol=2,byrow=FALSE) # Format bounds
results=c()
times=c()
pdf(paste("rgenoud",root,".pdf",sep=""))
op<-par(mfrow=c(arrdim,arrdim))
for(i in 1:arrnum){
	print(i)
	tim=system.time({
		z=genGuess(N,W,H); #plotCircles(z,W,H,numbers=TRUE) # Starting guess
		out=genoud(obj,3*N,starting.values=z,boundary.enforcement=2,Domains=doms, max.generations=10,print.level=0,control=list(maxit=300))# Optimise
		z=out$par; results=c(results,plotCircles(z,W,H,numbers=FALSE))# Result
	})
	times=c(times,as.numeric(tim[1]))
}
par(op)
op<-par(plt=c(0.125,0.95,0.13,0.9))
plot(times,results,pch=16)
hist(times)
hist(results)
par(op)
dev.off()

# DEoptim Differential optimisation test
results=c()
times=c()
pdf(paste("DEoptim",root,".pdf",sep=""))
op<-par(mfrow=c(arrdim,arrdim))
for(i in 1:arrnum){
	print(i)
	pop=matrix(0,nrow=NumPart,ncol=3*N)
	for(i in 1:NumPart) pop[i,]=genGuess(N,W,H); # Starting guesses
	tim=system.time({
		out=DEoptim(obj, low, up, DEoptim.control(trace=0,itermax=2000,NP=NumPart,initialpop=pop)) # Optimise
		z=as.numeric(out$optim$bestmem);  results=c(results,plotCircles(z,W,H,numbers=FALSE))# Result
	})
	times=c(times,as.numeric(tim[1]))
}
par(op)
op<-par(plt=c(0.125,0.95,0.13,0.9))
plot(times,results,pch=16)
hist(times)
hist(results)
par(op)
dev.off()