iSQUAREM {lazy.accel}R Documentation

iSQUAREM A function to be used to update parameters using inline SQUAREM:
Inline Squared Extrapolation Methods for accelerating fixed-point iterations

Description

This function must be generated using generate_iSQUAREM function.

Usage

iSQUAREM(param, ..., enforce_constraints = NULL, badness_of_fit = NULL,
  SQUAREM = 3, minalpha = -999, maxalpha = -1, bof_value = NULL,
  always = 0, reset1 = 1, reset2 = 2, print = 0, debug = 0)

Arguments

param

A parameter vector to be updated

...

Additional parameters to the functions enforce_constraints and badness_of_fit

enforce_constraints

A function to enfoce the constraints, if any, on the parameter vector.
This will be called as enforce_constraints(param, ...) from iSQUAREM function.

badness_of_fit

A function to evaluate the badness of fit of the current parameter vector. This will be called as badness_of_fit(param, ...) from iSQUAREM function.

SQUAREM

= 1 or 3

minalpha

The minimul value of alpha parameter of iSQUAREM

maxalpha

<= 1 The maximum value of alpha parameter of iSQUAREM

bof_value

current value of the badness of fit prior to iSQUAREM
Usually, it is not necessary to specify this.

always

= 1 to update parameters regardless of the badness of fit.

reset1

= 0 or 1 counter reset value after update

reset2

= 1 or 2 counter reset value after failure when always=0

print

= 1 to print the update process

debug

= 1 to print the generated function

Details

This is a dummy version of iSQUAREM function only to display its help.
The real one must be generated using generate_iSQUAREM as a closure in the function which uses the accelaration process.

How to use iSQUAREM:
Try always=1 with maxalpha=-1 or less first.
Changing to reset1=0 and reset2=1 or increasing nSQUAREM may help.
If it seems not working, use always=0 with maxalpha=1 or less.
Changing to reset1=0 and reset2=1 may help.
If all of the above fail, be patient and use SQUAREM=0.

When always=1 the badness of fit will not be evaluated in this function.
However, note that, even if always=0, the badness of fit evaluated in the calling function may not improve if iSQUAREM update is accepted with badness1 < badness where badness1 is the badness of fit evaluated at new update.
This happens because badness in this function is the value evaluated at the last SQUAREM update and not the one evaluated with the latest param value in the calling function.

The following is an example of a function which uses iSQUAREM.


my_algorithm <- function( data, param, etc ){

 # function definitions
 const <- function( param, additional_param_list ){
  # enforce constraints on param
  # All the variables except those listed as formal arguments
  # reffer to the ones defined in the environmend in which
  # this function is defined: namely environment(my_algorithm).
  return( param )
 } # end of const
 fit <- function( param, additional_param_list ){
  # calculate the badness of fit criterion of param
  # All the variables except those listed as formal arguments
  # reffer to the ones defined in the environmend in which
  # this function is defined: namely environment(my_algorithm).
  # This function can return a list whose first element is named as critval.
  return( crit )
 } # end of fit

 # initial value of parameter vector
 param <- initial_value
 converged <- 0

 ################ new block of code added 1 #########################
 if( SQUAREM ){
  # generate iSQUREM function as a closure
  iSQUAREM <- generate_iSQUAREM( param, print=1 )
  }
 ####################################################################

 for( iter in 1:maxiter ){

  # usual parameter update process satisfying constraints.
  param <- update_param( param )

  ################ new block of code added 2 #########################
  if( SQUAREM > 0  &  nSQUAREM <= iter ){
   # acceleration
   temp <- iSQUAREM( param, additional_param_list
                 , enforce_constraints=const, badness_of_fit=crit
                 , SQUAREM=SQUAREM )
   param <- temp$param
   critval <- temp$critval
  } # end of iSQUAREM
  ####################################################################

  # check convergence in terms of critval or changes of param value.
  if( converged ) break

 } # end of iterations

 return( param )

} # end of my_algorithm

Value

A list of
param The updated parameter value with constraints if any, or the input param as it is.
critval The value of badness_of_fit function when always=0 or NA.

References

R Varadhan and C Roland (2008), Simple and globally convergent numerical schemes for acceler- ating the convergence of any EM algorithm, Scandinavian Journal of Statistics, 35:335-353.
C Roland, R Varadhan, and CE Frangakis (2007), Squared polynomial extrapolation methods with cycling: an application to the positron emission tomography problem, Numerical Algorithms 44:159-172.

Examples

# Example No. 1:  Power method for eigen values/vectors.

power <- function( A, maxiter=500, epsd=1e-6
                 , SQUAREM=0, nSQUAREM=1, always=0 ){
 # power method for the largest eigen value and the associated vector

 n=nrow(A)

 # initial value
 eve <- rep(1,n)
 evep <- eve

 #################### block 1 for iSQUAREM #############################
 if( SQUAREM ){
  # generate iSQUAREM function
  iSQUAREM <- generate_iSQUAREM( eve, print=1 )
  }
 #################### block 1 for iSQUAREM #############################

 # main iteration loop
 for( iter in 1:maxiter ){

  eve <- A%*%eve
  eva <- sqrt(sum(eve*eve))
  eve <- eve / eva

  #################### block 2 for iSQUAREM #############################
  # use iSQUAREM function
  if( SQUAREM  &  iter >= nSQUAREM ){
   res <- iSQUAREM( eve, always=always, print=1 )
   eve <- res$param
  }
  #################### block 2 for iSQUAREM #############################

  # check convergence
  if( max(abs(abs(eve)-abs(evep))) <= epsd ) break

  evep <- eve

 } # end of main iteration

 ev <- eigen(A)
 evec_true <- (ev$vector)[,1]
 eva_true <- (ev$values)[1]
 Print( iter, SQUAREM, always )
 Print( eve, evec_true, eva, eva_true, fmt="8.5" )

} # end of power


# generate data
set.seed(1701)
n <- 50
A <- matrix(rnorm(n*n),n)
A <- t(A)%*%A
evec <- eigen(A)$vector

# no acceleration
power( A, SQUAREM=0 )
# with acceleration
power( A, SQUAREM=3, always=1 )




# Example No. 2:  Lower rank approximation of a matrix

pca_by_als <- function( Y, ndim=2, maxiter=100, eps=1e-7
                        , SQUAREM=0, nSQUAREM=1, maxalpha=-1, always=1
                        , reset1=1, reset2=2 ){
 # pca by als
 # Shin-ichi Mayekawa
 # 20160103
 #

 # This function finds F and A matrices of rank ndim such that they minimize
 #   sum( (Y-F%*%t(A))^2 )
 #

 get_rmse <- function( Y, F, A ){
  # root mean square errors as a badness of fit criterion
  rmse <- sqrt( sum( (Y-F%*%t(A))^2 )/nrow(Y)/ncol(Y) )
  return( rmse )
 } # end of get_rmse

 get_crit <- function ( param ){
  # badness of fit function for iSQUAREM
  #
  # Note that those variables except for "param" reffer to the ones
  # defined in the environmnet in which this function is defined.
  #

  # recover matrices from vector
  F <- matrix(param[1:(n*ndim)],n)
  A <- matrix(param[-(1:(n*ndim))],nvar)
  # evaluate rmse
  rmse <- sqrt( sum( (Y-F%*%t(A))^2 )/nrow(Y)/ncol(Y) )
  return( rmse )
 } # end of get_crit


 # matrix info
 n <- nrow(Y); nvar <- ncol(Y)

 # initial values
 F <- Y[,1:ndim,drop=0]
 A <- diag(ndim)
 A <- rbind(A,matrix(0,nvar-ndim,ndim))

 # previous values
 rmsep <- 99999
 Fp <- F; Ap <- A


 #################### block 1 for iSQUAREM #############################
 if( SQUAREM ){
  # parameter vector
  param <- c(c(F),c(A))
  # generate iSQUAREM function
  iSQUAREM <- generate_iSQUAREM( param, print=1 )
 }
 #################### block 1 for iSQUAREM #############################


 # main iteration loop
 for( iter in 1:maxiter ){

  # full conditional LSE
  A <- t( solve(t(F)%*%F)%*%t(F)%*%Y )
  F <- t( solve(t(A)%*%A)%*%t(A)%*%t(Y) )

  #################### block 2 for iSQUAREM #############################
  # use iSQUAREM function
  if( SQUAREM  &  iter >= nSQUAREM ){
   # parameter vector
   param <- c(c(F),c(A))
   param <- iSQUAREM( param, badness_of_fit=get_crit
                      , maxalpha=maxalpha, always=always
                      , reset1=reset1, reset2=reset2, print=1 )$param
   F <- matrix(param[1:(n*ndim)],n)
   A <- matrix(param[-(1:(n*ndim))],nvar)
  }
  #################### block 2 for iSQUAREM #############################


  # evaluate fit
  rmse <- get_rmse( Y, F, A )
  rmseimp <- rmsep-rmse
  maxad <- max(abs(F-Fp),abs(A-Ap))

  if( maxad <= eps ) break

  # print info
  Print( iter, rmse, rmseimp, maxad, fmt="3.0  8.6" )

  # next iteration
  rmsep <- rmse
  Fp <- F; Ap <- A

 } # end of main iteration loop

 cat( "Iteration terminated with ", iter, " iterations.\n" )


} # end of pca_by_als


# generate data
seed <- 1701
set.seed(seed)
n <- 100; nvar <- 9
Y <- matrix(rnorm(n*nvar),n)

# no acceleration
res <- pca_by_als( Y, 2, SQUAREM=0 )
# with acceleration 1
res <- pca_by_als( Y, 2, SQUAREM=3, always=1 )
# with acceleration 2
res <- pca_by_als( Y, 2, SQUAREM=3, always=0 )

[Package lazy.accel version 0.1.2 Index]