[R] adaptive optimization of mesh size

[baptiste Auguié]

DeaR list,


I'm running an external program that computes some electromagnetic  
response of a scattering body. The numerical scheme is based on a  
discretization with a characteristic mesh size "y". The smaller y is,  
the better the result (but obviously the computation will take longer).

A convergence study showed the error between the computed values and  
the exact solution of the problem to be a quadratic in y, with  
standard error increasing as y^3. I wrote the interface to the  
program in R, as it is much more user friendly and allows for post- 
processing analysis. Currently, it only runs with user-defined  
discretization parameter. I would like to implement an adaptive  
scheme [1] and provide the following improvements,

1) obtain an estimate of the error by fitting the result against a  
series of mesh sizes with the quadratic model, and extrapolate at y =  
0. (quite straight forward)

2) adapt "dynamically" the set of mesh sizes to fulfill a final  
accuracy condition, between a starting value (a rule-of thumb  
estimate is given by the problem values). The lower limit of y should  
also be constrained by the resources (again, an empirical rule  
dictates the computation time and memory usage).

I'm looking for advice on this second point (both on the technical  
aspect, and whether this is sound statistically):

- I can foresee that I should always start with a few y values before  
I can do any extrapolation, but how many of them? 3, 10? How could I  
know?

- once I have enough points (say, 10) to use the fitting procedure  
and get an estimate of the error, how should I decide the best  
location of the next y if the error is too important?

- in a practical implementation, I would use a while loop and append  
the successive values to a data.frame(y, value). However, this  
procedure will be run for different parameters (wavelengths,  
actually), so the set and number of y values may vary between one run  
and another. I think I'd be better off using a list with each new run  
having its own data.frame. Does this make sense?


Below are a few lines of code to illustrate the problem,


> program.result <- function(x, p){ # made up function that mimicks  
> the results of the real program
> 	
> 	y <- p[3]*x^2 + p[2]*x + p[1]
> 	y * (1 + rnorm(1, mean=0,  sd =  0.1 * y^3))
> }
>
>
> p0 <- c(0.1, 0.1, 2) # set of parameters
>
> ## user defined limits of the y parameter (log scale)
> limits <- c(0.1, 0.8)
> limits.log <- (10^limits)
> y.log <- seq(limits.log[1], limits.log[2], l=10)
>
> y <- log10(y.log)
>
> result <- sapply(y, function(x) program.result(x, p0)) # results of  
> the program
>
> #### fitting and extrapolation procedure ####
> library(gplots) # plot with CI
> plotCI(y,  result, y^3,  xlim=c(0, 1), ylim=c(0, 2)) # the data  
> with y^3 errors
>
> my.data <- data.frame(y = y,  value = result)
>
> fm <- lm(value ~ poly(y, degree=2, raw=TRUE), data = my.data ,   
> weights = 1/y^3)
> lines(y, predict(fm, data.frame(y=y)), col = 2)
>
> extrap <- summary(fm)$coefficients[1,] # intercept and error on it
> plotCI(0,extrap[1], 2 * extrap[2],  col = 2, add=T)
>
>
> ### my naive take on adaptive runs... ##
>
> objective <- 1e-3 # stop when the standard error of the  
> extrapolated value is smaller than this
>
> err <- extrap[2]
> my.color <- 3
> while (err > objective){
> 	
> 	new.value <- min(y)/2 # i don't know how to choose this optimally
> 	y <- c(new.value, y)
> 	new.result <- program.result(new.value, p0)
> 	result <- c(new.result, result)
> 	points(new.value, new.result, col= my.color)
> 	my.data <- data.frame(y = y,  value = result)
> 	fm <- lm(value ~ poly(y, degree=2, raw=TRUE), data = my.data ,   
> weights = 1/y^3)
> 	lines(y, predict(fm, data.frame(y=y)), col = my.color)
>
> 	extrap <- summary(fm)$coefficients[1,] # intercept and error on it
> 	err <- extrap[2]
> 	print(err)
> 	plotCI(0,extrap[1], 2 * err,  col = 2, add=T)
> 	my.color <- my.color + 1
> 	
> }
> err



Many thanks in advance for your comments,

baptiste


[1]: Yurkin et al., Convergence of the discrete dipole approximation.  
II. An extrapolation technique to increase the
accuracy. J. Opt. Soc. Am. A / Vol. 23, No. 10 / October 2006
_____________________________

Baptiste Auguié

Physics Department
University of Exeter
Stocker Road,
Exeter, Devon,
EX4 4QL, UK

Phone: +44 1392 264187

http://newton.ex.ac.uk/research/emag
http://projects.ex.ac.uk/atto