| iSQUAREM {lazy.irtx} | R Documentation |
This function must be generated using generate_iSQUAREM function.
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)
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. |
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
|
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 |
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
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.
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.
# 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 )