Nothing
R/squarem.R
In SQUAREM: Squared Extrapolation Methods for Accelerating EM-Like Monotone Algorithms
Defines functions
fpiter
cyclem2
cyclem1
squarem2
squarem1
squarem
Documented in
cyclem1
cyclem2
fpiter
squarem
squarem1
squarem2
squarem <- function(par, fixptfn, objfn, ... , control=list()) {
#
# The Master function for SQUAREM acceleration of "ANY" fixed-point iteration
# A wrapper that calls appropriate function below (squarem1, squarem2, cyclem1, cyclem2)
#
# Author: Ravi Varadhan, Johns Hopkins University
#
# See Varadhan and Roland (Scandinavian J Statistics 2008)
# See Roland, Varadhan and Frangakis (Applied Numerical Mathematics 2007)
#
# Last modified: June 25, 2010
# Last modified: August 09, 2010
# Last modified: December 31, 2010
#######################################################################
#
# par = starting value of parameter vector
# fixptfn = fixed-point iteration F(x)
# for which the solution: F(x*) = x* is sought
# objfn = underlying objective function which is minimized at x*
# method = 1, 2, or 3, indicating the type of steplength to be used
#
#
control.default <- list(K = 1, method=3, minimize=TRUE, square=TRUE, step.min0=1, step.max0=1, mstep=4, kr=1, objfn.inc=1,
tol=1.e-07, maxiter=1500, trace=FALSE, intermed=FALSE)
namc <- names(control)
if (!all(namc %in% names(control.default)))
stop("unknown names in control: ", namc[!(namc %in% names(control.default))])
ctrl <- modifyList(control.default, control)
if (ctrl$K > 1 & !(ctrl$method %in% c("rre", "mpe")) ) ctrl$method <- "rre"
if (ctrl$K == 1 & !(ctrl$method %in% c(1,2,3)) ) ctrl$method <- 3
if (!missing(objfn) ) {
if (ctrl$K == 1) sqobj <- squarem1(par, fixptfn, objfn, ... , control=ctrl)
else if (ctrl$K > 1 | ctrl$method %in% c("rre", "mpe")) sqobj <- cyclem1(par, fixptfn, objfn, ... , control=ctrl)
} else {
if (ctrl$K == 1) sqobj <- squarem2(par, fixptfn, ... , control=ctrl)
else if (ctrl$K > 1 | ctrl$method %in% c("rre", "mpe")) sqobj <- cyclem2(par, fixptfn, ... , control=ctrl)
}
return(sqobj)
}
#######################################################################
# Partially-monotone, globally-convergent SQUAREM acceleration
# Can accelerate EM, MM, and other fixed-point iterations
# Here we use an adaptive maximum step-size strategy
# This strategy is effective and represents a good trade-off between speed and stability
# See Varadhan and Roland (Scandinavian J Statistics 2008)
#
# Date: August 22, 2007
# Adaptive step-max strategy modified: May 2008
# Last modified: November 2008
# Last modified: October 2009
# Last modified: Feb 23 2010
#######################################################################
squarem1 <- function(par, fixptfn, objfn, ... , control=list()) {
# par = starting value of parameter vector
# fixptfn = fixed-point iteration F(x)
# for which the solution: F(x*) = x* is sought
# objfn = underlying objective function which is minimized at x*
#
#
control.default <- list(K=1, square=TRUE, minimize=TRUE, method=3, step.min0=1, step.max0=1, mstep=4, kr=1, objfn.inc=1,
tol=1.e-07, maxiter=1500, trace=FALSE, intermed=FALSE)
namc <- names(control)
if (!all(namc %in% names(control.default)))
stop("unknown names in control: ", namc[!(namc %in% names(control.default))])
ctrl <- modifyList(control.default, control)
#####
# method = 1, 2, or 3, indicating the type of steplength to be used
# A key parameter is `step.min0'. This can be either positive or negative if the eigenvalues of the Jacobian of fixed-point mapping are all positive;
# We set it to +1 for contraction mappings such as EM and MM algorithms
# Must be negative, e.g. -1, if the fixed-point mapping can have negative eiegnvalues
#####
# parameter "objfn.inc" dictates the amount of non-monotonicity in the objective function
# setting objfn.inc=0 would enforce monotonicity, whereas objfn.inc=Inf would be a non-monotonic scheme
# The defalut objfn.inc=1 would enforce monotonicity far from solution, but allows for non-monotonicity closer to solution
#
method <- ctrl$method
maxiter <- ctrl$maxiter
minimize <- ctrl$minimize
tol <- ctrl$tol
step.min <- ctrl$step.min0
step.max0 <- ctrl$step.max0
step.max <- ctrl$step.max0
mstep <- ctrl$mstep
objfn.inc <- ctrl$objfn.inc
trace <- ctrl$trace
intermed <- ctrl$intermed
if (trace) cat("Squarem-1 \n ")
if (missing(objfn)) stop("\n squarem2 should be used if objective function is not available \n\n")
iter <- 1
objval <- rep(NA,1)
p <- par
lold <- objfn(p, ...)
leval <- 1
if (trace) cat("Objective fn: ", lold, "\n")
feval <- 0
conv <- TRUE
p.inter <- c(p, lold)
while (feval < maxiter) {
extrap <- TRUE
p1 <- try(fixptfn(p, ...),silent=TRUE)
feval <- feval + 1
if (inherits(p1, "try-error") | any(is.nan(unlist(p1)))) stop("Error in function evaluation")
q1 <- p1 - p
sr2 <- crossprod(q1)
if (sqrt(sr2) < tol) {p <- p1; break}
p2 <- try(fixptfn(p1, ...),silent=TRUE)
feval <- feval + 1
if (inherits(p2, "try-error") | any(is.nan(unlist(p2)))) stop("Error in function evaluation")
q2 <- p2 - p1
sq2 <- sqrt(crossprod(q2))
if (sq2 < tol) {p <- p2; break}
sv2 <- crossprod(q2-q1)
srv <- crossprod(q1, q2-q1)
alpha <- switch(method, -srv/sv2, -sr2/srv, sqrt(sr2/sv2))
alpha <- max(step.min, min(step.max, alpha))
p.new <- p + 2*alpha*q1 + alpha^2*(q2-q1)
if (abs(alpha - 1) > 0.01 ) {
p.new <- try(fixptfn(p.new , ...),silent=TRUE) # stabilization step
feval <- feval + 1
}
if (inherits(p.new, "try-error") | any(is.nan(p.new))) {
p.new <- p2
lnew <- try(objfn(p2, ...), silent=TRUE)
leval <- leval + 1
if (alpha == step.max) step.max <- max(step.max0, step.max/mstep)
alpha <- 1
extrap <- FALSE
} else {
if (is.finite(objfn.inc)) {
lnew <- try(objfn(p.new, ...), silent=TRUE)
leval <- leval + 1
} else lnew <- lold
non.monotone <- if (minimize) lnew > lold + objfn.inc else lnew < lold - objfn.inc
if (inherits(lnew, "try-error") | is.nan(lnew) | non.monotone) {
p.new <- p2
lnew <- try(objfn(p2, ...), silent=TRUE)
leval <- leval + 1
if (alpha==step.max) step.max <- max(step.max0, step.max/mstep)
alpha <- 1
extrap <- FALSE
}
}
if (alpha == step.max) step.max <- mstep*step.max
if (step.min < 0 & alpha == step.min) step.min <- mstep*step.min
p <- p.new
if (!is.nan(lnew)) lold <- lnew
if (trace) cat("Objective fn: ", lnew, " Extrapolation: ", extrap, " Steplength: ", alpha, "\n")
if(intermed) p.inter <- rbind(p.inter, c(p, lnew))
iter <- iter+1
}
if (feval >= maxiter) conv <- FALSE
if (is.infinite(objfn.inc)) {
lold <- objfn(p, ...)
leval <- leval + 1
}
rownames(p.inter) <- NULL
if (!intermed) return(list(par=p, value.objfn=lold, iter= iter, fpevals=feval, objfevals = leval, convergence=conv))
else return(list(par=p, value.objfn=lold, iter= iter, fpevals=feval, objfevals = leval, convergence=conv, p.intermed=p.inter))
}
###################################################################
squarem2 <- function(par, fixptfn, ... , control=list() ) {
# par = starting value of parameter vector
# fixptfn = fixed-point iteration F(x)
# for which the solution: F(x*) = x* is sought
# method = 1, 2, or 3, indicating the type of steplength to be used
#
# See Varadhan and Roland (Scandinavian J Statistics 2008)
#
control.default <- list(K=1, square=TRUE, minimize=TRUE, method=3, step.min0=1, step.max0=1, mstep=4, kr=1, objfn.inc=1,
tol=1.e-07, maxiter=1500, trace=FALSE, intermed=FALSE)
namc <- names(control)
if (!all(namc %in% names(control.default)))
stop("unknown names in control: ", namc[!(namc %in% names(control.default))])
ctrl <- modifyList(control.default, control)
#####
# method = 1, 2, or 3, indicating the type of steplength to be used
# A key parameter is `step.min0'. This can be either positive or negative if the eigenvalues of the Jacobian of fixed-point mapping are all positive;
# We set it to +1 for contraction mappings such as EM and MM algorithms
# Must be negative, e.g. -1, if the fixed-point mapping can have negative eiegnvalues
#####
# Parameter "kr" dictates the amount of non-monotonicity in the fixed-point residual
# setting kr=0 would try to enforce monotonicity; kres=Inf would be a non-monotonic scheme
#
method <- ctrl$method
maxiter <- ctrl$maxiter
tol <- ctrl$tol
kr <- ctrl$kr
step.min <- ctrl$step.min0
step.max0 <- ctrl$step.max0
step.max <- ctrl$step.max0
mstep <- ctrl$mstep
trace <- ctrl$trace
if (trace) cat("Squarem-2 \n")
intermed <- ctrl$intermed
iter <- 1
feval <- 0
kount <- 0
conv <- TRUE
p.inter <- rep(NA, length(par)+1)
while (feval < maxiter) {
extrap <- TRUE
p1 <- try(fixptfn(par, ...),silent=TRUE)
feval <- feval + 1
if (inherits(p1, "try-error") | any(is.nan(unlist(p1)))) break
q1 <- p1 - par
sr2 <- crossprod(q1)
p.inter <- rbind(p.inter, c(par, sqrt(sr2)))
if (sqrt(sr2) < tol) {par <- p1; break}
p2 <- try(fixptfn(p1, ...),silent=TRUE)
feval <- feval + 1
if (inherits(p2, "try-error") | any(is.nan(unlist(p2)))) break
q2 <- p2 - p1
sq2 <- sqrt(crossprod(q2))
res <- sq2
if (sq2 < tol) {par <- p2; break}
sv2 <- crossprod(q2-q1)
srv <- crossprod(q1, q2-q1)
alpha <- switch(method, -srv/sv2, -sr2/srv, sqrt(sr2/sv2))
alpha <- max(step.min, min(step.max, alpha))
p.new <- par + 2*alpha*q1 + alpha^2*(q2-q1)
if (abs(alpha - 1) > 0.01 ) {
ptmp <- try(fixptfn(p.new, ...),silent=TRUE) # stabilization step
feval <- feval + 1
if (inherits(ptmp, "try-error") | any(is.nan(unlist(ptmp))) ) {
p.new <- p2
if (alpha == step.max) step.max <- max(step.max0, step.max/mstep)
alpha <- 1
extrap <- FALSE
} else {
res <- sqrt(crossprod(ptmp - p.new))
parnorm <- sqrt(crossprod(p2)/length(p2))
kres <- kr * (1 + parnorm) + sq2
p.new <- if (res <= kres) ptmp else p2
if (res > kres) {
if (alpha == step.max) step.max <- max(step.max0, step.max/mstep)
alpha <- 1
extrap <- FALSE
}
}
}
if (alpha == step.max) step.max <- mstep*step.max
if (step.min < 0 & alpha == step.min) step.min <- mstep*step.min
if (trace) cat("Residual: ", res, " Extrapolation: ", extrap, " Steplength: ", alpha, "\n")
par <- p.new
iter <- iter+1
}
p.inter <- p.inter[-1, ]
if (feval >= maxiter) conv <- FALSE
if (intermed) return(list(par=par, value.objfn=NA, iter = iter, fpevals=feval, objfevals = 0,
convergence=conv, p.intermed=p.inter))
else return(list(par=par, value.objfn=NA, iter = iter, fpevals=feval, objfevals = 0, convergence=conv))
}
#######################################################################
# Squared-Cycled vector extrapolation for accelerating vector sequences
# See Roland, Varadhan and Frangakis (Applie Numerical Mathematics 2007)
# Author: Ravi Varadhan, Johns Hopkins University
# Latest modofications:: Feb 23, 2010
#
###########################
cyclem1 <- function(par, fixptfn, objfn, control=list(), ...) {
#
# par = starting value of parameter vector
# fixptfn = fixed-point iteration F(x)
# for which the solution: F(x*) = x* is sought
# objfn = underlying objective function which is minimized at x*
#
control.default <- list(K=2, method="rre", minimize=TRUE, square=TRUE, step.min0=1, step.max0=1, mstep=4, kr=1, objfn.inc=1,
tol=1.e-07, maxiter=1500, trace=FALSE, intermed=FALSE)
namc <- names(control)
if (!all(namc %in% names(control.default)))
stop("unknown names in control: ", namc[!(namc %in% names(control.default))])
ctrl <- modifyList(control.default, control)
#
# method = reduced-rank or minimal-polynomial extrapolation
# K = order of extrapolation scheme; K=2,3,4 are typical choices
# square = a logical variable indicating whether or not "squaring" is used
method <- ctrl$method
K <- ctrl$K
square <- ctrl$square
tol <- ctrl$tol
maxiter <- ctrl$maxiter
objfn.inc <- ctrl$objfn.inc
trace <- ctrl$trace
if (trace) cat("Cyclem-1 \n")
if (missing(objfn)) stop("\n cyclem2 should be used if objective function is not available \n\n")
if ( K > length(par)) {
cat("K is too large. Decrease it. \n")
return()
}
iter <- 0
lold <- objfn(par, ...)
leval <- 1
ptmp <- fixptfn(par,...)
res <- sqrt(c(crossprod(ptmp - par)))
par <- ptmp
feval <- 1
conv <- TRUE
while (feval < maxiter & res > tol) {
iter <- iter+1
extrap <- TRUE
p <- matrix(NA, nrow=length(par), ncol=K+2)
U <- matrix(NA, nrow=length(par), ncol=K+1)
p[,1] <- par
for (i in 2:(K+2)){
p[,i] <- try(fixptfn(p[,i-1], ...), silent=TRUE)
if (inherits(p[,i], "try-error") | any(is.nan(unlist(p[,i]))) | any(!is.numeric(p[,i]))) stop("Error in function evaluation \n")
U[,i-1] <- p[, i] - p[, i-1]
}
feval <- feval + K + 1
sqKp1 <- sqrt(c(crossprod(U[, K+1])))
if (method == "rre"){
coef <- try(solve(qr(crossprod(U), LAPACK=TRUE, tol=1.e-14), rep(1,K+1)), silent=TRUE)
if (inherits(coef, "try-error") | any(is.nan(coef)) ) {
extrap <- FALSE
} else {
if (abs(sum(coef)) < 1.e-07) extrap <- FALSE
coef <- coef/sum(coef)
}
}
if (method == "mpe"){
coef <- try(solve(qr(U[,-(K+1)], LAPACK=TRUE, tol=1.e-14), -U[,K+1]), silent=TRUE)
if (inherits(coef, "try-error") | any(is.nan(coef))) {
extrap <- FALSE
} else {
coef <- c(coef, 1)
if (abs(sum(coef)) < 1.e-07) extrap <- FALSE
coef <- coef/sum(coef)
}
}
if (!extrap) {
par <- p[, ncol(p)]
res <- sqKp1
iter <- iter + 1
if (trace) cat(" objective fn: ", lold, "extrapolation: ", extrap, "\n")
next
}
if (!square) pnew <- c(p[,-(K+2)] %*% coef)
if (square) {
pnew <- rep(0, length(par))
if (K > 1) {
for (i in 1:(K-1)) {
p <- try(cbind(p, fixptfn(p[,i+K+1],...)), silent=TRUE)
if (inherits(p, "try-error")| any(is.nan(unlist(p))) | any(!is.numeric(p))) stop("Error in function evaluation \n")
}
feval <- feval + K - 1
}
for (i in 0:K){
for (j in 0:K){
pnew <- pnew + coef[i+1]* coef[j+1] * p[,i+j+1]
} }
sqKp1 <- sqrt(c(crossprod(p[, 2*K] - p[, 2*K-1])))
}
# Perform one stabilization step between cycles
if (extrap) {
ptmp <- try(fixptfn(pnew, ...), silent=TRUE)
res <- sqrt(c(crossprod(ptmp - pnew)))
feval <- feval + 1
}
if (inherits(ptmp, "try-error") | any(is.nan(unlist(ptmp))) ) {
pnew <- p[, ncol(p)]
res <- sqKp1
extrap <- FALSE
} else pnew <- ptmp
lnew <- try(objfn(pnew, ...), silent=TRUE)
leval <- leval + 1
non.monotone <- if(ctrl$minimize) lnew > lold+objfn.inc else lnew < lold-objfn.inc
if (inherits(lnew, "try-error") | is.nan(lnew) | non.monotone) {
pnew <- p[, ncol(p)]
res <- sqKp1
extrap <- FALSE
if (trace) cat(" objective fn: ", lold, "extrapolation: ", extrap, "\n")
} else lold <- lnew
par <- pnew
if (trace) cat(" objective fn: ", lold, "extrapolation: ", extrap, "\n")
} # Main loop complete
if (feval >= maxiter) conv <- FALSE
return(list(par=par, value.objfn=lold, iter = iter, fpevals=feval, objfevals = leval, convergence=conv))
}
###################################################################
# Squared-Cycled vector extrapolation for accelerating vector sequences
# See Roland, Varadhan and Frangakis (Applie Numerical Mathematics 2007)
# Author: Ravi Varadhan, Johns Hopkins University
# Latest modofications:: Feb 23, 2010
#
###########################
cyclem2 <- function(par, fixptfn, control=list(), ...) {
# par = starting value of parameter vector
# fixptfn = fixed-point iteration F(x)
# for which the solution: F(x*) = x* is sought
# method = reduced-rank or minimal-polynomial extrapolation
#
control.default <- list(K=2, method="rre", minimize=TRUE, square=TRUE, step.min0=1, step.max0=1, mstep=4, kr=1, objfn.inc=1,
tol=1.e-07, maxiter=1500, trace=FALSE, intermed=FALSE)
namc <- names(control)
if (!all(namc %in% names(control.default)))
stop("unknown names in control: ", namc[!(namc %in% names(control.default))])
ctrl <- modifyList(control.default, control)
#
# method = reduced-rank ("rre") or minimal-polynomial ("mpe") extrapolation
# K = order of extrapolation scheme; K=2,3,4 are typical choices
# square = a logical variable indicating whether or not "squaring" is used
method <- ctrl$method
K <- ctrl$K
square <- ctrl$square
tol <- ctrl$tol
maxiter <- ctrl$maxiter
kr <- ctrl$kr
trace <- ctrl$trace
if (trace) cat("Cyclem-2 \n")
if ( K > length(par)) {
cat("K is too large. Decrease it. \n")
return()
}
iter <- 1
ptmp <- fixptfn(par,...)
res <- sqrt(c(crossprod(ptmp - par)))
par <- ptmp
conv <- TRUE
feval <- 1
while (feval < maxiter & res > tol) {
extrap <- TRUE
p <- matrix(NA, nrow=length(par), ncol=K+2)
U <- matrix(NA, nrow=length(par), ncol=K+1)
p[,1] <- par
for (i in 2:(K+2)){
p[,i] <- try(fixptfn(p[,i-1], ...), silent=TRUE)
if (inherits(p[,i], "try-error") | any(is.nan(unlist(p[,i]))) | any(!is.numeric(p[,i]))) stop("Error in function evaluation")
U[,i-1] <- p[, i] - p[, i-1]
}
feval <- feval + K + 1
sqKp1 <- sqrt(c(crossprod(U[, K+1])))
if (method == "rre"){
coef <- try(solve(qr(crossprod(U), LAPACK=TRUE, tol=1.e-14), rep(1,K+1)), silent=TRUE)
if (inherits(coef, "try-error") | any(is.nan(coef))) {
extrap <- FALSE
} else {
if (abs(sum(coef)) < 1.e-07) extrap <- FALSE
coef <- coef/sum(coef)
}
}
if (method == "mpe"){
coef <- try(solve(qr(U[,-(K+1)], LAPACK=TRUE, tol=1.e-14), -U[,K+1]), silent=TRUE)
if (inherits(coef, "try-error") | any(is.nan(coef))) {
extrap <- FALSE
} else {
coef <- c(coef, 1)
if (abs(sum(coef)) < 1.e-07) extrap <- FALSE
coef <- coef/sum(coef)
}
}
if (!extrap) {
par <- p[, ncol(p)]
res <- sqKp1
iter <- iter + 1
if (trace) cat(" residual: ", res, "extrapolation: ", extrap, "\n")
next
}
if (!square) pnew <- c(p[,-(K+2)] %*% coef)
if (square) {
pnew <- rep(0, length(par))
if (K > 1) {
for (i in 1:(K-1)) {
p <- try(cbind(p, fixptfn(p[,i+K+1],...)), silent=TRUE)
if (inherits(p, "try-error") | any(is.nan(unlist(p))) | any(!is.numeric(p))) stop("Error in function evaluation")
}
feval <- feval + K - 1
}
for (i in 0:K){
for (j in 0:K){
pnew <- pnew + coef[i+1]* coef[j+1] * p[,i+j+1]
} }
sqKp1 <- sqrt(c(crossprod(p[, 2*K] - p[, 2*K-1])))
}
# Perform one stabilization step between cycles
ptemp <- try(fixptfn(pnew, ...), silent=TRUE)
feval <- feval + 1
if(inherits(ptemp, "try-error") | any(is.nan(ptemp)) ){
par <- p[, ncol(p)]
res <- sqKp1
extrap <- FALSE
} else {
res <- sqrt(c(crossprod(ptemp - pnew)))
parnorm <- as.numeric(sqrt(crossprod(p[, ncol(p)])/length(par)))
if (res <= (kr * (1 + parnorm) + sqKp1)) {
par <- ptemp
extrap <- TRUE
} else {
par <- p[, ncol(p)]
res <- sqKp1
extrap <- FALSE
}
}
iter <- iter+1
if (trace) cat(" residual: ", res, "extrapolation: ", extrap, "\n")
} # Main loop complete
if (feval >= maxiter) conv <- FALSE
return(list(par=par, value.objfn=NA, iter = iter, fpevals=feval, objfevals = 0, convergence=conv))
}
###################################################################
# The EM algorithm
#
fpiter <- function(par, fixptfn, objfn=NULL, control=list( ), ...){
control.default <- list(tol=1.e-07, maxiter=5000, trace=FALSE, intermed=FALSE)
namc <- names(control)
if (!all(namc %in% names(control.default)))
stop("unknown names in control: ", namc[!(namc %in% names(control.default))])
ctrl <- modifyList(control.default, control)
#
# method = reduced-rank ("rre") or minimal-polynomial ("mpe") extrapolation
# K = order of extrapolation scheme; K=2,3,4 are typical choices
# square = a logical variable indicating whether or not "squaring" is used
tol <- ctrl$tol
maxiter <- ctrl$maxiter
trace <- ctrl$trace
intermed <- ctrl$intermed
if (trace) cat("fpiter \n")
iter <- 1
resid <- rep(NA,1)
objeval <- 0
conv <- FALSE
if(intermed) p.inter <- rep(NA, 1+length(par))
while (iter < maxiter) {
p.new <- fixptfn(par, ...)
res <- sqrt(crossprod(p.new - par))
if(intermed) p.inter <- rbind(p.inter, c(par, res))
if ( res < tol) {conv <- TRUE; break}
if (trace) {
if (!is.null(objfn)) {cat("Iter: ", iter, "Objective fn: ",objfn(par, ...), "\n"); objeval <- objeval + 1}
else cat("Iter: ", iter, "Residual: ",res, "\n")
}
par <- p.new
iter <- iter+1
}
if(intermed) p.inter <- p.inter[-1,]
loglik.best <- if (!is.null(objfn)) objfn(par, ...) else NA
if(!intermed) return(list(par=par, value.objfn=loglik.best, fpevals=iter, objfevals = objeval,
convergence=conv))
else return(list(par=par, value.objfn=loglik.best, fpevals=iter, objfevals = objeval,
p.intermed=p.inter, convergence=conv))
}
###################################################################
Try the SQUAREM package in your browser
Any scripts or data that you put into this service are public.
SQUAREM documentation built on Jan. 16, 2021, 5:38 p.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 fpiter cyclem2 cyclem1 squarem2 squarem1 squarem
Documented in cyclem1 cyclem2 fpiter squarem squarem1 squarem2
squarem <- function(par, fixptfn, objfn, ... , control=list()) {
#
# The Master function for SQUAREM acceleration of "ANY" fixed-point iteration
# A wrapper that calls appropriate function below (squarem1, squarem2, cyclem1, cyclem2)
#
# Author: Ravi Varadhan, Johns Hopkins University
#
# See Varadhan and Roland (Scandinavian J Statistics 2008)
# See Roland, Varadhan and Frangakis (Applied Numerical Mathematics 2007)
#
# Last modified: June 25, 2010
# Last modified: August 09, 2010
# Last modified: December 31, 2010
#######################################################################
#
# par = starting value of parameter vector
# fixptfn = fixed-point iteration F(x)
# for which the solution: F(x*) = x* is sought
# objfn = underlying objective function which is minimized at x*
# method = 1, 2, or 3, indicating the type of steplength to be used
#
#
control.default <- list(K = 1, method=3, minimize=TRUE, square=TRUE, step.min0=1, step.max0=1, mstep=4, kr=1, objfn.inc=1,
tol=1.e-07, maxiter=1500, trace=FALSE, intermed=FALSE)
namc <- names(control)
if (!all(namc %in% names(control.default)))
stop("unknown names in control: ", namc[!(namc %in% names(control.default))])
ctrl <- modifyList(control.default, control)
if (ctrl$K > 1 & !(ctrl$method %in% c("rre", "mpe")) ) ctrl$method <- "rre"
if (ctrl$K == 1 & !(ctrl$method %in% c(1,2,3)) ) ctrl$method <- 3
if (!missing(objfn) ) {
if (ctrl$K == 1) sqobj <- squarem1(par, fixptfn, objfn, ... , control=ctrl)
else if (ctrl$K > 1 | ctrl$method %in% c("rre", "mpe")) sqobj <- cyclem1(par, fixptfn, objfn, ... , control=ctrl)
} else {
if (ctrl$K == 1) sqobj <- squarem2(par, fixptfn, ... , control=ctrl)
else if (ctrl$K > 1 | ctrl$method %in% c("rre", "mpe")) sqobj <- cyclem2(par, fixptfn, ... , control=ctrl)
}
return(sqobj)
}
#######################################################################
# Partially-monotone, globally-convergent SQUAREM acceleration
# Can accelerate EM, MM, and other fixed-point iterations
# Here we use an adaptive maximum step-size strategy
# This strategy is effective and represents a good trade-off between speed and stability
# See Varadhan and Roland (Scandinavian J Statistics 2008)
#
# Date: August 22, 2007
# Adaptive step-max strategy modified: May 2008
# Last modified: November 2008
# Last modified: October 2009
# Last modified: Feb 23 2010
#######################################################################
squarem1 <- function(par, fixptfn, objfn, ... , control=list()) {
# par = starting value of parameter vector
# fixptfn = fixed-point iteration F(x)
# for which the solution: F(x*) = x* is sought
# objfn = underlying objective function which is minimized at x*
#
#
control.default <- list(K=1, square=TRUE, minimize=TRUE, method=3, step.min0=1, step.max0=1, mstep=4, kr=1, objfn.inc=1,
tol=1.e-07, maxiter=1500, trace=FALSE, intermed=FALSE)
namc <- names(control)
if (!all(namc %in% names(control.default)))
stop("unknown names in control: ", namc[!(namc %in% names(control.default))])
ctrl <- modifyList(control.default, control)
#####
# method = 1, 2, or 3, indicating the type of steplength to be used
# A key parameter is `step.min0'. This can be either positive or negative if the eigenvalues of the Jacobian of fixed-point mapping are all positive;
# We set it to +1 for contraction mappings such as EM and MM algorithms
# Must be negative, e.g. -1, if the fixed-point mapping can have negative eiegnvalues
#####
# parameter "objfn.inc" dictates the amount of non-monotonicity in the objective function
# setting objfn.inc=0 would enforce monotonicity, whereas objfn.inc=Inf would be a non-monotonic scheme
# The defalut objfn.inc=1 would enforce monotonicity far from solution, but allows for non-monotonicity closer to solution
#
method <- ctrl$method
maxiter <- ctrl$maxiter
minimize <- ctrl$minimize
tol <- ctrl$tol
step.min <- ctrl$step.min0
step.max0 <- ctrl$step.max0
step.max <- ctrl$step.max0
mstep <- ctrl$mstep
objfn.inc <- ctrl$objfn.inc
trace <- ctrl$trace
intermed <- ctrl$intermed
if (trace) cat("Squarem-1 \n ")
if (missing(objfn)) stop("\n squarem2 should be used if objective function is not available \n\n")
iter <- 1
objval <- rep(NA,1)
p <- par
lold <- objfn(p, ...)
leval <- 1
if (trace) cat("Objective fn: ", lold, "\n")
feval <- 0
conv <- TRUE
p.inter <- c(p, lold)
while (feval < maxiter) {
extrap <- TRUE
p1 <- try(fixptfn(p, ...),silent=TRUE)
feval <- feval + 1
if (inherits(p1, "try-error") | any(is.nan(unlist(p1)))) stop("Error in function evaluation")
q1 <- p1 - p
sr2 <- crossprod(q1)
if (sqrt(sr2) < tol) {p <- p1; break}
p2 <- try(fixptfn(p1, ...),silent=TRUE)
feval <- feval + 1
if (inherits(p2, "try-error") | any(is.nan(unlist(p2)))) stop("Error in function evaluation")
q2 <- p2 - p1
sq2 <- sqrt(crossprod(q2))
if (sq2 < tol) {p <- p2; break}
sv2 <- crossprod(q2-q1)
srv <- crossprod(q1, q2-q1)
alpha <- switch(method, -srv/sv2, -sr2/srv, sqrt(sr2/sv2))
alpha <- max(step.min, min(step.max, alpha))
p.new <- p + 2*alpha*q1 + alpha^2*(q2-q1)
if (abs(alpha - 1) > 0.01 ) {
p.new <- try(fixptfn(p.new , ...),silent=TRUE) # stabilization step
feval <- feval + 1
}
if (inherits(p.new, "try-error") | any(is.nan(p.new))) {
p.new <- p2
lnew <- try(objfn(p2, ...), silent=TRUE)
leval <- leval + 1
if (alpha == step.max) step.max <- max(step.max0, step.max/mstep)
alpha <- 1
extrap <- FALSE
} else {
if (is.finite(objfn.inc)) {
lnew <- try(objfn(p.new, ...), silent=TRUE)
leval <- leval + 1
} else lnew <- lold
non.monotone <- if (minimize) lnew > lold + objfn.inc else lnew < lold - objfn.inc
if (inherits(lnew, "try-error") | is.nan(lnew) | non.monotone) {
p.new <- p2
lnew <- try(objfn(p2, ...), silent=TRUE)
leval <- leval + 1
if (alpha==step.max) step.max <- max(step.max0, step.max/mstep)
alpha <- 1
extrap <- FALSE
}
}
if (alpha == step.max) step.max <- mstep*step.max
if (step.min < 0 & alpha == step.min) step.min <- mstep*step.min
p <- p.new
if (!is.nan(lnew)) lold <- lnew
if (trace) cat("Objective fn: ", lnew, " Extrapolation: ", extrap, " Steplength: ", alpha, "\n")
if(intermed) p.inter <- rbind(p.inter, c(p, lnew))
iter <- iter+1
}
if (feval >= maxiter) conv <- FALSE
if (is.infinite(objfn.inc)) {
lold <- objfn(p, ...)
leval <- leval + 1
}
rownames(p.inter) <- NULL
if (!intermed) return(list(par=p, value.objfn=lold, iter= iter, fpevals=feval, objfevals = leval, convergence=conv))
else return(list(par=p, value.objfn=lold, iter= iter, fpevals=feval, objfevals = leval, convergence=conv, p.intermed=p.inter))
}
###################################################################
squarem2 <- function(par, fixptfn, ... , control=list() ) {
# par = starting value of parameter vector
# fixptfn = fixed-point iteration F(x)
# for which the solution: F(x*) = x* is sought
# method = 1, 2, or 3, indicating the type of steplength to be used
#
# See Varadhan and Roland (Scandinavian J Statistics 2008)
#
control.default <- list(K=1, square=TRUE, minimize=TRUE, method=3, step.min0=1, step.max0=1, mstep=4, kr=1, objfn.inc=1,
tol=1.e-07, maxiter=1500, trace=FALSE, intermed=FALSE)
namc <- names(control)
if (!all(namc %in% names(control.default)))
stop("unknown names in control: ", namc[!(namc %in% names(control.default))])
ctrl <- modifyList(control.default, control)
#####
# method = 1, 2, or 3, indicating the type of steplength to be used
# A key parameter is `step.min0'. This can be either positive or negative if the eigenvalues of the Jacobian of fixed-point mapping are all positive;
# We set it to +1 for contraction mappings such as EM and MM algorithms
# Must be negative, e.g. -1, if the fixed-point mapping can have negative eiegnvalues
#####
# Parameter "kr" dictates the amount of non-monotonicity in the fixed-point residual
# setting kr=0 would try to enforce monotonicity; kres=Inf would be a non-monotonic scheme
#
method <- ctrl$method
maxiter <- ctrl$maxiter
tol <- ctrl$tol
kr <- ctrl$kr
step.min <- ctrl$step.min0
step.max0 <- ctrl$step.max0
step.max <- ctrl$step.max0
mstep <- ctrl$mstep
trace <- ctrl$trace
if (trace) cat("Squarem-2 \n")
intermed <- ctrl$intermed
iter <- 1
feval <- 0
kount <- 0
conv <- TRUE
p.inter <- rep(NA, length(par)+1)
while (feval < maxiter) {
extrap <- TRUE
p1 <- try(fixptfn(par, ...),silent=TRUE)
feval <- feval + 1
if (inherits(p1, "try-error") | any(is.nan(unlist(p1)))) break
q1 <- p1 - par
sr2 <- crossprod(q1)
p.inter <- rbind(p.inter, c(par, sqrt(sr2)))
if (sqrt(sr2) < tol) {par <- p1; break}
p2 <- try(fixptfn(p1, ...),silent=TRUE)
feval <- feval + 1
if (inherits(p2, "try-error") | any(is.nan(unlist(p2)))) break
q2 <- p2 - p1
sq2 <- sqrt(crossprod(q2))
res <- sq2
if (sq2 < tol) {par <- p2; break}
sv2 <- crossprod(q2-q1)
srv <- crossprod(q1, q2-q1)
alpha <- switch(method, -srv/sv2, -sr2/srv, sqrt(sr2/sv2))
alpha <- max(step.min, min(step.max, alpha))
p.new <- par + 2*alpha*q1 + alpha^2*(q2-q1)
if (abs(alpha - 1) > 0.01 ) {
ptmp <- try(fixptfn(p.new, ...),silent=TRUE) # stabilization step
feval <- feval + 1
if (inherits(ptmp, "try-error") | any(is.nan(unlist(ptmp))) ) {
p.new <- p2
if (alpha == step.max) step.max <- max(step.max0, step.max/mstep)
alpha <- 1
extrap <- FALSE
} else {
res <- sqrt(crossprod(ptmp - p.new))
parnorm <- sqrt(crossprod(p2)/length(p2))
kres <- kr * (1 + parnorm) + sq2
p.new <- if (res <= kres) ptmp else p2
if (res > kres) {
if (alpha == step.max) step.max <- max(step.max0, step.max/mstep)
alpha <- 1
extrap <- FALSE
}
}
}
if (alpha == step.max) step.max <- mstep*step.max
if (step.min < 0 & alpha == step.min) step.min <- mstep*step.min
if (trace) cat("Residual: ", res, " Extrapolation: ", extrap, " Steplength: ", alpha, "\n")
par <- p.new
iter <- iter+1
}
p.inter <- p.inter[-1, ]
if (feval >= maxiter) conv <- FALSE
if (intermed) return(list(par=par, value.objfn=NA, iter = iter, fpevals=feval, objfevals = 0,
convergence=conv, p.intermed=p.inter))
else return(list(par=par, value.objfn=NA, iter = iter, fpevals=feval, objfevals = 0, convergence=conv))
}
#######################################################################
# Squared-Cycled vector extrapolation for accelerating vector sequences
# See Roland, Varadhan and Frangakis (Applie Numerical Mathematics 2007)
# Author: Ravi Varadhan, Johns Hopkins University
# Latest modofications:: Feb 23, 2010
#
###########################
cyclem1 <- function(par, fixptfn, objfn, control=list(), ...) {
#
# par = starting value of parameter vector
# fixptfn = fixed-point iteration F(x)
# for which the solution: F(x*) = x* is sought
# objfn = underlying objective function which is minimized at x*
#
control.default <- list(K=2, method="rre", minimize=TRUE, square=TRUE, step.min0=1, step.max0=1, mstep=4, kr=1, objfn.inc=1,
tol=1.e-07, maxiter=1500, trace=FALSE, intermed=FALSE)
namc <- names(control)
if (!all(namc %in% names(control.default)))
stop("unknown names in control: ", namc[!(namc %in% names(control.default))])
ctrl <- modifyList(control.default, control)
#
# method = reduced-rank or minimal-polynomial extrapolation
# K = order of extrapolation scheme; K=2,3,4 are typical choices
# square = a logical variable indicating whether or not "squaring" is used
method <- ctrl$method
K <- ctrl$K
square <- ctrl$square
tol <- ctrl$tol
maxiter <- ctrl$maxiter
objfn.inc <- ctrl$objfn.inc
trace <- ctrl$trace
if (trace) cat("Cyclem-1 \n")
if (missing(objfn)) stop("\n cyclem2 should be used if objective function is not available \n\n")
if ( K > length(par)) {
cat("K is too large. Decrease it. \n")
return()
}
iter <- 0
lold <- objfn(par, ...)
leval <- 1
ptmp <- fixptfn(par,...)
res <- sqrt(c(crossprod(ptmp - par)))
par <- ptmp
feval <- 1
conv <- TRUE
while (feval < maxiter & res > tol) {
iter <- iter+1
extrap <- TRUE
p <- matrix(NA, nrow=length(par), ncol=K+2)
U <- matrix(NA, nrow=length(par), ncol=K+1)
p[,1] <- par
for (i in 2:(K+2)){
p[,i] <- try(fixptfn(p[,i-1], ...), silent=TRUE)
if (inherits(p[,i], "try-error") | any(is.nan(unlist(p[,i]))) | any(!is.numeric(p[,i]))) stop("Error in function evaluation \n")
U[,i-1] <- p[, i] - p[, i-1]
}
feval <- feval + K + 1
sqKp1 <- sqrt(c(crossprod(U[, K+1])))
if (method == "rre"){
coef <- try(solve(qr(crossprod(U), LAPACK=TRUE, tol=1.e-14), rep(1,K+1)), silent=TRUE)
if (inherits(coef, "try-error") | any(is.nan(coef)) ) {
extrap <- FALSE
} else {
if (abs(sum(coef)) < 1.e-07) extrap <- FALSE
coef <- coef/sum(coef)
}
}
if (method == "mpe"){
coef <- try(solve(qr(U[,-(K+1)], LAPACK=TRUE, tol=1.e-14), -U[,K+1]), silent=TRUE)
if (inherits(coef, "try-error") | any(is.nan(coef))) {
extrap <- FALSE
} else {
coef <- c(coef, 1)
if (abs(sum(coef)) < 1.e-07) extrap <- FALSE
coef <- coef/sum(coef)
}
}
if (!extrap) {
par <- p[, ncol(p)]
res <- sqKp1
iter <- iter + 1
if (trace) cat(" objective fn: ", lold, "extrapolation: ", extrap, "\n")
next
}
if (!square) pnew <- c(p[,-(K+2)] %*% coef)
if (square) {
pnew <- rep(0, length(par))
if (K > 1) {
for (i in 1:(K-1)) {
p <- try(cbind(p, fixptfn(p[,i+K+1],...)), silent=TRUE)
if (inherits(p, "try-error")| any(is.nan(unlist(p))) | any(!is.numeric(p))) stop("Error in function evaluation \n")
}
feval <- feval + K - 1
}
for (i in 0:K){
for (j in 0:K){
pnew <- pnew + coef[i+1]* coef[j+1] * p[,i+j+1]
} }
sqKp1 <- sqrt(c(crossprod(p[, 2*K] - p[, 2*K-1])))
}
# Perform one stabilization step between cycles
if (extrap) {
ptmp <- try(fixptfn(pnew, ...), silent=TRUE)
res <- sqrt(c(crossprod(ptmp - pnew)))
feval <- feval + 1
}
if (inherits(ptmp, "try-error") | any(is.nan(unlist(ptmp))) ) {
pnew <- p[, ncol(p)]
res <- sqKp1
extrap <- FALSE
} else pnew <- ptmp
lnew <- try(objfn(pnew, ...), silent=TRUE)
leval <- leval + 1
non.monotone <- if(ctrl$minimize) lnew > lold+objfn.inc else lnew < lold-objfn.inc
if (inherits(lnew, "try-error") | is.nan(lnew) | non.monotone) {
pnew <- p[, ncol(p)]
res <- sqKp1
extrap <- FALSE
if (trace) cat(" objective fn: ", lold, "extrapolation: ", extrap, "\n")
} else lold <- lnew
par <- pnew
if (trace) cat(" objective fn: ", lold, "extrapolation: ", extrap, "\n")
} # Main loop complete
if (feval >= maxiter) conv <- FALSE
return(list(par=par, value.objfn=lold, iter = iter, fpevals=feval, objfevals = leval, convergence=conv))
}
###################################################################
# Squared-Cycled vector extrapolation for accelerating vector sequences
# See Roland, Varadhan and Frangakis (Applie Numerical Mathematics 2007)
# Author: Ravi Varadhan, Johns Hopkins University
# Latest modofications:: Feb 23, 2010
#
###########################
cyclem2 <- function(par, fixptfn, control=list(), ...) {
# par = starting value of parameter vector
# fixptfn = fixed-point iteration F(x)
# for which the solution: F(x*) = x* is sought
# method = reduced-rank or minimal-polynomial extrapolation
#
control.default <- list(K=2, method="rre", minimize=TRUE, square=TRUE, step.min0=1, step.max0=1, mstep=4, kr=1, objfn.inc=1,
tol=1.e-07, maxiter=1500, trace=FALSE, intermed=FALSE)
namc <- names(control)
if (!all(namc %in% names(control.default)))
stop("unknown names in control: ", namc[!(namc %in% names(control.default))])
ctrl <- modifyList(control.default, control)
#
# method = reduced-rank ("rre") or minimal-polynomial ("mpe") extrapolation
# K = order of extrapolation scheme; K=2,3,4 are typical choices
# square = a logical variable indicating whether or not "squaring" is used
method <- ctrl$method
K <- ctrl$K
square <- ctrl$square
tol <- ctrl$tol
maxiter <- ctrl$maxiter
kr <- ctrl$kr
trace <- ctrl$trace
if (trace) cat("Cyclem-2 \n")
if ( K > length(par)) {
cat("K is too large. Decrease it. \n")
return()
}
iter <- 1
ptmp <- fixptfn(par,...)
res <- sqrt(c(crossprod(ptmp - par)))
par <- ptmp
conv <- TRUE
feval <- 1
while (feval < maxiter & res > tol) {
extrap <- TRUE
p <- matrix(NA, nrow=length(par), ncol=K+2)
U <- matrix(NA, nrow=length(par), ncol=K+1)
p[,1] <- par
for (i in 2:(K+2)){
p[,i] <- try(fixptfn(p[,i-1], ...), silent=TRUE)
if (inherits(p[,i], "try-error") | any(is.nan(unlist(p[,i]))) | any(!is.numeric(p[,i]))) stop("Error in function evaluation")
U[,i-1] <- p[, i] - p[, i-1]
}
feval <- feval + K + 1
sqKp1 <- sqrt(c(crossprod(U[, K+1])))
if (method == "rre"){
coef <- try(solve(qr(crossprod(U), LAPACK=TRUE, tol=1.e-14), rep(1,K+1)), silent=TRUE)
if (inherits(coef, "try-error") | any(is.nan(coef))) {
extrap <- FALSE
} else {
if (abs(sum(coef)) < 1.e-07) extrap <- FALSE
coef <- coef/sum(coef)
}
}
if (method == "mpe"){
coef <- try(solve(qr(U[,-(K+1)], LAPACK=TRUE, tol=1.e-14), -U[,K+1]), silent=TRUE)
if (inherits(coef, "try-error") | any(is.nan(coef))) {
extrap <- FALSE
} else {
coef <- c(coef, 1)
if (abs(sum(coef)) < 1.e-07) extrap <- FALSE
coef <- coef/sum(coef)
}
}
if (!extrap) {
par <- p[, ncol(p)]
res <- sqKp1
iter <- iter + 1
if (trace) cat(" residual: ", res, "extrapolation: ", extrap, "\n")
next
}
if (!square) pnew <- c(p[,-(K+2)] %*% coef)
if (square) {
pnew <- rep(0, length(par))
if (K > 1) {
for (i in 1:(K-1)) {
p <- try(cbind(p, fixptfn(p[,i+K+1],...)), silent=TRUE)
if (inherits(p, "try-error") | any(is.nan(unlist(p))) | any(!is.numeric(p))) stop("Error in function evaluation")
}
feval <- feval + K - 1
}
for (i in 0:K){
for (j in 0:K){
pnew <- pnew + coef[i+1]* coef[j+1] * p[,i+j+1]
} }
sqKp1 <- sqrt(c(crossprod(p[, 2*K] - p[, 2*K-1])))
}
# Perform one stabilization step between cycles
ptemp <- try(fixptfn(pnew, ...), silent=TRUE)
feval <- feval + 1
if(inherits(ptemp, "try-error") | any(is.nan(ptemp)) ){
par <- p[, ncol(p)]
res <- sqKp1
extrap <- FALSE
} else {
res <- sqrt(c(crossprod(ptemp - pnew)))
parnorm <- as.numeric(sqrt(crossprod(p[, ncol(p)])/length(par)))
if (res <= (kr * (1 + parnorm) + sqKp1)) {
par <- ptemp
extrap <- TRUE
} else {
par <- p[, ncol(p)]
res <- sqKp1
extrap <- FALSE
}
}
iter <- iter+1
if (trace) cat(" residual: ", res, "extrapolation: ", extrap, "\n")
} # Main loop complete
if (feval >= maxiter) conv <- FALSE
return(list(par=par, value.objfn=NA, iter = iter, fpevals=feval, objfevals = 0, convergence=conv))
}
###################################################################
# The EM algorithm
#
fpiter <- function(par, fixptfn, objfn=NULL, control=list( ), ...){
control.default <- list(tol=1.e-07, maxiter=5000, trace=FALSE, intermed=FALSE)
namc <- names(control)
if (!all(namc %in% names(control.default)))
stop("unknown names in control: ", namc[!(namc %in% names(control.default))])
ctrl <- modifyList(control.default, control)
#
# method = reduced-rank ("rre") or minimal-polynomial ("mpe") extrapolation
# K = order of extrapolation scheme; K=2,3,4 are typical choices
# square = a logical variable indicating whether or not "squaring" is used
tol <- ctrl$tol
maxiter <- ctrl$maxiter
trace <- ctrl$trace
intermed <- ctrl$intermed
if (trace) cat("fpiter \n")
iter <- 1
resid <- rep(NA,1)
objeval <- 0
conv <- FALSE
if(intermed) p.inter <- rep(NA, 1+length(par))
while (iter < maxiter) {
p.new <- fixptfn(par, ...)
res <- sqrt(crossprod(p.new - par))
if(intermed) p.inter <- rbind(p.inter, c(par, res))
if ( res < tol) {conv <- TRUE; break}
if (trace) {
if (!is.null(objfn)) {cat("Iter: ", iter, "Objective fn: ",objfn(par, ...), "\n"); objeval <- objeval + 1}
else cat("Iter: ", iter, "Residual: ",res, "\n")
}
par <- p.new
iter <- iter+1
}
if(intermed) p.inter <- p.inter[-1,]
loglik.best <- if (!is.null(objfn)) objfn(par, ...) else NA
if(!intermed) return(list(par=par, value.objfn=loglik.best, fpevals=iter, objfevals = objeval,
convergence=conv))
else return(list(par=par, value.objfn=loglik.best, fpevals=iter, objfevals = objeval,
p.intermed=p.inter, convergence=conv))
}
###################################################################
Try the SQUAREM 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.