Nothing
R/Rtnmin-package.R
In Rtnmin: Truncated Newton Function Minimization with Bounds Constraints
Defines functions
initpc
lin1
modlnp
msolve
ndia3
step1
stpmax
ztime
cnvtst <- function (alpha, pnorm, xnorm,
dif, ftest, gnorm, gtp, f, flast, g, ipivot, accrcy) {
## ---------------------------------------------------------
## test for convergence
## ---------------------------------------------------------
## set up
## ---------------------------------------------------------
conv <- 0;
eps <- .Machine$double.eps
toleps <- sqrt(accrcy) + sqrt(eps);
rtleps <- accrcy + eps;
imax <- 0;
ltest <- (flast - f <= -0.5*gtp);
## ---------------------------------------------------------
## if anti-zigzag test satisfied, test multipliers;
## if appropriate, modify active set
## ---------------------------------------------------------
if ( ! ltest) {
ind <- which( (ipivot != 0) & (ipivot != 2))
if ( length(ind) > 0 ) { ## how to ensure ind not empty
t <- -sum(ipivot[ind]*g[ind])
cmax <- min(t) # ?? why [cmax, imax] = min(t) ??
imax <- cmax
if (cmax >= 0) { imax <- 0 }
}
}
if (imax != 0) {
ipivot[ ind[imax] ] <- 0;
flast <- f;
} else {
conv = ( ( (alpha*pnorm < toleps*(1 + xnorm) )
&& (abs(dif) < rtleps*ftest)
&& (gnorm < accrcy^(1/3)*ftest) ) ||
( gnorm < .01*sqrt(accrcy)*ftest) )
}
result <- list(conv=conv, flast1=flast, ipivot1=ipivot)
}
crash <- function (x, low, up) {
##---------------------------------------------------------
## this initializes the constraint information, and
## ensures that the initial point satisfies
## low <= x <= up.
## the constraints are checked for consistency.
##---------------------------------------------------------
ierror <- 0
if (any(low > up)) { ierror = - max(which(low > up)) }
# above is check on error in bounds specification
xnew <- pmax (low, x) # force params into bounds
xnew <- pmin (up, xnew) # No diagnostic! ??
# we output revised parameters
n <- length(x)
ind <- which(low == x)
ipivot <- rep(0, n) ## zeros(size(x)) in Matlab
if (length(ind) > 0) { ipivot[ind] <- -1 }
ind <- which(x == up)
if (length(ind) > 0) { ipivot[ind] <- 1 }
ind <- which(low == up)
if (length(ind) > 0) { ipivot[ind] <- 2 }
list(ipivot=ipivot, ierror=ierror, xnew=xnew)
}
gtims <- function (v, x, g, accrcy, xnorm, sfun, ...) {
##---------------------------------------------------------
## compute the product of the Hessian times the vector v;
## store result in the vector gv
## (finite-difference version)
##---------------------------------------------------------
## cat("v, x, g:")
## print(v)
## print(x)
## print(g)
delta <- sqrt(accrcy)*(1 + xnorm)/sqrt(sum(v^2))
hg <- x + delta*v
tresult<-sfun(hg, ...)
## cat("tresult: ")
## print(tresult)
gv <- attr(tresult, "gradient")
gv <- (gv - g)/delta
gv
}
initpc <- function(d, upd1, ireset) {
##---------------------------------------------------------
## initialize the diagonal preconditioner (d -- a vector!)
## ---------------------------------------------------------
## global hyk sk yk sr yr yksk yrsr ?? do we have these
# global vectors hyk sk yk sr yr & scalars yksk yrsr
## ---------------------------------------------------------
#% cat("In initpc -- sk and hyk if upd1 false: ", upd1,"\n")
#% print(envjn$sk)
if (upd1) {
td <- d
} else {
if (ireset) {
envjn$hyk <- d * envjn$sk # vector * vector by element??
sds <- as.numeric(crossprod(envjn$sk, envjn$hyk)) # matrix multiply
if (all(envjn$hyk == 0)) { cat("INITPC: envjn$hyk = 0 \n") }
if (sds == 0) { cat("INITPC: sds = 0 \n") }
td <- d - d*d*envjn$sk*envjn$sk/sds + envjn$yk*envjn$yk/envjn$yksk
# by element
} else {
envjn$hyk <- d * envjn$sr # vector * vector
sds <- as.numeric(crossprod(envjn$sr, envjn$hyk))
srds <- as.numeric(crossprod(envjn$sk, envjn$hyk))
yrsk <- as.numeric(crossprod(envjn$yr, envjn$sk))
envjn$hyk <- d*envjn$sk - envjn$hyk*srds/sds + envjn$yr*yrsk/envjn$yrsr
td <- d - d*d*envjn$sr*envjn$sr/sds+envjn$yr*envjn$yr/envjn$yrsr
sds <- as.numeric(crossprod(envjn$sk, envjn$hyk))
td <- td - envjn$hyk*envjn$hyk/sds + envjn$yk*envjn$yk/envjn$yksk
}
}
#% cat("td:")
#% print(td)
ans<-list(td=td)
}
lin1 <- function(p, x, f, alpha, g, sfun, ...){
## ---------------------------------------------------------
## line search (naive)
## ---------------------------------------------------------
## set up
## ---------------------------------------------------------
# cat("lin1: alpha=",alpha," p:\n")
# print(p)
if (is.null(alpha)) alpha <- 0
ierror <- 3
xnew <- x
fnew <- f
gnew <- g
maxit <- 15
if (alpha == 0) {
ierror <- 0
maxit <- 1
}
alpha1 <- alpha
## ---------------------------------------------------------
## line search
## ---------------------------------------------------------
for (itcnt in 1:maxit) {
xt <- x + alpha1*p
fg <- sfun(xt, ...) # Note: added dots 140902
ft<-fg
gt<- attr(fg,"gradient") # may simplify later
if (ft < f) {
ierror <- 0
xnew <- xt
fnew <- ft
gnew <- gt
## cat("about to break in lin1\n")
break
}
alpha1 <- alpha1 / 2
}
if (ierror == 3) { alpha1 <- 0 }
nf1 <- itcnt
## never used in SGN code
## if (nargout == 7) { ## ?? what is nargout?
## dfdp <- as.numeric(crossprod(gt, p) )
## varargout{1} <- dfdp ## ?? what is varargout
## end
result<-list(xnew=xnew, fnew=fnew, gnew=gnew, nf1=nf1,
ierror=ierror, alpha1=alpha1)
}
lmqnbc <- function (x, sfun, lower, upper, maxit, maxfun, stepmx, accrcy, trace, ...) {
## ---------------------------------------------------------
## This is a bounds-constrained truncated-newton method.
## The truncated-newton method is preconditioned by a
## limited-memory quasi-newton method (computed by
## this routine) with a diagonal scaling (routine ndia3).
## For further details, see routine tnbc.
## ---------------------------------------------------------
## global hyk sk yk sr yr yksk yrsr
## ---------------------------------------------------------
## check that initial x is feasible and that the bounds
## are consistent
## ---------------------------------------------------------
n <- length(x)
# JN: Define globals here
gtn<-list(yrsr=0, yksk=0, yr = rep(0, n), yk = rep(0, n),
sr = rep(0, n), sk = rep(0, n),
hg=rep(0,n), hyk=rep(0,n), hyr=rep(0,n) )
envjn<<-list2env(gtn) # Note this important tool for globals in Rtnmin
# end globals
## [ipivot, ierror, x] = crash(x, lower, upper)
cat("lmqnbc -- crout:")
crout<-crash(x, lower, upper)
## print(crout)
ierror <- crout$ierror
ipivot <- crout$ipivot
x<-crout$xnew # in case x changed by bounds
f = 0
g = rep(0,n)
if (ierror != 0) {
stop('LMQNBC: terminating (no feasible point)')
## return
}
## ---------------------------------------------------------
## initialize variables, parameters, and constants
## ---------------------------------------------------------
options(digits=5)
## cat("stempmx, accrcy, maxfun:",stepmx, accrcy, maxfun, "\n")
cat(' it nf cg f |g|\n')
eps <- .Machine$double.eps
upd1 <- TRUE
ncg <- 0
conv <- FALSE
xnorm <- max(abs(x))
ierror <- 0
if ( (stepmx < sqrt(accrcy)) || (maxfun < 1) ) {
ierror <- -1
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror,
nfngr=ncg)
cat("Exiting lmqnbc - almgn:")
print(almqn)
return(almqn) # return
}
## ---------------------------------------------------------
## compute initial function value and related information
## ---------------------------------------------------------
# cat("Try initial fn\n")
fg<- sfun(x, ...)
nf <- 1
nit <- 0
g<- attr(fg,"gradient")
f<-fg
flast <- f
gnorm <- max(abs(g)) ## norm(g,'inf')
## ---------------------------------------------------------
## Test if Lagrange multipliers are non-negative.
## Because the constraints are only bounds, the Lagrange
## multipliers are components of the gradient.
## Then form the projected gradient.
## ---------------------------------------------------------
ind <- which((ipivot != 2) &
(as.numeric(crossprod(ipivot,g)) >0 ) )
if (length(ind) > 0) {
ipivot[ind] <- rep(0, length(ind))
}
g <- ztime (g, ipivot)
gnorm <- max(abs(g))
cat(nit,"\t", nf,"\t", ncg,"\t", f," ", gnorm,"\n")
## print(x)
## print(g)
## ---------------------------------------------------------
## check if the initial point is a local minimum.
## ---------------------------------------------------------
ftest <- 1 + abs(f)
if (gnorm < .01*sqrt(eps)*ftest) {
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror,
nfngr=ncg)
return(almqn)
}
## cat("before set initial, flast=",flast,"\n")
## ---------------------------------------------------------
## set initial values to other parameters
## ---------------------------------------------------------
icycle <- n-1
ireset <- 0
bounds <- TRUE
difnew <- 0
epsred <- .05
fkeep <- f
d <- rep(1,n)
## ---------------------------------------------------------
## ..........main iterative loop..........
## ---------------------------------------------------------
## compute the search direction
## ---------------------------------------------------------
argvec <- c(accrcy, gnorm, xnorm)
##[p, gtp, ncg1, d] <- ...
## modlnp (d, x, g, maxit, upd1, ireset, bounds, ipivot, ## argvec, sfun)
mres <- modlnp (d, x, g, maxit, upd1, ireset,
bounds, ipivot, argvec, sfun, ...)
ncg1 <- mres$ncg1
gtp <- mres$gtp
d <- mres$dnew
p <- mres$p
## cat("p:")
## print(p)
## tmp<-readline("cont.")
ncg <- ncg + ncg1
while (!conv) {
oldg <- g
pnorm <- max(abs(p)) # norm2(p, 'inf')
oldf <- f
## ---------------------------------------------------------
## line search
## ---------------------------------------------------------
pe <- pnorm + eps
spe <- stpmax (stepmx, pe, x, p, ipivot, lower, upper)
## cat("spe=",spe,"\n")
alpha <- step1 (f, gtp, spe)
alpha0 <- alpha
## cat("alpha0=",alpha0,"\n")
## [x_new, f_new, g_new, nf1, ierror, alpha] <- lin1 (p, x,
## f, alpha0, g, sfun)
reslin <- lin1 (p, x, f, alpha0, g, sfun, ...)
ierror <- reslin$ierror
alpha <- reslin$alpha1
## ---------------------------------------------------------
if ((alpha == 0) && (alpha0 != 0) || (ierror == 3)){
cat('Error in Line Search\n')
cat(' ierror = ', ierror, "\n")
cat(' alpha = ',alpha, "\n")
cat(' alpha0 = ', alpha0, "\n")
cat(' gtp = ', gtp, "\n")
## ############################
cat(' |g| = ', norm2(g), "\n")
cat(' |p| = ', norm2(p), "\n")
tmp <- readline('Hit any key to continue')
## ############################
}
## #######################
x <- reslin$xnew # need fixup
f <- reslin$fnew
g <- reslin$gnew
nf1 <- reslin$nf1
nf <- nf + nf1
nit <- nit + 1
## ---------------------------------------------------------
## update active set, if appropriate
## ---------------------------------------------------------
newcon <- FALSE
## cat("Check active set - alpha, spe, f:",alpha, spe,f,"\n")
if (abs(alpha-spe) <= 10*eps) {
newcon <- TRUE
ierror <- 0
## cat("flast:",flast,"\n")
# [ipivot, flast] <- modz (x, p, ipivot, lower, upper,
# flast, f, alpha)
## cat("ipivot before and after update by modz:\n")
## print(ipivot)
modzres<-modz(x, p, ipivot, lower, upper,
flast, f, alpha)
ipivot <- modzres$ipivot1
flast<-modzres$flast1
## print(ipivot)
}
if (ierror == 3) {
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror,
nfngr=ncg)
cat("Error updating active set -- almqn:")
print(almqn)
return(almqn)
}
## ---------------------------------------------------------
## stop if more than maxfun evaluations have been made
## ---------------------------------------------------------
if (nf > maxfun) {
ierror <- 2
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror,
nfngr=ncg)
cat("Too many function evaluations -- almqn:")
print(almqn)
return(almqn)
}
## ---------------------------------------------------------
## set up for convergence and resetting tests
## ---------------------------------------------------------
difold <- difnew
difnew <- oldf - f # scalar
if (icycle == 1) {
if (difnew > 2*difold) { epsred <- 2*epsred }
if (difnew < .5*difold) { epsred <- .5*epsred }
}
gv <- ztime (g, ipivot)
gnorm <- max(abs(gv))
ftest <- 1 + abs(f)
xnorm <- max(abs(x))
############### DISPLAY ##############
cat(nit,"\t", nf,"\t", ncg,"\t", f," ", gnorm,"\n")
## print(x)
## print(g)
## ---------------------------------------------------------
## test for convergence
## ---------------------------------------------------------
## [conv, flast, ipivot] <- cnvtst (alpha, pnorm, xnorm, ...
## difnew, ftest, gnorm, gtp, f, flast, g, ...
## ipivot, accrcy)
## cat("before cnvtst, flast=",flast,"\n")
ctres <- cnvtst (alpha, pnorm, xnorm,
difnew, ftest, gnorm, gtp, f, flast, g,
ipivot, accrcy)
conv <- ctres$conv
flast <- ctres$flast1
ipivot <- ctres$ipivot1
## cat("after cnvtst - conv, flast, ipivot:", conv, flast,"\n")
## print(ipivot)
if (conv) {
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror,
nfngr=ncg)
return(almqn)
}
g <- ztime (g, ipivot)
## ---------------------------------------------------------
## modify data for LMQN preconditioner
## ---------------------------------------------------------
if (! newcon) { ## 0 is FALSE and 1 is TRUE supposedly
envjn$yk <- g - oldg
envjn$sk <- alpha*p
envjn$yksk <- as.numeric(crossprod(envjn$yk, envjn$sk))
## cat("yksk:",envjn$yksk,"\n")
ireset <- ( (icycle == n-1) |
(difnew < epsred*(fkeep-f)) )
if (! ireset) {
envjn$yrsr <-
as.numeric(crossprod(envjn$yr,envjn$sr))
ireset <- (envjn$yrsr <= 0)
}
upd1 <- (envjn$yksk <= 0)
}
## cat("newcon, upd1:", newcon, upd1,"\n")
## ---------------------------------------------------------
## compute the search direction
## ---------------------------------------------------------
argvec <- c(accrcy, gnorm, xnorm)
## [p, gtp, ncg1, d] <- ...
## modlnp (d, x, g, maxit, upd1, ireset, bounds,
## ipivot, argvec, sfun)
mres <- modlnp (d, x, g, maxit, upd1, ireset,
bounds, ipivot, argvec, sfun, ...)
ncg1 <- mres$ncg1
gtp <- mres$gtp
d <- mres$dnew
p <- mres$p
## cat("New p:")
## print(p)
## tmp<-readline("cont.")
ncg <- ncg + ncg1
## ---------------------------------------------------------
## update LMQN preconditioner
## ---------------------------------------------------------
if (! newcon) {
if (ireset) {
envjn$sr <- envjn$sk
envjn$yr <- envjn$yk
fkeep <- f
icycle <- 1
} else {
envjn$sr <- envjn$sr + envjn$sk
envjn$yr <- envjn$yr + envjn$yk
icycle <- icycle + 1
}
}
} # end while !conv
}
lmqn <- function (x, sfun, maxit, maxfun, stepmx, accrcy, trace, ...) {
## ---------------------------------------------------------
## truncated-newton method for unconstrained minimization
## (customized version)
## ---------------------------------------------------------
# global vectors hyk sk yk sr yr & scalars yksk yrsr
## ---------------------------------------------------------
## set up
## ---------------------------------------------------------
## format compact
## format short e
n<-length(x)
# JN: Define globals here. Is it necessary to set up.
gtn<-list(yrsr=0, yksk=0, yr = rep(0, n), yk = rep(0, n), sr = rep(0, n), sk = rep(0, n))
envjn<<-list2env(gtn)
# end globals
eps <- .Machine$double.eps
upd1 <- 1
ncg <- 0
xnorm <- max(abs(x)) # norm(x,'inf')
ierror <- 0
if (stepmx < sqrt(accrcy) || maxfun < 1) {
ierror <- -1
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror,
nfngr=ncg)
return(almqn)
}
## ---------------------------------------------------------
## compute initial function value and related information
## ---------------------------------------------------------
fg <- sfun(x, ...)
#% print(fg)
g<- attr(fg, "gradient")
f<-fg
gnorm <- max(abs(g)) ## norm(g,'inf')
nf <- 1
nit <- 0
if (trace) cat("Itn ",nit," ",nf," ",ncg, " ",f, " ", gnorm,"\n")
## ---------------------------------------------------------
## check for small gradient at the starting point.
## ---------------------------------------------------------
ftest <- 1 + abs(f)
if (gnorm < .01*sqrt(eps)*ftest) {
ierror <- 0
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror, nfngr=ncg)
return(almqn)
}
## ---------------------------------------------------------
## set initial values to other parameters
## ---------------------------------------------------------
n <- length(x)
icycle <- n-1
toleps <- sqrt(accrcy) + sqrt(eps)
rtleps <- accrcy + eps
difnew <- 0
epsred <- .05
fkeep <- f
conv <- FALSE
ireset <- 0
ipivot <- 0
## ---------------------------------------------------------
## initialize diagonal preconditioner to the identity
## ---------------------------------------------------------
d <- rep(1,n) # as a vector
## ---------------------------------------------------------
## ..........main iterative loop..........
## ---------------------------------------------------------
## compute search direction
## ---------------------------------------------------------
argvec <- c(accrcy, gnorm, xnorm)
mres <- modlnp (d, x, g, maxit, upd1, ireset, bounds=FALSE, ipivot, argvec, sfun, ...)
p <- mres$p
# cat("p from first call to modlnp\n")
# print(p)
# tmp<-readline("cont.")
gtp <- mres$gtp
ncg1<-mres$ncg1
d <- mres$dnew
#% cat("d from modlnp:")
#% print(d)
ncg <- ncg + ncg1
while ( ! conv) {
#% cat("top while ")
#% tmp <- readline("Top of iteration")
oldg <- g
pnorm <- max(abs(p)) # norm(p,'inf')
oldf <- f
## ---------------------------------------------------------
## line search
## ---------------------------------------------------------
pe <- pnorm + eps
spe <- stepmx/pe
#% cat("gtp, spe:", gtp, spe,"\n")
alpha0 <- step1 (f, gtp, spe)
reslin <- lin1 (p, x, f, alpha0, g, sfun, ...)
## [x, f, g, nf1, ierror, alpha] <-
x <- reslin$xnew # need fixup
f <- reslin$fnew
g <- reslin$gnew
nf1 <- reslin$nf1
ierror <- reslin$ierror
alpha <- reslin$alpha1
# cat("after lin1, alpha=",alpha,"\n")
# tmp<-readline("cont.")
nf <- nf + nf1
## ---------------------------------------------------------
nit <- nit + 1
gnorm <- max(abs(g)) # norm(g,'inf')
### Display info
if (trace) cat("Itn ",nit," ",nf," ",ncg, " ",f, " ", gnorm,"\n")
if (ierror == 3) {
if (length(ncg) == 0) { ncg <- 0 } # ?? is.null(ncg)??
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror, nfngr=ncg)
return(almqn)
}
## ---------------------------------------------------------
## stop if more than maxfun evalutations have been made
## ---------------------------------------------------------
if (nf >= maxfun) {
ierror <- 2
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror, nfngr=ncg)
return(almqn)
}
## ---------------------------------------------------------
## set up for convergence and resetting tests
## ---------------------------------------------------------
ftest <- 1 + abs(f)
xnorm <- max(abs(x)) # norm(x,'inf')
difold <- difnew
difnew <- oldf - f
envjn$yk <- g - oldg
envjn$sk <- alpha*p
if (icycle == 1) {
if (difnew > 2*difold) { epsred <- 2*epsred }
if (difnew < 0.5*difold) { epsred <- 0.5*epsred }
}
## ---------------------------------------------------------
## convergence test
## ---------------------------------------------------------
conv <- ( ( (alpha*pnorm < toleps*(1 + xnorm)) &&
(abs(difnew) < rtleps*ftest) &&
(gnorm < accrcy^(1/3)*ftest) )||
( gnorm < .01*sqrt(accrcy)*ftest ) )
if (conv) {
ierror <- 0
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror, nfngr=ncg)
return(almqn)
}
## ---------------------------------------------------------
## update lmqn preconditioner
## ---------------------------------------------------------
envjn$yksk <- as.numeric(crossprod(envjn$yk, envjn$sk) )
ireset <- ((icycle == n-1) || (difnew < epsred*(fkeep-f)) )
if ( ! ireset) {
envjn$yrsr <- as.numeric(crossprod(envjn$yr, envjn$sr))
ireset <- (envjn$yrsr <= 0)
}
upd1 <- (envjn$yksk <= 0)
## ---------------------------------------------------------
## compute search direction
## ---------------------------------------------------------
argvec <- c(accrcy, gnorm, xnorm)
## cat("ireset =", ireset," upd1=",upd1,"\n")
#% cat("New d to modlnp:")
#% print(d)
## [p, gtp, ncg1, d] <- ...
mres <- modlnp (d, x, g, maxit, upd1, ireset, 0, ipivot, argvec, sfun, ...)
p <- mres$p
gtp <- mres$gtp
ncg1<-mres$ncg1
d <- mres$dnew
#% cat("returned dnew, envjn$john: ", envjn$john)
#% print(d)
ncg <- ncg + ncg1
## ---------------------------------------------------------
## store information for lmqn preconditioner
## ---------------------------------------------------------
if (ireset) {
envjn$sr <- envjn$sk
envjn$yr <- envjn$yk
fkeep <- f
icycle <- 1
} else {
envjn$sr <- envjn$sr + envjn$sk
envjn$yr <- envjn$yr + envjn$yk
icycle <- icycle + 1
}
} # end while
## [xstar, f, g, ierror] = ..
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror, nfngr=ncg)
} # end lmqn
modlnp <- function(d, x, g, maxit, upd1, ireset, bounds,
ipivot, argvec, sfun, ...) {
##---------------------------------------------------------
## this routine performs a preconditioned conjugate-gradient
## iteration to solve the Newton equations for a search
## direction for a truncated-newton algorithm.
## When the value of the quadratic model is sufficiently
## reduced, the iteration is terminated.
##---------------------------------------------------------
## parameters
##
## p - computed search direction
## g - current gradient
## gv,gz1,v - scratch vectors
## r - residual
## d - diagonal preconditoning matrix
## feval - value of quadratic function
##------------------------------------------------------------
## initialization
##------------------------------------------------------------
if (is.null(d)) stop("Null d")
## print(x)
## print(g)
## cat(bounds,"\n")
accrcy <- argvec[1]
gnorm <- argvec[2]
xnorm <- argvec[3]
if (maxit == 0) {
p <- -g
gtp <- as.numeric(crossprod(p, g))
ncg1 <- 1
dnew <- d
if (sqrt(sum(p^2))==0) {
## cat("MODLNP 01: |p| = 0\n")
## pause(1)
}
#% cat("modlnp - dout01:")
#% print(dnew)
result<-list(p=p, gtp=gtp, ncg1=ncg1, dnew=dnew)
return(result)
}
first <- 1
tol <- 1e-6 ######### was 1.d-12 #########
qold <- 0
ainit <- initpc (d, upd1, ireset)
dnew<-ainit$td
#% cat("after initpc dnew:")
#% print(dnew)
r <- -g
v <- rep(0, length(r))
p <- v
gtp <- as.numeric(crossprod(p, g))
rho <- rep(0, maxit+1)
beta <- rep(0, maxit)
v.gv <- rep(0, maxit)
rho[[1]] <- as.numeric(crossprod(r))
##------------------------------------------------------------
## main iteration (conjugate-gradient iterations for Ax = b)
##------------------------------------------------------------
ind <- 0
ncg1 <- 0
for (k in 1:maxit) {
ncg1 <- ncg1 + 1
if (bounds) { r <- ztime(r, ipivot) }
amsolve <- msolve (r, upd1, ireset, first, d)
zk<-amsolve$y
if (bounds) { zk <- ztime (zk, ipivot) }
#% cat("r:")
#% print(r)
#% cat("zk:")
#% print(zk)
rz <- as.numeric(crossprod(r,zk))
if (rz/gnorm < tol) {
ind <- 80
if (sqrt(sum(p^2))==0) {
p <- -g
gtp <- as.numeric(crossprod(p, g))
}
#% cat("modlnp - dout - ind80:")
#% print(dnew)
result<-list(p=p, gtp=gtp, ncg1=ncg1, dnew=dnew)
return(result)
}
if (k > 1) {
beta[k] <- rz/rzold
} else {
beta[k] <- 0
}
#% cat("beta[",k,"] =",beta[k],"\n")
v <- zk + beta[k]*v
if (bounds) { v <- ztime( v, ipivot) }
## cat("about to call gtims\n")
gv <- gtims(v, x, g, accrcy, xnorm, sfun, ...)
#% cat("After gtims: gv=")
#% print(gv)
## cat(bounds,"\n")
if (bounds) { gv <- ztime (gv, ipivot) }
#% cat("gv, v:")
## cat("gv=", gv, "\n")
#% print(gv)
#% print(v)
v.gv[[k]] <- as.numeric(crossprod(v, gv)) ## ?? v.gv has underscore!! ??
#% cat("v.gv[[",k,"]]=",v.gv[[k]],"\n")
if (v.gv[[k]]/gnorm < tol) {
ind <- 50
if (sqrt(sum(p^2))==0) {
p <- -g ## Mod SGN 140912
if (bounds) {
p <- ztime(p, ipivot)
}
gtp <- crossprod(p, g)
## cat("MODLNP 03: |p| = 0 \n")
## pause(1)
}
#% cat("modlnp - dout03:")
#% print(dnew)
result<-list(p=p, gtp=gtp, ncg1=ncg1, dnew=dnew)
return(result)
}
dnew <- ndia3(dnew, v, gv, r, v.gv[[k]]) # NB need to return something but it doesn't get used
#% cat("after ndia3 below dout03, dnew:")
#% print(dnew)
##------------------------------------------------------------
## compute current solution and related vectors
##------------------------------------------------------------
alpha <- rz / v.gv[[k]]
p <- p + alpha* v
r <- r - alpha*gv
rho[[k+1]] <- as.numeric(crossprod(r))
##------------------------------------------------------------
## test for convergence
##------------------------------------------------------------
gtp <- as.numeric(crossprod(p,g))
pr <- as.numeric(crossprod(r,p))
q <- (gtp + pr) / 2
qtest <- k * (1 - qold/q)
if (qtest <= 0.5) {
if (sqrt(sum(p^2))==0) {
## cat("MODLNP 04: |p| = 0\n")
## pause(1)
}
#% cat("modlnp - dout04:")
#% print(dnew)
result<-list(p=p, gtp=gtp, ncg1=ncg1, dnew=dnew)
return(result)
}
qold <- q
##------------------------------------------------------------
## perform cautionary test
##------------------------------------------------------------
if (gtp > 0) {
ind <- 40
if (sqrt(sum(p^2))==0) {
## cat("MODLNP 05: |p| = 0 \n")
## pause(1)
}
#% cat("modlnp - dout05:")
#% print(dnew)
result<-list(p=p, gtp=gtp, ncg1=ncg1, dnew=dnew)
return(result)
}
rzold <- rz
} ## end loop over k
k <- k-1
##------------------------------------------------------------
## terminate algorithm
##------------------------------------------------------------
if (ind == 40) {
p <- p - alpha*v
}
##------------------------------------------------------------
if (ind == 50 && k <= 1) {
amsolve <- msolve (g, upd1, ireset, first, d)
p<-amsolve$y
p <- -p
if (bounds) { p <- ztime (p, ipivot) }
}
##------------------------------------------------------------
if (ind == 80 && k <= 1) {
p <- -g
if (bounds) { p <- ztime (p, ipivot) }
}
##------------------------------------------------------------
## store new diagonal preconditioner
##------------------------------------------------------------
gtp <- as.numeric(crossprod(p, g))
ncg1 <- k + 1
if (sqrt(sum(p^2))==0) {
## cat("MODLNP 06: |p| = 0 \n")
## pause(1)
}
## cat("modlnp - dout06:")
#% print(dnew)
result<-list(p=p, gtp=gtp, ncg1=ncg1, dnew=dnew)
return(result)
}
modz <- function (x, p, ipivot, low, up, flast, f, alpha){
##---------------------------------------------------------------------
## update the constraint matrix if a new constraint is encountered
##---------------------------------------------------------------------
eps <- .Machine$double.eps
indl <- which(ipivot == 0 & p < 0)
if (length(indl) > 0) {
toll <- 10 * eps * (abs(low[indl]) + 1)
hitl <- which(x[indl]-low[indl] <= toll)
if (length(hitl) > 0) {
flast <- f
ipivot[indl[hitl]] <- -1
}
}
##---------------------------------------------------------------------
indu <- which((ipivot == 0) & (p > 0));
if (length(indu) > 0) {
tolu <- 10 * eps * (abs( up[indu]) + 1)
hitu <- which(up[indu]-x[indu] <= tolu)
if (length(hitu) > 0) {
flast <- f
ipivot[indu[hitu]] <- 1
}
}
##---------------------------------------------------------------------
flast1 <- flast
ipivot1 <- ipivot
list(flast1 = flast, ipivot1=ipivot)
}
msolve <- function(g, upd1, ireset, first, d) {
##---------------------------------------------------------
## This routine acts as a preconditioning step for the
## linear conjugate-gradient routine. It is also the
## method of computing the search direction from the
## gradient for the non-linear conjugate-gradient code.
## It represents a two-step self-scaled bfgs formula.
##---------------------------------------------------------
#% cat("msolve, envjn$john=",envjn$john,"\n")
envjn$john<-"msolve"
if (upd1) {
#% cat("upd1 is TRUE\n")
y <- g/d
} else {
gsk <- as.numeric(crossprod(g, envjn$sk))
if (ireset) {
envjn$hg <- g/d
if (first) {
envjn$hyk <- envjn$yk/d
ykhyk <- as.numeric(crossprod(envjn$yk, envjn$hyk))
}
ghyk <- as.numeric(crossprod(g, envjn$hyk))
y <- ssbfgs(envjn$sk,envjn$hg,envjn$hyk,envjn$yksk,ykhyk,gsk,ghyk)
} else {
envjn$hg <- g/d
if (first) {
envjn$hyk <- envjn$yk/d
envjn$hyr <- envjn$yr/d
envjn$yksr <- as.numeric(crossprod(envjn$yk, envjn$sr))
ykhyr <- as.numeric(crossprod(envjn$yk, envjn$hyr))
}
gsr <- as.numeric(crossprod(g, envjn$sr))
ghyr <- as.numeric(crossprod(g, envjn$hyr))
if (first) {
yrhyr <- as.numeric(crossprod(envjn$yr, envjn$hyr))
}
envjn$hg <- ssbfgs(envjn$sr,envjn$hg,envjn$hyr,envjn$yrsr,yrhyr,gsr,ghyr)
if (first) {
envjn$hyk <- ssbfgs(envjn$sr,envjn$hyk,envjn$hyr,envjn$yrsr,yrhyr,envjn$yksr,ykhyr)
}
ykhyk <- as.numeric(crossprod(envjn$hyk, envjn$yk))
ghyk <- as.numeric(crossprod(envjn$hyk, g))
y <- ssbfgs(envjn$sk,envjn$hg,envjn$hyk,envjn$yksk,ykhyk,gsk,ghyk)
}
}
amsolve<-list(y=y) ## returns y
}
ndia3<-function(e, v, gv, r, vgv){
##---------------------------------------------------------
## update the preconditioning matrix based on a diagonal
## version of the bfgs quasi-newton update.
##---------------------------------------------------------
tol <- 1e-6
vr <- as.numeric(crossprod(v,r))
#% cat("ndia3: vgv, vr, then e:",vgv, vr,"\n")
#% print(e)
if (abs(vr)>tol && abs(vgv)>tol) {
e <- e - (r*r)/vr + (gv*gv)/vgv # CAUTION! May be crossprod
ind <- which(e < tol)
if (length(ind)>0){
e[ind] <- 1
}
}
#% print(e)
e # Need to return the object?? Possibly not!
}
norm2 <- function (x) {
##---------------------------------------------------------
## compute the 2 norm of vector x
##---------------------------------------------------------
res<-sqrt(as.numeric(crossprod(x)))
}
ssbfgs <- function (s, hv, hy, ys, yhy, vs, vhy){
##---------------------------------------------------------
## self-scaled bfgs quasi-Newton update
## (used by preconditioner)
##---------------------------------------------------------
delta <- (1 + yhy/ys)*vs/ys - vhy/ys
beta <- -vs/ys
z <- hv + delta*s + beta*hy
}
step1 <- function(f, gtp, smax) {
##---------------------------------------------------------
## step1 returns the length of the initial step to be
## taken along the vector p in the next linear search.
##---------------------------------------------------------
## [fm is supposed to be an estimate of the optimal function value]
eps<-.Machine$double.eps
fm <- 0
d <- abs(f-fm)
alpha <- min(1, smax)
if ((2*d <= (-gtp)) && (d >= eps)) {
alpha = min(-2*d/gtp, smax)
}
alpha # need to ensure returned value
}
stpmax <- function(stepmx, pe, x, p, ipivot, low, up) {
##------------------------------------------------
## compute the maximum allowable step length
## (spe is the standard (unconstrained) max step)
##------------------------------------------------
## cat("stpmax: stepmx, pe:",stepmx, pe,"\n")
spe <- stepmx / pe
al <- spe
au <- spe
##------------------------------------------------
indl <- which(ipivot==0 & p < 0)
if ( length(indl)>0) {
tl <- low[indl] - x[indl]
al <- min(tl/p[indl])
}
##------------------------------------------------
indu <- which(ipivot==0 & p > 0)
if (length(indu)>0) {
tu <- up[indu] - x[indu]
au <- min(tu/p[indu])
}
##------------------------------------------------
spe <- min(c(spe, al, au))
}
ztime<-function(x,ipivot) {
## ---------------------------------------------------------
## this routine multiplies the vector x by the
## constraint matrix z
## ---------------------------------------------------------
ind <- which(ipivot!=0);
x[ind] <- 0
x1 <- x;
}
Try the Rtnmin package in your browser
Any scripts or data that you put into this service are public.
Rtnmin documentation built on May 2, 2019, 7:22 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions initpc lin1 modlnp msolve ndia3 step1 stpmax ztime
cnvtst <- function (alpha, pnorm, xnorm,
dif, ftest, gnorm, gtp, f, flast, g, ipivot, accrcy) {
## ---------------------------------------------------------
## test for convergence
## ---------------------------------------------------------
## set up
## ---------------------------------------------------------
conv <- 0;
eps <- .Machine$double.eps
toleps <- sqrt(accrcy) + sqrt(eps);
rtleps <- accrcy + eps;
imax <- 0;
ltest <- (flast - f <= -0.5*gtp);
## ---------------------------------------------------------
## if anti-zigzag test satisfied, test multipliers;
## if appropriate, modify active set
## ---------------------------------------------------------
if ( ! ltest) {
ind <- which( (ipivot != 0) & (ipivot != 2))
if ( length(ind) > 0 ) { ## how to ensure ind not empty
t <- -sum(ipivot[ind]*g[ind])
cmax <- min(t) # ?? why [cmax, imax] = min(t) ??
imax <- cmax
if (cmax >= 0) { imax <- 0 }
}
}
if (imax != 0) {
ipivot[ ind[imax] ] <- 0;
flast <- f;
} else {
conv = ( ( (alpha*pnorm < toleps*(1 + xnorm) )
&& (abs(dif) < rtleps*ftest)
&& (gnorm < accrcy^(1/3)*ftest) ) ||
( gnorm < .01*sqrt(accrcy)*ftest) )
}
result <- list(conv=conv, flast1=flast, ipivot1=ipivot)
}
crash <- function (x, low, up) {
##---------------------------------------------------------
## this initializes the constraint information, and
## ensures that the initial point satisfies
## low <= x <= up.
## the constraints are checked for consistency.
##---------------------------------------------------------
ierror <- 0
if (any(low > up)) { ierror = - max(which(low > up)) }
# above is check on error in bounds specification
xnew <- pmax (low, x) # force params into bounds
xnew <- pmin (up, xnew) # No diagnostic! ??
# we output revised parameters
n <- length(x)
ind <- which(low == x)
ipivot <- rep(0, n) ## zeros(size(x)) in Matlab
if (length(ind) > 0) { ipivot[ind] <- -1 }
ind <- which(x == up)
if (length(ind) > 0) { ipivot[ind] <- 1 }
ind <- which(low == up)
if (length(ind) > 0) { ipivot[ind] <- 2 }
list(ipivot=ipivot, ierror=ierror, xnew=xnew)
}
gtims <- function (v, x, g, accrcy, xnorm, sfun, ...) {
##---------------------------------------------------------
## compute the product of the Hessian times the vector v;
## store result in the vector gv
## (finite-difference version)
##---------------------------------------------------------
## cat("v, x, g:")
## print(v)
## print(x)
## print(g)
delta <- sqrt(accrcy)*(1 + xnorm)/sqrt(sum(v^2))
hg <- x + delta*v
tresult<-sfun(hg, ...)
## cat("tresult: ")
## print(tresult)
gv <- attr(tresult, "gradient")
gv <- (gv - g)/delta
gv
}
initpc <- function(d, upd1, ireset) {
##---------------------------------------------------------
## initialize the diagonal preconditioner (d -- a vector!)
## ---------------------------------------------------------
## global hyk sk yk sr yr yksk yrsr ?? do we have these
# global vectors hyk sk yk sr yr & scalars yksk yrsr
## ---------------------------------------------------------
#% cat("In initpc -- sk and hyk if upd1 false: ", upd1,"\n")
#% print(envjn$sk)
if (upd1) {
td <- d
} else {
if (ireset) {
envjn$hyk <- d * envjn$sk # vector * vector by element??
sds <- as.numeric(crossprod(envjn$sk, envjn$hyk)) # matrix multiply
if (all(envjn$hyk == 0)) { cat("INITPC: envjn$hyk = 0 \n") }
if (sds == 0) { cat("INITPC: sds = 0 \n") }
td <- d - d*d*envjn$sk*envjn$sk/sds + envjn$yk*envjn$yk/envjn$yksk
# by element
} else {
envjn$hyk <- d * envjn$sr # vector * vector
sds <- as.numeric(crossprod(envjn$sr, envjn$hyk))
srds <- as.numeric(crossprod(envjn$sk, envjn$hyk))
yrsk <- as.numeric(crossprod(envjn$yr, envjn$sk))
envjn$hyk <- d*envjn$sk - envjn$hyk*srds/sds + envjn$yr*yrsk/envjn$yrsr
td <- d - d*d*envjn$sr*envjn$sr/sds+envjn$yr*envjn$yr/envjn$yrsr
sds <- as.numeric(crossprod(envjn$sk, envjn$hyk))
td <- td - envjn$hyk*envjn$hyk/sds + envjn$yk*envjn$yk/envjn$yksk
}
}
#% cat("td:")
#% print(td)
ans<-list(td=td)
}
lin1 <- function(p, x, f, alpha, g, sfun, ...){
## ---------------------------------------------------------
## line search (naive)
## ---------------------------------------------------------
## set up
## ---------------------------------------------------------
# cat("lin1: alpha=",alpha," p:\n")
# print(p)
if (is.null(alpha)) alpha <- 0
ierror <- 3
xnew <- x
fnew <- f
gnew <- g
maxit <- 15
if (alpha == 0) {
ierror <- 0
maxit <- 1
}
alpha1 <- alpha
## ---------------------------------------------------------
## line search
## ---------------------------------------------------------
for (itcnt in 1:maxit) {
xt <- x + alpha1*p
fg <- sfun(xt, ...) # Note: added dots 140902
ft<-fg
gt<- attr(fg,"gradient") # may simplify later
if (ft < f) {
ierror <- 0
xnew <- xt
fnew <- ft
gnew <- gt
## cat("about to break in lin1\n")
break
}
alpha1 <- alpha1 / 2
}
if (ierror == 3) { alpha1 <- 0 }
nf1 <- itcnt
## never used in SGN code
## if (nargout == 7) { ## ?? what is nargout?
## dfdp <- as.numeric(crossprod(gt, p) )
## varargout{1} <- dfdp ## ?? what is varargout
## end
result<-list(xnew=xnew, fnew=fnew, gnew=gnew, nf1=nf1,
ierror=ierror, alpha1=alpha1)
}
lmqnbc <- function (x, sfun, lower, upper, maxit, maxfun, stepmx, accrcy, trace, ...) {
## ---------------------------------------------------------
## This is a bounds-constrained truncated-newton method.
## The truncated-newton method is preconditioned by a
## limited-memory quasi-newton method (computed by
## this routine) with a diagonal scaling (routine ndia3).
## For further details, see routine tnbc.
## ---------------------------------------------------------
## global hyk sk yk sr yr yksk yrsr
## ---------------------------------------------------------
## check that initial x is feasible and that the bounds
## are consistent
## ---------------------------------------------------------
n <- length(x)
# JN: Define globals here
gtn<-list(yrsr=0, yksk=0, yr = rep(0, n), yk = rep(0, n),
sr = rep(0, n), sk = rep(0, n),
hg=rep(0,n), hyk=rep(0,n), hyr=rep(0,n) )
envjn<<-list2env(gtn) # Note this important tool for globals in Rtnmin
# end globals
## [ipivot, ierror, x] = crash(x, lower, upper)
cat("lmqnbc -- crout:")
crout<-crash(x, lower, upper)
## print(crout)
ierror <- crout$ierror
ipivot <- crout$ipivot
x<-crout$xnew # in case x changed by bounds
f = 0
g = rep(0,n)
if (ierror != 0) {
stop('LMQNBC: terminating (no feasible point)')
## return
}
## ---------------------------------------------------------
## initialize variables, parameters, and constants
## ---------------------------------------------------------
options(digits=5)
## cat("stempmx, accrcy, maxfun:",stepmx, accrcy, maxfun, "\n")
cat(' it nf cg f |g|\n')
eps <- .Machine$double.eps
upd1 <- TRUE
ncg <- 0
conv <- FALSE
xnorm <- max(abs(x))
ierror <- 0
if ( (stepmx < sqrt(accrcy)) || (maxfun < 1) ) {
ierror <- -1
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror,
nfngr=ncg)
cat("Exiting lmqnbc - almgn:")
print(almqn)
return(almqn) # return
}
## ---------------------------------------------------------
## compute initial function value and related information
## ---------------------------------------------------------
# cat("Try initial fn\n")
fg<- sfun(x, ...)
nf <- 1
nit <- 0
g<- attr(fg,"gradient")
f<-fg
flast <- f
gnorm <- max(abs(g)) ## norm(g,'inf')
## ---------------------------------------------------------
## Test if Lagrange multipliers are non-negative.
## Because the constraints are only bounds, the Lagrange
## multipliers are components of the gradient.
## Then form the projected gradient.
## ---------------------------------------------------------
ind <- which((ipivot != 2) &
(as.numeric(crossprod(ipivot,g)) >0 ) )
if (length(ind) > 0) {
ipivot[ind] <- rep(0, length(ind))
}
g <- ztime (g, ipivot)
gnorm <- max(abs(g))
cat(nit,"\t", nf,"\t", ncg,"\t", f," ", gnorm,"\n")
## print(x)
## print(g)
## ---------------------------------------------------------
## check if the initial point is a local minimum.
## ---------------------------------------------------------
ftest <- 1 + abs(f)
if (gnorm < .01*sqrt(eps)*ftest) {
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror,
nfngr=ncg)
return(almqn)
}
## cat("before set initial, flast=",flast,"\n")
## ---------------------------------------------------------
## set initial values to other parameters
## ---------------------------------------------------------
icycle <- n-1
ireset <- 0
bounds <- TRUE
difnew <- 0
epsred <- .05
fkeep <- f
d <- rep(1,n)
## ---------------------------------------------------------
## ..........main iterative loop..........
## ---------------------------------------------------------
## compute the search direction
## ---------------------------------------------------------
argvec <- c(accrcy, gnorm, xnorm)
##[p, gtp, ncg1, d] <- ...
## modlnp (d, x, g, maxit, upd1, ireset, bounds, ipivot, ## argvec, sfun)
mres <- modlnp (d, x, g, maxit, upd1, ireset,
bounds, ipivot, argvec, sfun, ...)
ncg1 <- mres$ncg1
gtp <- mres$gtp
d <- mres$dnew
p <- mres$p
## cat("p:")
## print(p)
## tmp<-readline("cont.")
ncg <- ncg + ncg1
while (!conv) {
oldg <- g
pnorm <- max(abs(p)) # norm2(p, 'inf')
oldf <- f
## ---------------------------------------------------------
## line search
## ---------------------------------------------------------
pe <- pnorm + eps
spe <- stpmax (stepmx, pe, x, p, ipivot, lower, upper)
## cat("spe=",spe,"\n")
alpha <- step1 (f, gtp, spe)
alpha0 <- alpha
## cat("alpha0=",alpha0,"\n")
## [x_new, f_new, g_new, nf1, ierror, alpha] <- lin1 (p, x,
## f, alpha0, g, sfun)
reslin <- lin1 (p, x, f, alpha0, g, sfun, ...)
ierror <- reslin$ierror
alpha <- reslin$alpha1
## ---------------------------------------------------------
if ((alpha == 0) && (alpha0 != 0) || (ierror == 3)){
cat('Error in Line Search\n')
cat(' ierror = ', ierror, "\n")
cat(' alpha = ',alpha, "\n")
cat(' alpha0 = ', alpha0, "\n")
cat(' gtp = ', gtp, "\n")
## ############################
cat(' |g| = ', norm2(g), "\n")
cat(' |p| = ', norm2(p), "\n")
tmp <- readline('Hit any key to continue')
## ############################
}
## #######################
x <- reslin$xnew # need fixup
f <- reslin$fnew
g <- reslin$gnew
nf1 <- reslin$nf1
nf <- nf + nf1
nit <- nit + 1
## ---------------------------------------------------------
## update active set, if appropriate
## ---------------------------------------------------------
newcon <- FALSE
## cat("Check active set - alpha, spe, f:",alpha, spe,f,"\n")
if (abs(alpha-spe) <= 10*eps) {
newcon <- TRUE
ierror <- 0
## cat("flast:",flast,"\n")
# [ipivot, flast] <- modz (x, p, ipivot, lower, upper,
# flast, f, alpha)
## cat("ipivot before and after update by modz:\n")
## print(ipivot)
modzres<-modz(x, p, ipivot, lower, upper,
flast, f, alpha)
ipivot <- modzres$ipivot1
flast<-modzres$flast1
## print(ipivot)
}
if (ierror == 3) {
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror,
nfngr=ncg)
cat("Error updating active set -- almqn:")
print(almqn)
return(almqn)
}
## ---------------------------------------------------------
## stop if more than maxfun evaluations have been made
## ---------------------------------------------------------
if (nf > maxfun) {
ierror <- 2
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror,
nfngr=ncg)
cat("Too many function evaluations -- almqn:")
print(almqn)
return(almqn)
}
## ---------------------------------------------------------
## set up for convergence and resetting tests
## ---------------------------------------------------------
difold <- difnew
difnew <- oldf - f # scalar
if (icycle == 1) {
if (difnew > 2*difold) { epsred <- 2*epsred }
if (difnew < .5*difold) { epsred <- .5*epsred }
}
gv <- ztime (g, ipivot)
gnorm <- max(abs(gv))
ftest <- 1 + abs(f)
xnorm <- max(abs(x))
############### DISPLAY ##############
cat(nit,"\t", nf,"\t", ncg,"\t", f," ", gnorm,"\n")
## print(x)
## print(g)
## ---------------------------------------------------------
## test for convergence
## ---------------------------------------------------------
## [conv, flast, ipivot] <- cnvtst (alpha, pnorm, xnorm, ...
## difnew, ftest, gnorm, gtp, f, flast, g, ...
## ipivot, accrcy)
## cat("before cnvtst, flast=",flast,"\n")
ctres <- cnvtst (alpha, pnorm, xnorm,
difnew, ftest, gnorm, gtp, f, flast, g,
ipivot, accrcy)
conv <- ctres$conv
flast <- ctres$flast1
ipivot <- ctres$ipivot1
## cat("after cnvtst - conv, flast, ipivot:", conv, flast,"\n")
## print(ipivot)
if (conv) {
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror,
nfngr=ncg)
return(almqn)
}
g <- ztime (g, ipivot)
## ---------------------------------------------------------
## modify data for LMQN preconditioner
## ---------------------------------------------------------
if (! newcon) { ## 0 is FALSE and 1 is TRUE supposedly
envjn$yk <- g - oldg
envjn$sk <- alpha*p
envjn$yksk <- as.numeric(crossprod(envjn$yk, envjn$sk))
## cat("yksk:",envjn$yksk,"\n")
ireset <- ( (icycle == n-1) |
(difnew < epsred*(fkeep-f)) )
if (! ireset) {
envjn$yrsr <-
as.numeric(crossprod(envjn$yr,envjn$sr))
ireset <- (envjn$yrsr <= 0)
}
upd1 <- (envjn$yksk <= 0)
}
## cat("newcon, upd1:", newcon, upd1,"\n")
## ---------------------------------------------------------
## compute the search direction
## ---------------------------------------------------------
argvec <- c(accrcy, gnorm, xnorm)
## [p, gtp, ncg1, d] <- ...
## modlnp (d, x, g, maxit, upd1, ireset, bounds,
## ipivot, argvec, sfun)
mres <- modlnp (d, x, g, maxit, upd1, ireset,
bounds, ipivot, argvec, sfun, ...)
ncg1 <- mres$ncg1
gtp <- mres$gtp
d <- mres$dnew
p <- mres$p
## cat("New p:")
## print(p)
## tmp<-readline("cont.")
ncg <- ncg + ncg1
## ---------------------------------------------------------
## update LMQN preconditioner
## ---------------------------------------------------------
if (! newcon) {
if (ireset) {
envjn$sr <- envjn$sk
envjn$yr <- envjn$yk
fkeep <- f
icycle <- 1
} else {
envjn$sr <- envjn$sr + envjn$sk
envjn$yr <- envjn$yr + envjn$yk
icycle <- icycle + 1
}
}
} # end while !conv
}
lmqn <- function (x, sfun, maxit, maxfun, stepmx, accrcy, trace, ...) {
## ---------------------------------------------------------
## truncated-newton method for unconstrained minimization
## (customized version)
## ---------------------------------------------------------
# global vectors hyk sk yk sr yr & scalars yksk yrsr
## ---------------------------------------------------------
## set up
## ---------------------------------------------------------
## format compact
## format short e
n<-length(x)
# JN: Define globals here. Is it necessary to set up.
gtn<-list(yrsr=0, yksk=0, yr = rep(0, n), yk = rep(0, n), sr = rep(0, n), sk = rep(0, n))
envjn<<-list2env(gtn)
# end globals
eps <- .Machine$double.eps
upd1 <- 1
ncg <- 0
xnorm <- max(abs(x)) # norm(x,'inf')
ierror <- 0
if (stepmx < sqrt(accrcy) || maxfun < 1) {
ierror <- -1
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror,
nfngr=ncg)
return(almqn)
}
## ---------------------------------------------------------
## compute initial function value and related information
## ---------------------------------------------------------
fg <- sfun(x, ...)
#% print(fg)
g<- attr(fg, "gradient")
f<-fg
gnorm <- max(abs(g)) ## norm(g,'inf')
nf <- 1
nit <- 0
if (trace) cat("Itn ",nit," ",nf," ",ncg, " ",f, " ", gnorm,"\n")
## ---------------------------------------------------------
## check for small gradient at the starting point.
## ---------------------------------------------------------
ftest <- 1 + abs(f)
if (gnorm < .01*sqrt(eps)*ftest) {
ierror <- 0
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror, nfngr=ncg)
return(almqn)
}
## ---------------------------------------------------------
## set initial values to other parameters
## ---------------------------------------------------------
n <- length(x)
icycle <- n-1
toleps <- sqrt(accrcy) + sqrt(eps)
rtleps <- accrcy + eps
difnew <- 0
epsred <- .05
fkeep <- f
conv <- FALSE
ireset <- 0
ipivot <- 0
## ---------------------------------------------------------
## initialize diagonal preconditioner to the identity
## ---------------------------------------------------------
d <- rep(1,n) # as a vector
## ---------------------------------------------------------
## ..........main iterative loop..........
## ---------------------------------------------------------
## compute search direction
## ---------------------------------------------------------
argvec <- c(accrcy, gnorm, xnorm)
mres <- modlnp (d, x, g, maxit, upd1, ireset, bounds=FALSE, ipivot, argvec, sfun, ...)
p <- mres$p
# cat("p from first call to modlnp\n")
# print(p)
# tmp<-readline("cont.")
gtp <- mres$gtp
ncg1<-mres$ncg1
d <- mres$dnew
#% cat("d from modlnp:")
#% print(d)
ncg <- ncg + ncg1
while ( ! conv) {
#% cat("top while ")
#% tmp <- readline("Top of iteration")
oldg <- g
pnorm <- max(abs(p)) # norm(p,'inf')
oldf <- f
## ---------------------------------------------------------
## line search
## ---------------------------------------------------------
pe <- pnorm + eps
spe <- stepmx/pe
#% cat("gtp, spe:", gtp, spe,"\n")
alpha0 <- step1 (f, gtp, spe)
reslin <- lin1 (p, x, f, alpha0, g, sfun, ...)
## [x, f, g, nf1, ierror, alpha] <-
x <- reslin$xnew # need fixup
f <- reslin$fnew
g <- reslin$gnew
nf1 <- reslin$nf1
ierror <- reslin$ierror
alpha <- reslin$alpha1
# cat("after lin1, alpha=",alpha,"\n")
# tmp<-readline("cont.")
nf <- nf + nf1
## ---------------------------------------------------------
nit <- nit + 1
gnorm <- max(abs(g)) # norm(g,'inf')
### Display info
if (trace) cat("Itn ",nit," ",nf," ",ncg, " ",f, " ", gnorm,"\n")
if (ierror == 3) {
if (length(ncg) == 0) { ncg <- 0 } # ?? is.null(ncg)??
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror, nfngr=ncg)
return(almqn)
}
## ---------------------------------------------------------
## stop if more than maxfun evalutations have been made
## ---------------------------------------------------------
if (nf >= maxfun) {
ierror <- 2
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror, nfngr=ncg)
return(almqn)
}
## ---------------------------------------------------------
## set up for convergence and resetting tests
## ---------------------------------------------------------
ftest <- 1 + abs(f)
xnorm <- max(abs(x)) # norm(x,'inf')
difold <- difnew
difnew <- oldf - f
envjn$yk <- g - oldg
envjn$sk <- alpha*p
if (icycle == 1) {
if (difnew > 2*difold) { epsred <- 2*epsred }
if (difnew < 0.5*difold) { epsred <- 0.5*epsred }
}
## ---------------------------------------------------------
## convergence test
## ---------------------------------------------------------
conv <- ( ( (alpha*pnorm < toleps*(1 + xnorm)) &&
(abs(difnew) < rtleps*ftest) &&
(gnorm < accrcy^(1/3)*ftest) )||
( gnorm < .01*sqrt(accrcy)*ftest ) )
if (conv) {
ierror <- 0
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror, nfngr=ncg)
return(almqn)
}
## ---------------------------------------------------------
## update lmqn preconditioner
## ---------------------------------------------------------
envjn$yksk <- as.numeric(crossprod(envjn$yk, envjn$sk) )
ireset <- ((icycle == n-1) || (difnew < epsred*(fkeep-f)) )
if ( ! ireset) {
envjn$yrsr <- as.numeric(crossprod(envjn$yr, envjn$sr))
ireset <- (envjn$yrsr <= 0)
}
upd1 <- (envjn$yksk <= 0)
## ---------------------------------------------------------
## compute search direction
## ---------------------------------------------------------
argvec <- c(accrcy, gnorm, xnorm)
## cat("ireset =", ireset," upd1=",upd1,"\n")
#% cat("New d to modlnp:")
#% print(d)
## [p, gtp, ncg1, d] <- ...
mres <- modlnp (d, x, g, maxit, upd1, ireset, 0, ipivot, argvec, sfun, ...)
p <- mres$p
gtp <- mres$gtp
ncg1<-mres$ncg1
d <- mres$dnew
#% cat("returned dnew, envjn$john: ", envjn$john)
#% print(d)
ncg <- ncg + ncg1
## ---------------------------------------------------------
## store information for lmqn preconditioner
## ---------------------------------------------------------
if (ireset) {
envjn$sr <- envjn$sk
envjn$yr <- envjn$yk
fkeep <- f
icycle <- 1
} else {
envjn$sr <- envjn$sr + envjn$sk
envjn$yr <- envjn$yr + envjn$yk
icycle <- icycle + 1
}
} # end while
## [xstar, f, g, ierror] = ..
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror, nfngr=ncg)
} # end lmqn
modlnp <- function(d, x, g, maxit, upd1, ireset, bounds,
ipivot, argvec, sfun, ...) {
##---------------------------------------------------------
## this routine performs a preconditioned conjugate-gradient
## iteration to solve the Newton equations for a search
## direction for a truncated-newton algorithm.
## When the value of the quadratic model is sufficiently
## reduced, the iteration is terminated.
##---------------------------------------------------------
## parameters
##
## p - computed search direction
## g - current gradient
## gv,gz1,v - scratch vectors
## r - residual
## d - diagonal preconditoning matrix
## feval - value of quadratic function
##------------------------------------------------------------
## initialization
##------------------------------------------------------------
if (is.null(d)) stop("Null d")
## print(x)
## print(g)
## cat(bounds,"\n")
accrcy <- argvec[1]
gnorm <- argvec[2]
xnorm <- argvec[3]
if (maxit == 0) {
p <- -g
gtp <- as.numeric(crossprod(p, g))
ncg1 <- 1
dnew <- d
if (sqrt(sum(p^2))==0) {
## cat("MODLNP 01: |p| = 0\n")
## pause(1)
}
#% cat("modlnp - dout01:")
#% print(dnew)
result<-list(p=p, gtp=gtp, ncg1=ncg1, dnew=dnew)
return(result)
}
first <- 1
tol <- 1e-6 ######### was 1.d-12 #########
qold <- 0
ainit <- initpc (d, upd1, ireset)
dnew<-ainit$td
#% cat("after initpc dnew:")
#% print(dnew)
r <- -g
v <- rep(0, length(r))
p <- v
gtp <- as.numeric(crossprod(p, g))
rho <- rep(0, maxit+1)
beta <- rep(0, maxit)
v.gv <- rep(0, maxit)
rho[[1]] <- as.numeric(crossprod(r))
##------------------------------------------------------------
## main iteration (conjugate-gradient iterations for Ax = b)
##------------------------------------------------------------
ind <- 0
ncg1 <- 0
for (k in 1:maxit) {
ncg1 <- ncg1 + 1
if (bounds) { r <- ztime(r, ipivot) }
amsolve <- msolve (r, upd1, ireset, first, d)
zk<-amsolve$y
if (bounds) { zk <- ztime (zk, ipivot) }
#% cat("r:")
#% print(r)
#% cat("zk:")
#% print(zk)
rz <- as.numeric(crossprod(r,zk))
if (rz/gnorm < tol) {
ind <- 80
if (sqrt(sum(p^2))==0) {
p <- -g
gtp <- as.numeric(crossprod(p, g))
}
#% cat("modlnp - dout - ind80:")
#% print(dnew)
result<-list(p=p, gtp=gtp, ncg1=ncg1, dnew=dnew)
return(result)
}
if (k > 1) {
beta[k] <- rz/rzold
} else {
beta[k] <- 0
}
#% cat("beta[",k,"] =",beta[k],"\n")
v <- zk + beta[k]*v
if (bounds) { v <- ztime( v, ipivot) }
## cat("about to call gtims\n")
gv <- gtims(v, x, g, accrcy, xnorm, sfun, ...)
#% cat("After gtims: gv=")
#% print(gv)
## cat(bounds,"\n")
if (bounds) { gv <- ztime (gv, ipivot) }
#% cat("gv, v:")
## cat("gv=", gv, "\n")
#% print(gv)
#% print(v)
v.gv[[k]] <- as.numeric(crossprod(v, gv)) ## ?? v.gv has underscore!! ??
#% cat("v.gv[[",k,"]]=",v.gv[[k]],"\n")
if (v.gv[[k]]/gnorm < tol) {
ind <- 50
if (sqrt(sum(p^2))==0) {
p <- -g ## Mod SGN 140912
if (bounds) {
p <- ztime(p, ipivot)
}
gtp <- crossprod(p, g)
## cat("MODLNP 03: |p| = 0 \n")
## pause(1)
}
#% cat("modlnp - dout03:")
#% print(dnew)
result<-list(p=p, gtp=gtp, ncg1=ncg1, dnew=dnew)
return(result)
}
dnew <- ndia3(dnew, v, gv, r, v.gv[[k]]) # NB need to return something but it doesn't get used
#% cat("after ndia3 below dout03, dnew:")
#% print(dnew)
##------------------------------------------------------------
## compute current solution and related vectors
##------------------------------------------------------------
alpha <- rz / v.gv[[k]]
p <- p + alpha* v
r <- r - alpha*gv
rho[[k+1]] <- as.numeric(crossprod(r))
##------------------------------------------------------------
## test for convergence
##------------------------------------------------------------
gtp <- as.numeric(crossprod(p,g))
pr <- as.numeric(crossprod(r,p))
q <- (gtp + pr) / 2
qtest <- k * (1 - qold/q)
if (qtest <= 0.5) {
if (sqrt(sum(p^2))==0) {
## cat("MODLNP 04: |p| = 0\n")
## pause(1)
}
#% cat("modlnp - dout04:")
#% print(dnew)
result<-list(p=p, gtp=gtp, ncg1=ncg1, dnew=dnew)
return(result)
}
qold <- q
##------------------------------------------------------------
## perform cautionary test
##------------------------------------------------------------
if (gtp > 0) {
ind <- 40
if (sqrt(sum(p^2))==0) {
## cat("MODLNP 05: |p| = 0 \n")
## pause(1)
}
#% cat("modlnp - dout05:")
#% print(dnew)
result<-list(p=p, gtp=gtp, ncg1=ncg1, dnew=dnew)
return(result)
}
rzold <- rz
} ## end loop over k
k <- k-1
##------------------------------------------------------------
## terminate algorithm
##------------------------------------------------------------
if (ind == 40) {
p <- p - alpha*v
}
##------------------------------------------------------------
if (ind == 50 && k <= 1) {
amsolve <- msolve (g, upd1, ireset, first, d)
p<-amsolve$y
p <- -p
if (bounds) { p <- ztime (p, ipivot) }
}
##------------------------------------------------------------
if (ind == 80 && k <= 1) {
p <- -g
if (bounds) { p <- ztime (p, ipivot) }
}
##------------------------------------------------------------
## store new diagonal preconditioner
##------------------------------------------------------------
gtp <- as.numeric(crossprod(p, g))
ncg1 <- k + 1
if (sqrt(sum(p^2))==0) {
## cat("MODLNP 06: |p| = 0 \n")
## pause(1)
}
## cat("modlnp - dout06:")
#% print(dnew)
result<-list(p=p, gtp=gtp, ncg1=ncg1, dnew=dnew)
return(result)
}
modz <- function (x, p, ipivot, low, up, flast, f, alpha){
##---------------------------------------------------------------------
## update the constraint matrix if a new constraint is encountered
##---------------------------------------------------------------------
eps <- .Machine$double.eps
indl <- which(ipivot == 0 & p < 0)
if (length(indl) > 0) {
toll <- 10 * eps * (abs(low[indl]) + 1)
hitl <- which(x[indl]-low[indl] <= toll)
if (length(hitl) > 0) {
flast <- f
ipivot[indl[hitl]] <- -1
}
}
##---------------------------------------------------------------------
indu <- which((ipivot == 0) & (p > 0));
if (length(indu) > 0) {
tolu <- 10 * eps * (abs( up[indu]) + 1)
hitu <- which(up[indu]-x[indu] <= tolu)
if (length(hitu) > 0) {
flast <- f
ipivot[indu[hitu]] <- 1
}
}
##---------------------------------------------------------------------
flast1 <- flast
ipivot1 <- ipivot
list(flast1 = flast, ipivot1=ipivot)
}
msolve <- function(g, upd1, ireset, first, d) {
##---------------------------------------------------------
## This routine acts as a preconditioning step for the
## linear conjugate-gradient routine. It is also the
## method of computing the search direction from the
## gradient for the non-linear conjugate-gradient code.
## It represents a two-step self-scaled bfgs formula.
##---------------------------------------------------------
#% cat("msolve, envjn$john=",envjn$john,"\n")
envjn$john<-"msolve"
if (upd1) {
#% cat("upd1 is TRUE\n")
y <- g/d
} else {
gsk <- as.numeric(crossprod(g, envjn$sk))
if (ireset) {
envjn$hg <- g/d
if (first) {
envjn$hyk <- envjn$yk/d
ykhyk <- as.numeric(crossprod(envjn$yk, envjn$hyk))
}
ghyk <- as.numeric(crossprod(g, envjn$hyk))
y <- ssbfgs(envjn$sk,envjn$hg,envjn$hyk,envjn$yksk,ykhyk,gsk,ghyk)
} else {
envjn$hg <- g/d
if (first) {
envjn$hyk <- envjn$yk/d
envjn$hyr <- envjn$yr/d
envjn$yksr <- as.numeric(crossprod(envjn$yk, envjn$sr))
ykhyr <- as.numeric(crossprod(envjn$yk, envjn$hyr))
}
gsr <- as.numeric(crossprod(g, envjn$sr))
ghyr <- as.numeric(crossprod(g, envjn$hyr))
if (first) {
yrhyr <- as.numeric(crossprod(envjn$yr, envjn$hyr))
}
envjn$hg <- ssbfgs(envjn$sr,envjn$hg,envjn$hyr,envjn$yrsr,yrhyr,gsr,ghyr)
if (first) {
envjn$hyk <- ssbfgs(envjn$sr,envjn$hyk,envjn$hyr,envjn$yrsr,yrhyr,envjn$yksr,ykhyr)
}
ykhyk <- as.numeric(crossprod(envjn$hyk, envjn$yk))
ghyk <- as.numeric(crossprod(envjn$hyk, g))
y <- ssbfgs(envjn$sk,envjn$hg,envjn$hyk,envjn$yksk,ykhyk,gsk,ghyk)
}
}
amsolve<-list(y=y) ## returns y
}
ndia3<-function(e, v, gv, r, vgv){
##---------------------------------------------------------
## update the preconditioning matrix based on a diagonal
## version of the bfgs quasi-newton update.
##---------------------------------------------------------
tol <- 1e-6
vr <- as.numeric(crossprod(v,r))
#% cat("ndia3: vgv, vr, then e:",vgv, vr,"\n")
#% print(e)
if (abs(vr)>tol && abs(vgv)>tol) {
e <- e - (r*r)/vr + (gv*gv)/vgv # CAUTION! May be crossprod
ind <- which(e < tol)
if (length(ind)>0){
e[ind] <- 1
}
}
#% print(e)
e # Need to return the object?? Possibly not!
}
norm2 <- function (x) {
##---------------------------------------------------------
## compute the 2 norm of vector x
##---------------------------------------------------------
res<-sqrt(as.numeric(crossprod(x)))
}
ssbfgs <- function (s, hv, hy, ys, yhy, vs, vhy){
##---------------------------------------------------------
## self-scaled bfgs quasi-Newton update
## (used by preconditioner)
##---------------------------------------------------------
delta <- (1 + yhy/ys)*vs/ys - vhy/ys
beta <- -vs/ys
z <- hv + delta*s + beta*hy
}
step1 <- function(f, gtp, smax) {
##---------------------------------------------------------
## step1 returns the length of the initial step to be
## taken along the vector p in the next linear search.
##---------------------------------------------------------
## [fm is supposed to be an estimate of the optimal function value]
eps<-.Machine$double.eps
fm <- 0
d <- abs(f-fm)
alpha <- min(1, smax)
if ((2*d <= (-gtp)) && (d >= eps)) {
alpha = min(-2*d/gtp, smax)
}
alpha # need to ensure returned value
}
stpmax <- function(stepmx, pe, x, p, ipivot, low, up) {
##------------------------------------------------
## compute the maximum allowable step length
## (spe is the standard (unconstrained) max step)
##------------------------------------------------
## cat("stpmax: stepmx, pe:",stepmx, pe,"\n")
spe <- stepmx / pe
al <- spe
au <- spe
##------------------------------------------------
indl <- which(ipivot==0 & p < 0)
if ( length(indl)>0) {
tl <- low[indl] - x[indl]
al <- min(tl/p[indl])
}
##------------------------------------------------
indu <- which(ipivot==0 & p > 0)
if (length(indu)>0) {
tu <- up[indu] - x[indu]
au <- min(tu/p[indu])
}
##------------------------------------------------
spe <- min(c(spe, al, au))
}
ztime<-function(x,ipivot) {
## ---------------------------------------------------------
## this routine multiplies the vector x by the
## constraint matrix z
## ---------------------------------------------------------
ind <- which(ipivot!=0);
x[ind] <- 0
x1 <- x;
}
Try the Rtnmin package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.