diff --git a/DESCRIPTION b/DESCRIPTION index f509677..e007590 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: SQUAREM -Version: 2020.5 -Date: 2020-10-21 +Version: 2021.1 +Date: 2021-01-12 Title: Squared Extrapolation Methods for Accelerating EM-Like Monotone Algorithms Description: Algorithms for accelerating the convergence of slow, @@ -19,5 +19,5 @@ Maintainer: Ravi Varadhan URL: https://coah.jhu.edu/people/Faculty_personal_Pages/Varadhan.html Repository: CRAN NeedsCompilation: no -Packaged: 2020-10-21 20:02:12 UTC; rvaradhan -Date/Publication: 2020-10-21 21:00:02 UTC +Packaged: 2021-01-12 23:59:02 UTC; rvaradhan +Date/Publication: 2021-01-13 06:40:10 UTC diff --git a/MD5 b/MD5 index bea7815..08e9322 100644 --- a/MD5 +++ b/MD5 @@ -1,7 +1,7 @@ -8bb2006ce2485548e5d1b9c39d88678d *DESCRIPTION +f650729f3c5ac0de565ad6d1874de336 *DESCRIPTION a62319b6c4374fd15963266db928a1fc *NAMESPACE -fbff7916cf82707eab9bd5350b9503c2 *NEWS -4ead5ae1cb6525899c0cb75b1e0e7339 *R/squarem.R +481ab651d51595187da4c07ffac78abe *NEWS +1d24cb6e023eb43df0368d8cbc49354f *R/squarem.R 2ce900a9bc6751f7beb84b79942bbac3 *build/vignette.rds 4eb525977b2cd840efcd384cada97dff *demo/00Index 99297f5664d508013f345afcecf8d90c *demo/factoranalysis.R @@ -13,8 +13,10 @@ ac216a9a91116475a6b28b933f34e2b4 *inst/CITATION 4df21965b20b54fd14215afe50bff3e4 *inst/doc/SQUAREM.R 8baa991f5faaaabb04be226d2aac8db0 *inst/doc/SQUAREM.Rnw 87ec2678b267ecd99d38567c6937cd25 *inst/doc/SQUAREM.pdf -028e8abac54c45e3d66ff7d5a54e5ca7 *man/fpiter.Rd +5e1af6c79774307501a7e4e864ad1f5c *man/fpiter.Rd dc86f820449b08ca2a684572e599dd55 *man/internal.Rd -1734baa4947c57dc975c3456d325435c *man/squarem.Rd +f894e4e0dc269079cfacf0fe1703b423 *man/squarem.Rd 1f6e68e9322bd524f4db846a970670f0 *tests/Hasselblad1969.R 8baa991f5faaaabb04be226d2aac8db0 *vignettes/SQUAREM.Rnw +50fb1d088187a510f7156630ee9d0a3e *vignettes/SQUAREM.pdf +29124e75e51c284de8359c8c9cec6bbd *vignettes/SQUAREM.tex diff --git a/NEWS b/NEWS index 2a82b31..dd419b2 100644 --- a/NEWS +++ b/NEWS @@ -29,9 +29,15 @@ Log of changes to SQUAREM - Modified the code to output updated parameters from EM algorithm immediately before breaking for convergence o August 19, 2020 - - Fixed a couple of lines in squarem2() to output updated parameters from EM algorithm immediately before breaking for convergence - thanks to Dan Schaid (Mayo) for pointing this out + - Fixed a couple of lines in squarem2() to output updated parameters from EM algorithm + immediately before breaking for convergence - thanks to Dan Schaid (Mayo) for pointing this out o October 21, 2020 - Added the option to either "minimize"" or "maximize" the objective function - Control parameter `minimize' can be set to TRUE (for minimization) or FALSE (for maximization) + o January 12, 2021 + - Added the option to output the residuals when objective funcion is NOT given + - also added this option to `fpiter' + - In order to use this option, the user has to set `intermed=TRUE' + - this option is not available for higher-order squarem algorithms \ No newline at end of file diff --git a/R/squarem.R b/R/squarem.R index 1670a46..0a0b9e4 100644 --- a/R/squarem.R +++ b/R/squarem.R @@ -175,7 +175,7 @@ while (feval < maxiter) { 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.inter=p.inter)) + 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() ) { @@ -214,11 +214,13 @@ ctrl <- modifyList(control.default, control) 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 @@ -227,6 +229,7 @@ while (feval < maxiter) { 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) @@ -273,9 +276,12 @@ while (feval < maxiter) { par <- p.new iter <- iter+1 } - if (feval >= maxiter) conv <- FALSE +p.inter <- p.inter[-1, ] +if (feval >= maxiter) conv <- FALSE -return(list(par=par, value.objfn=NA, iter = iter, fpevals=feval, objfevals = 0, convergence=conv)) +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)) } ####################################################################### @@ -564,7 +570,7 @@ return(list(par=par, value.objfn=NA, iter = iter, fpevals=feval, objfevals = 0, # fpiter <- function(par, fixptfn, objfn=NULL, control=list( ), ...){ -control.default <- list(tol=1.e-07, maxiter=5000, trace=FALSE) +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))]) @@ -578,7 +584,8 @@ ctrl <- modifyList(control.default, control) tol <- ctrl$tol maxiter <- ctrl$maxiter trace <- ctrl$trace - + intermed <- ctrl$intermed + if (trace) cat("fpiter \n") iter <- 1 @@ -586,10 +593,13 @@ 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} @@ -600,10 +610,13 @@ if (trace) { par <- p.new iter <- iter+1 } - +if(intermed) p.inter <- p.inter[-1,] loglik.best <- if (!is.null(objfn)) objfn(par, ...) else NA -return(list(par=par, value.objfn=loglik.best, fpevals=iter, objfevals = objeval, convergence=conv)) +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)) } ################################################################### diff --git a/man/fpiter.Rd b/man/fpiter.Rd index 4e66d0d..82030a0 100644 --- a/man/fpiter.Rd +++ b/man/fpiter.Rd @@ -50,21 +50,25 @@ \code{control} is list of control parameters for the algorithm. \describe{ -\code{control = list(tol = 1.e-07, maxiter = 1500, trace = FALSE)} +\code{control = list(tol = 1.e-07, maxiter = 1500, trace = FALSE, intermed=FALSE)} - \code{tol}{ A small, positive scalar that determines when iterations + \code{tol}{ - a small, positive scalar that determines when iterations should be terminated. Iteration is terminated when \eqn{||x_k - F(x_k)|| \leq tol}{abs(x[k] - F(x[k]) <= tol}. Default is \code{1.e-07}.} - \code{maxiter}{ An integer denoting the maximum limit on the number of + \code{maxiter}{ - an integer denoting the maximum limit on the number of evaluations of \code{fixptfn}, \eqn{F}{F}. Default is \code{1500}.} - \code{trace}{ A logical variable denoting whether some of the intermediate + \code{trace}{ - a logical variable denoting whether some of the intermediate results of iterations should be displayed to the user. Default is \code{FALSE}.} + + \code{intermed}}{ - a logical variable denoting whether the intermediate + results of iterations should be returned. If set to \code{TRUE}, the function + will return a matrix where each row corresponds to parameters at each iteration, + along with the corresponding fixed-point residual value. Default is \code{FALSE}.} } -} \seealso{ \code{\link{squarem}} @@ -101,7 +105,7 @@ poissmix.dat <- data.frame(death=0:9, freq=c(162,267,271,185,111,61,27,8,3,1)) y <- poissmix.dat$freq tol <- 1.e-08 -# Use a preset seed so the example is reproducable. +# Use a preset seed so the example is reproducible. require("setRNG") old.seed <- setRNG(list(kind="Mersenne-Twister", normal.kind="Inversion", seed=54321)) diff --git a/man/squarem.Rd b/man/squarem.Rd index 9bab986..634187e 100644 --- a/man/squarem.Rd +++ b/man/squarem.Rd @@ -47,7 +47,7 @@ Fixed-Point Iterations} indicates successful convergence, whereas \code{1} denotes failure to converge.} -\item{p.inter}{A matrix where each row corresponds to parameters at each iteration, +\item{p.intermed}{A matrix where each row corresponds to parameters at each iteration, along with the corresponding log-likelihood value. This object is returned only when the control parameter \code{intermed} is set to \code{TRUE}. It is not returned when \code{objfn} is not specified.} @@ -129,10 +129,9 @@ Default values of \code{control} are: results of iterations should be displayed to the user. Default is \code{FALSE}.} \item{\code{intermed}}{A logical variable denoting whether the intermediate - results of iterations should be returned. If set to \code{TRUE}, the function + results of iterat ions should be returned. If set to \code{TRUE}, the function will return a matrix where each row corresponds to parameters at each iteration, - along with the corresponding log-likelihood value. This option is inactive when - \code{objfn} is not specified. + along with the corresponding log-likelihood value. When the code{objfn} is not specified it will return the fixed-point residual instead of the objective function values. Default is \code{FALSE}.} } } diff --git a/vignettes/SQUAREM.pdf b/vignettes/SQUAREM.pdf new file mode 100644 index 0000000..624d19a Binary files /dev/null and b/vignettes/SQUAREM.pdf differ diff --git a/vignettes/SQUAREM.tex b/vignettes/SQUAREM.tex new file mode 100644 index 0000000..ffba3d5 --- /dev/null +++ b/vignettes/SQUAREM.tex @@ -0,0 +1,55 @@ +\documentclass[12pt]{article} + +\usepackage[margin=1in]{geometry} +\usepackage{amsmath, amssymb, amsfonts} +%\usepackage{natbib} +\usepackage{graphicx} +\usepackage{color} %% red, green, and blue (for screen display) and cyan, magenta, and yellow +\definecolor{Navy}{rgb}{0,0,0.8} +\usepackage{hyperref} +\hypersetup{colorlinks=true, urlcolor={Navy}, linkcolor={Navy}, citecolor={Navy}} + +\parskip 7.2pt + +\newcommand{\compresslist}{% +%\setlength{\itemsep}{1pt}% +\setlength{\itemsep}{0pt}% +\setlength{\parskip}{0pt}% +\setlength{\parsep}{0pt}% +} + +\newcommand{\pb}{\mathbb{P}} +\newcommand{\E}{\mathbb{E}} +\newcommand{\V}{\mathbb{V}} +\newcommand{\C}{\mathbb{C}} +\newcommand{\bea}{\begin{align*}} +\newcommand{\eea}{\end{align*}} +\newcommand{\beq}{\begin{equation}} +\newcommand{\eeq}{\end{equation}} +\newcommand{\be}{\begin{enumerate}} +\newcommand{\ee}{\end{enumerate}} +\newcommand{\bi}{\begin{itemize}} +\newcommand{\ei}{\end{itemize}} +\renewcommand{\baselinestretch}{1} + +\title{\texttt{SQUAREM}: Accelerating the Convergence of EM, MM and Other Fixed-Point Algorithms} +\author{Ravi Varadhan} +\date{} +\usepackage{Sweave} +\begin{document} + +%\VignetteIndexEntry{SQUAREM Tutorial} +%\VignetteDepends{setRNG} +%\VignetteKeywords{EM algorithm, fixed-point iteration, acceleration, extrapolation} +%\VignettePackage{SQUAREM} + +\maketitle + +\section{Overview of SQUAREM} +''SQUAREM'' is a package intended for accelerating slowly-convergent contraction mappings. It can be used for accelerating the convergence of slow, linearly convergent contraction mappings such as the EM (expectation-maximization) algorithm, MM (majorize and minimize) algorithm, and other nonlinear fixed-point iterations such as the power method for finding the dominant eigenvector. It uses a novel approach callled squared extrapolation method (SQUAREM) that was proposed in Varadhan and Roland (Scandinavian Journal of Statistics, 35: 335-353), and also in Roland, Vardhan, and Frangakis (Numerical Algorithms, 44: 159-172). + +The functions in this package are made available with: + +\begin{Schunk} +\begin{Sinput} +> library("SQUAREM") \ No newline at end of file