go.renorm {medea}R Documentation

Descent by renormalisation in space of measures under linear constraints.

Description

Minimization (maximization) of an objective function defined on the space of measures under linear constraints (currently not more than 3 constraints are supported). Examples of linear constaints on a measure: fixed total mass, fixed barycentre, etc. If the only constraint is the fixed mass, use the function go.steep of the library mefista instead.

Usage

med<-go.renorm(mu,obj.func,grad.func,
                 obj.func.name=NA,grad.func.name=NA,
                 minim=TRUE,value=NA,grad=NA,tol=0.01,
                 st.size=0.1*sum(mu), No.steps=100,space=NA,
                 details=FALSE,visualize=FALSE,smallest.st.size=0.000001,
                 draw.style=4,xgrid=NA,ygrid=NA,
                 tmp.file=NA,alpha=NA,base=0.5,
                 uni=FALSE,dim=NA,true.m=NA,nlevels=15,ncolors=0,...)

med<-go.renorm(med,obj.func.name=med$obj.func$name,
                  grad.func.name=med$grad.func$name,
                  minim=med$minim,value=med$value,
                  grad=med$grad,tol=med$param$tol,st.size=med$st.size,
                  No.steps=100,space=med$param$space,
                  details=FALSE,visualize=FALSE,smallest.st.size=0.000001,
                  draw.style=4,xgrid=med$param$xgrid,ygrid=med$param$ygrid,
                  tmp.file=med$param$file,alpha=med$param$alpha,
                  base=med$param$base,uni=med$param$uni,
                  true.m=med$true.m,nlevels=15,ncolors=0,...)

Arguments

mu A vector of the measure weights of all the points of the state space.
med An object of class medea - the result of a previous application of method go.renorm.
obj.func The objective function to minimize (maximize). It should have the form <obj.func.name>(mu,space,<other parameters>,...), where mu is a vector representing the measure, space is the state space as below. Part ... should be present if you wish to submit additional parameters to the plots. If the other parameters are present, they should be listed in ... part of the go.renorm() call.
grad.func The gradient function corresponding to F=obj.func, i.e. a function of the same form <grad.func.name>(mu,space,<other parameters>,...) as the obj.func above. It should return a vector grad(x) such that

lim_{t->+0} [F(mu+t*eta)-F(mu)]/t=int_{space} grad(mu,x) eta(dx)

for any measure eta. For discrete spaces mu=(mu_1,...,mu_n) and

grad(mu,x_i)=(d/dmu_i)F(mu).

The function should return Inf if not defined for some mu (see Example below).

obj.func.name The name of the object function. Unless you want it to appear in the summary differently, you need not set this parameter: the name will be retrieved automatically.
grad.func.name The name of the gradient function, respectively. Retrieved automatically.
minim Set this to FALSE for maximization of the objective function.
fval The initial value of obj.func, if known.
grad The initial value of grad.func, if known.
tol Tolerance: the gradient function is considered to be a constant on a set A if its variation on A is at most tol times the total variation of the gradient function. Equally, the objective function is considered to achieve the known true.m value (see below) if it is within tol error from it.
st.size The initial step size: half of the total variation of the increment measure.
No.steps Maximal number of steps to go.
space The state space of the measures. A vector or a matrix giving the coordinates of the points in each row of the latter. If not supplied it is defined as expand.grid(xgrid, ygrid) if the latter are available, otherwise it is considered to be 1:length(mu).
details If TRUE additional information on values of the objective function, the current step size and the variation of the gradient at each step will be printed. Otherwise, the output is a sequence of dots corresponding to steps, "+" and "-" corresponding to increase or decrease of the step size, respectively.
visualize If TRUE, the graph of the obtained measure and the gradient at each step will be plotted.
smallest.st.size The program is stopped when the value of st.size becomes less than this value.
draw.style Values: 0 - linear, 1 - image(), 2 - contour(), 3 - image() and contour superposed, 4 - persp(). Additional parameters to plot like nlevels - the number of levels, may be supplied in ... part of the call. If dimension dim of the space is not 2 this parameter is automatically set to 0 - the measure and the gradient will be shown on the plot as a linear sequence of values.
xgrid, ygrid In dimension 2 the state space can be defined as expand.grid(xgrid,ygrid) - the Cartesian product of the vectors. 2D plotting procedures will then be used by default.
tmp.file A character string (name of the file) to store the currently obtained measure.
base The coefficient the step size is multiplied by when the latter is decreased. Should be between 0 and 1.
uni logical: see section Details below.
dim Dimension of the state space: length(space[1,]).
true.m The "true minimal" (maximal) value of obj.func if known.
nlevels The number of level lines to draw in 2D drawing procedures.
ncolors The number of colors to use in 2D drawing procedures. You may run out of color allocated memory if this parameter is too large. If 0, the default color palette is used.
constr The function that calculates constraints. It should have the form <constr.func.name>(x,<other parameters>,...), where x is a point in space. This function returns the vector of dimension equal to the number of contraints. Some typical examples of constraints are in the library, for example, bary - constraint on the total mass and the barycentre, bary.sq - constraint on the total mass, the barycentre and the second moment
cval The vector of dimension equal to the number of constraints and providing values of constraints. This vector should match the values of constraints on the initial measure mu. If they do not match ERROR is returned. If only mu is specified, cval are calculated automatically. If only cval are given and mu is NA, the program tries to find an initial measure to satisfy the imposed constraints and fails if they are contradictory.
d.mat The design matrix that, if not specified, is calculated by applying constr function to all elements of space.
c.con The matrix obtained by inverting constraints. It specifies masses of points (their number k if equal to the number of constraints) such that the corresponding atomic measure satisfies the imposed constraints and is non-negative. It also creates an auxiliary matrix by applying expand.grid to the space repeated k times. As the last row it returns the number of the element of this auxiliary matrix. Another auxiliary programme by.div calculates the positions of the corresponding points.

Details

The following algorithm is implemented.

Incremental measures of total variation 2*st.size.cur chosen so that to minimize the integral of the gradient are added to the current measure mu while obj.func diminishes. Then a new gradient is calculated. If m increments were made at the previous step then the current step size st.size.cur is increased m times. If no increments were possible at the current step size then the latter is reduced by base<1. When constructing the increment measure, the algorithm adds positive masses to maintain constraints and minimise the gradient function. The obtained measure is then renormalised to satisfy the imposed constraints.

The program runs until at least one of the following conditions is satisfied: the maximal count of steps is achieved; the step size becomes too small; the current measure satisfies the necessary optimality condition: the gradient function is minimal (maximal) and constant on the support of the optimal measure; the objective function equals the known true.m value within error tol. The latter is important, e.g., in the case of equation solving: F(mu)=C is equivalent to minimization of (F(mu)-C)^2 with true.m=0.

In computationally extensive cases it is advisable to set tmp.file to a name of the file to store the currently obtained measure. You then will be able to retrieve it with mu<-dget("<name.of.tmp.file>") and reuse it if the program crushes or needs to be interrupted.

To avoid warnings during visualisation it is advisable to set par(err=-1).

Value

Returns an object of class medea with the following named components:

mes the current measure.
obj.func the list with 2 components:
obj.func$name the name of the objective function,
obj.func$code the objective function itself.
grad.func the list with 2 components:
grad.func$name the name of the gradient function,
grad.func$code the gradient function itself.
minim TRUE for minimization problem, FALSE otherwise.
fval the value of the objective function.
grad the value of the gradient function.
var.grad the total variation of the gradient on the support of the measure.
optim TRUE if optimum is attained, FALSE otherwise.
st.size the last value of the step size (the total variation of the increment measure).
atoms the formatted character string giving the atoms of the measure. Use cat(<medea.object>$atoms) to see them.
No.atoms the number of the measure's atoms.
tot.steps the total number of steps made during the last call.
true.m the value of the true minimum (maximum), if known.
d.mat the value of constraints applied to the space.
c.con the result of inversion of the constraints matrix.
param the list with the following components:
param$tol the tolerance,
param$init.st.size the initial step size,
param$smallest.st.size the smallest step size,
param$uni the logical value of parameter uni,
param$space the matrix of the state space,
param$xgrid vector xgrid,
param$ygrid vector ygrid,
param$dim the dimension of the state space,
param$constr the function that determines constraints.
param$cval the values of the constraints.
param$base the parameter of the step size decrease,
param$file the name of the file used to store the obtained measure.
other the other parameters given in the ... part of the call.

Author(s)

Ilya Molchanov, E-mail: ilya@stat.unibe.ch
Sergei Zuyev, E-mail: sergei@stams.strath.ac.uk

References

[1] I. Molchanov and S. Zuyev. "Steepest descent algorithms in space of measures". Statistics and Computing, 12, 115-123, 2002.
http://www.cx.unibe.ch/~ilya
http://www.stams.strath.ac.uk/~sergei

See Also

summary.medea, plot.medea and go.steep of the library mefista.

Examples

# Optimal design of observation points x_i in linear model
# y(x_i)=\sum_i f(x_i)*beta
# f is a vector function (row) and beta is a column-vector of unknown coeffs.
# Example from Atkinson & Donev "Optimum Experimental Design",
# Oxford 1992 (p.119): cubic regression through the origin.
# We are looking for an optimal design having the total mass 1
# and the barycenter 0.7
#
library(medea)
info.matrix <- function(mu,space,func,...)
{ # Information matrix
        return(t(func(space)) %*% diag(mu) %*% func(space))
}

detIM.func<-function(mu,space,...)
  { # The object function: determinant of the info matrix
       return(det(info.matrix(mu,space,...))) 
  }

detIM.grad <- function(mu,space,func,...)
  { # The gradient function
    inf.mat <- info.matrix(mu,space,func,...)
    if(det(inf.mat)>1e-10)
        res <- diag(func(space) %*% solve(inf.mat) %*% t(func(space)))
    else
        res <- rep(Inf,length=length(mu)) # exclude singular matrix case
    return( res)
  }

 x<-seq(0,1,by=0.05) # State space
 f1<-function(x) {return(x)}
 f2<-function(x) {return(x^2)}
 f3<-function(x) {return(x^3)}
 f<-function(x) {return(cbind(f1(x),f2(x),f3(x)))}

 # And here we go.renorm:
 ## Not run: 
res<-go.renorm(mu=NA,detIM.func,detIM.grad,minim=FALSE,tol=0.01,space=x,
      No.steps=10,details=TRUE,visualize=TRUE,func=f,constr=bary,
      cval=c(1,0.7))
 res
 res<-go.renorm(res,visualize=TRUE)
 # ATTENTION: the next command performs a few hundred iterations
 # to arive to the optimal solution:
 # while(!res$opt) res<-go.renorm(res)
 summary(res)
 plot(res)
## End(Not run)

[Package medea version 1.01 Index]