go.renorm {medea} | R Documentation |
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.
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,...)
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 |
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. |
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)
.
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. |
Ilya Molchanov, E-mail: ilya@stat.unibe.ch
Sergei Zuyev, E-mail: sergei@stams.strath.ac.uk
[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
summary.medea
,
plot.medea
and
go.steep
of the library mefista
.
# 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)